<div dir="ltr"><div class="gmail_default" style="font-family:verdana,sans-serif">Hi Fabian,</div><div class="gmail_default" style="font-family:verdana,sans-serif"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif">Looking at the FreeBSD code, I see there are other changes that differ from newlib's version starting from:</div><div class="gmail_default" style="font-family:verdana,sans-serif"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif"> else if(hx < 0x39000000) { /* when |x|<2**-14 */ <=== in newlib this is 2** - 23<br> if(huge+x>one) return one+x;/* trigger inexact */<br> }<br> else k = 0;<br><br> /* x is now in primary range */<br> t = x*x;<br> if(k >= -125)<br> SET_FLOAT_WORD(twopk,((u_int32_t)(0x7f+k))<<23);<br> else<br> SET_FLOAT_WORD(twopk,((u_int32_t)(0x7f+(k+100)))<<23);<br> c = x - t*(P1+t*P2);<br> if(k==0) return one-((x*c)/(c-(float)2.0)-x);<br> else y = one-((lo-(x*c)/((float)2.0-c))-hi);<br> if(k >= -125) {<br> if(k==128) return y*2.0F*0x1p127F;<br> return y*twopk;<br> } else {<br> return y*twopk*twom100;<br> }</div><div class="gmail_default" style="font-family:verdana,sans-serif"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif">Would you like to port over those changes as well? Some of them are from 18 years ago and there are a few from 6 years ago (ones with 0x7f).<br></div><div class="gmail_default" style="font-family:verdana,sans-serif"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif">-- Jeff J.<br></div></div><br><div class="gmail_quote gmail_quote_container"><div dir="ltr" class="gmail_attr">On Tue, Jan 28, 2025 at 5:52 AM Fabian Schriever <<a href="mailto:fabian.schriever@gtd-gmbh.de">fabian.schriever@gtd-gmbh.de</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">that was used doesn't work normally here, since we want to be able to<br>
multiply `hi' by the exponent of x _exactly_, and the exponent of x has<br>
more than 7 significant bits for most denormal x's, so the multiplication<br>
was not always exact despite a cloned comment claiming that it was. (The<br>
comment is correct in the double precision case -- with the normal 33+53<br>
bit decomposition the exponent can have 20 significant bits and the extra<br>
bit for denormals is only the 11th.)<br>
<br>
Fixing this had little or no effect for denormals (I think because<br>
more precision is inherently lost for denormals than is lost by roundoff<br>
errors in the multiplication).<br>
<br>
The fix is to reduce the precision of the decomposition to 16+24 bits.<br>
Due to 2 bugs in the old deomposition and numerical accidents, reducing<br>
the precision actually increased the precision of hi+lo. The old hi+lo<br>
had about 39 bits instead of at least 41 like it should have had.<br>
There were off-by-1-bit errors in each of hi and lo, apparently due<br>
to mistranslation from the double precision hi and lo. The correct<br>
16 bit hi happens to give about 19 bits of precision, so the correct<br>
hi+lo gives about 43 bits instead of at least 40. The end result is<br>
that expf() is now perfectly rounded (to nearest) except in 52561 cases<br>
instead of except in 67027 cases, and the maximum error is 0.5013 ulps<br>
instead of 0.5023 ulps.<br>
<br>
Commit and message are from FreeBSD.<br>
Reference: <a href="https://github.com/freebsd/freebsd-src/commit/06d803185574a0ba37434b1e98b322e54cc32a4e" rel="noreferrer" target="_blank">https://github.com/freebsd/freebsd-src/commit/06d803185574a0ba37434b1e98b322e54cc32a4e</a><br>
Original Author: Bruce Evans<br>
<br>
We have taken the liberty to also add the f suffix to all float<br>
constants to ensure they are not doubles.<br>
---<br>
newlib/libm/math/ef_exp.c | 34 +++++++++++++++++-----------------<br>
1 file changed, 17 insertions(+), 17 deletions(-)<br>
<br>
diff --git a/newlib/libm/math/ef_exp.c b/newlib/libm/math/ef_exp.c<br>
index d6e25bfcd..85e1cf9d3 100644<br>
--- a/newlib/libm/math/ef_exp.c<br>
+++ b/newlib/libm/math/ef_exp.c<br>
@@ -26,20 +26,20 @@ static const float<br>
#else<br>
static float<br>
#endif<br>
-one = 1.0,<br>
-halF[2] = {0.5,-0.5,},<br>
-huge = 1.0e+30,<br>
-twom100 = 7.8886090522e-31, /* 2**-100=0x0d800000 */<br>
-ln2HI[2] ={ 6.9313812256e-01, /* 0x3f317180 */<br>
- -6.9313812256e-01,}, /* 0xbf317180 */<br>
-ln2LO[2] ={ 9.0580006145e-06, /* 0x3717f7d1 */<br>
- -9.0580006145e-06,}, /* 0xb717f7d1 */<br>
-invln2 = 1.4426950216e+00, /* 0x3fb8aa3b */<br>
-P1 = 1.6666667163e-01, /* 0x3e2aaaab */<br>
-P2 = -2.7777778450e-03, /* 0xbb360b61 */<br>
-P3 = 6.6137559770e-05, /* 0x388ab355 */<br>
-P4 = -1.6533901999e-06, /* 0xb5ddea0e */<br>
-P5 = 4.1381369442e-08; /* 0x3331bb4c */<br>
+one = 1.0f,<br>
+halF[2] = {0.5f,-0.5f,},<br>
+huge = 1.0e+30f,<br>
+twom100 = 7.8886090522e-31f, /* 2**-100=0x0d800000 */<br>
+ln2HI[2] ={ 6.9314575195e-01f, /* 0x3f317200 */<br>
+ -6.9314575195e-01f,}, /* 0xbf317200 */<br>
+ln2LO[2] ={ 1.4286067653e-06f, /* 0x35bfbe8e */<br>
+ -1.4286067653e-06f,}, /* 0xb5bfbe8e */<br>
+invln2 = 1.4426950216e+00f, /* 0x3fb8aa3b */<br>
+P1 = 1.6666667163e-01f, /* 0x3e2aaaab */<br>
+P2 = -2.7777778450e-03f, /* 0xbb360b61 */<br>
+P3 = 6.6137559770e-05f, /* 0x388ab355 */<br>
+P4 = -1.6533901999e-06f, /* 0xb5ddea0e */<br>
+P5 = 4.1381369442e-08f; /* 0x3331bb4c */<br>
<br>
#ifdef __STDC__<br>
float __ieee754_expf(float x) /* default IEEE double exp */<br>
@@ -60,7 +60,7 @@ P5 = 4.1381369442e-08; /* 0x3331bb4c */<br>
if(FLT_UWORD_IS_NAN(hx))<br>
return x+x; /* NaN */<br>
if(FLT_UWORD_IS_INFINITE(hx))<br>
- return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */<br>
+ return (xsb==0)? x:0.0f; /* exp(+-inf)={inf,0} */<br>
if(sx > FLT_UWORD_LOG_MAX)<br>
return __math_oflowf(0); /* overflow */<br>
if(sx < 0 && hx > FLT_UWORD_LOG_MIN)<br>
@@ -85,8 +85,8 @@ P5 = 4.1381369442e-08; /* 0x3331bb4c */<br>
/* x is now in primary range */<br>
t = x*x;<br>
c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));<br>
- if(k==0) return one-((x*c)/(c-(float)2.0)-x); <br>
- else y = one-((lo-(x*c)/((float)2.0-c))-hi);<br>
+ if(k==0) return one-((x*c)/(c-2.0f)-x); <br>
+ else y = one-((lo-(x*c)/(2.0f-c))-hi);<br>
if(k >= -125) {<br>
__uint32_t hy;<br>
GET_FLOAT_WORD(hy,y);<br>
-- <br>
2.33.0.windows.1<br>
<br>
<br>
</blockquote></div>