d/freebsd-skq
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Revert rev 1.8, which causes small (e.g. 2 ulp) errors for some

Browse Source 
inputs.  The trouble with replacing two floats with a double is that
the latter has 6 extra bits of precision, which actually hurts
accuracy in many cases.  All of the constants are optimal when float
arithmetic is used, and would need to be recomputed to do this right.

Noticed by:	bde (ucbtest)
This commit is contained in:
David Schultz 2005-02-24 06:32:13 +00:00
parent eaadc155a8
commit aa28340df9
1 changed files with 13 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
lib/msun/src/e_expf.c
View File

@ -27,24 +27,27 @@ huge = 1.0e+30,
twom100 = 7.8886090522e-31, /* 2**-100=0x0d800000 */
o_threshold= 8.8721679688e+01, /* 0x42b17180 */
u_threshold= -1.0397208405e+02, /* 0xc2cff1b5 */
ln2HI[2] ={ 6.9313812256e-01, /* 0x3f317180 */
-6.9313812256e-01,}, /* 0xbf317180 */
ln2LO[2] ={ 9.0580006145e-06, /* 0x3717f7d1 */
-9.0580006145e-06,}, /* 0xb717f7d1 */
invln2 = 1.4426950216e+00, /* 0x3fb8aa3b */
P1 = 1.6666667163e-01, /* 0x3e2aaaab */
P2 = -2.7777778450e-03, /* 0xbb360b61 */
P3 = 6.6137559770e-05, /* 0x388ab355 */
P4 = -1.6533901999e-06, /* 0xb5ddea0e */
P5 = 4.1381369442e-08; /* 0x3331bb4c */
double ln2[2] = { 6.93147180369123816490e-01, -6.93147180369123816490e-01 };
float
__ieee754_expf(float x) /* IEEE float exp */
__ieee754_expf(float x) /* default IEEE double exp */
{
float y,c,t;
float y,hi=0.0,lo=0.0,c,t;
int32_t k=0,xsb;
u_int32_t hx;
GET_FLOAT_WORD(hx,x);
xsb = (hx>>31)&1; /* sign bit of x */
hx &= 0x7fffffff; /* |x| */
hx &= 0x7fffffff; /* high word of |x| */
/* filter out non-finite argument */
if(hx >= 0x42b17218) { /* if |x|>=88.721... */
@ -59,12 +62,14 @@ __ieee754_expf(float x) /* IEEE float exp */
/* argument reduction */
if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */
if(hx < 0x3F851592) { /* and |x| < 1.5 ln2 */
x = x-ln2[xsb]; k = 1-xsb-xsb;
hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
} else {
k = invln2*x+halF[xsb];
t = k;
x = x - t*ln2[0];
hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */
lo = t*ln2LO[0];
}
x = hi - lo;
}
else if(hx < 0x31800000) { /* when |x|<2**-28 */
if(huge+x>one) return one+x;/* trigger inexact */
@ -74,8 +79,8 @@ __ieee754_expf(float x) /* IEEE float exp */
/* x is now in primary range */
t = x*x;
c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
y = one-(((double)x*c)/(c-2.0)-x);
if(k==0) return y;
if(k==0) return one-((x*c)/(c-(float)2.0)-x);
else y = one-((lo-(x*c)/((float)2.0-c))-hi);
if(k >= -125) {
u_int32_t hy;
GET_FLOAT_WORD(hy,y);