diff options
-rw-r--r-- | lib/msun/ld128/s_exp2l.c | 26 | ||||
-rw-r--r-- | lib/msun/ld80/s_exp2l.c | 25 | ||||
-rw-r--r-- | lib/msun/src/e_exp.c | 17 | ||||
-rw-r--r-- | lib/msun/src/e_expf.c | 17 | ||||
-rw-r--r-- | lib/msun/src/s_exp2.c | 18 | ||||
-rw-r--r-- | lib/msun/src/s_exp2f.c | 9 |
6 files changed, 57 insertions, 55 deletions
diff --git a/lib/msun/ld128/s_exp2l.c b/lib/msun/ld128/s_exp2l.c index 28a6935..a94b892 100644 --- a/lib/msun/ld128/s_exp2l.c +++ b/lib/msun/ld128/s_exp2l.c @@ -356,8 +356,8 @@ static const float eps[TBLSIZE] = { long double exp2l(long double x) { - union IEEEl2bits u; - long double r, t, z; + union IEEEl2bits u, v; + long double r, t, twopk, twopkp10000, z; uint32_t hx, ix, i0; int k; @@ -403,6 +403,15 @@ exp2l(long double x) i0 = i0 & (TBLSIZE - 1); u.e -= redux; z = x - u.e; + v.xbits.manh = 0; + v.xbits.manl = 0; + if (k >= LDBL_MIN_EXP) { + v.xbits.expsign = LDBL_MAX_EXP - 1 + k; + twopk = v.e; + } else { + v.xbits.expsign = LDBL_MAX_EXP - 1 + k + 10000; + twopkp10000 = v.e; + } /* Compute r = exp2(y) = exp2t[i0] * p(z - eps[i]). */ t = tbl[i0]; /* exp2t[i0] */ @@ -412,15 +421,10 @@ exp2l(long double x) /* Scale by 2**k. */ if(k >= LDBL_MIN_EXP) { - if (k != 0) { - u.e = r; - u.xbits.expsign += k; - r = u.e; - } - return (r); + if (k == LDBL_MAX_EXP) + return (r * 2.0 * 0x1p16383L); + return (r * twopk); } else { - u.e = r; - u.xbits.expsign += k + 10000; - return (u.e * twom10000); + return (r * twopkp10000 * twom10000); } } diff --git a/lib/msun/ld80/s_exp2l.c b/lib/msun/ld80/s_exp2l.c index 10c7da7..8a0dbe2 100644 --- a/lib/msun/ld80/s_exp2l.c +++ b/lib/msun/ld80/s_exp2l.c @@ -214,8 +214,8 @@ static const double tbl[TBLSIZE * 2] = { long double exp2l(long double x) { - union IEEEl2bits u; - long double r, z; + union IEEEl2bits u, v; + long double r, twopk, twopkp10000, z; uint32_t hx, ix, i0; int k; @@ -267,6 +267,14 @@ exp2l(long double x) i0 = (i0 & (TBLSIZE - 1)) << 1; u.e -= redux; z = x - u.e; + v.xbits.man = 1ULL << 63; + if (k >= LDBL_MIN_EXP) { + v.xbits.expsign = LDBL_MAX_EXP - 1 + k; + twopk = v.e; + } else { + v.xbits.expsign = LDBL_MAX_EXP - 1 + k + 10000; + twopkp10000 = v.e; + } /* Compute r = exp2l(y) = exp2lt[i0] * p(z). */ long double t_hi = tbl[i0]; @@ -277,15 +285,10 @@ exp2l(long double x) /* Scale by 2**k. */ if (k >= LDBL_MIN_EXP) { - if (k != 0) { - u.e = r; - u.xbits.expsign += k; - return (u.e); - } - return (r); + if (k == LDBL_MAX_EXP) + return (r * 2.0 * 0x1p16383L); + return (r * twopk); } else { - u.e = r; - u.xbits.expsign += k + 10000; - return (u.e * twom10000); + return (r * twopkp10000 * twom10000); } } diff --git a/lib/msun/src/e_exp.c b/lib/msun/src/e_exp.c index 60b8b2a..0809c66 100644 --- a/lib/msun/src/e_exp.c +++ b/lib/msun/src/e_exp.c @@ -102,7 +102,7 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ double __ieee754_exp(double x) /* default IEEE double exp */ { - double y,hi=0.0,lo=0.0,c,t; + double y,hi=0.0,lo=0.0,c,t,twopk; int32_t k=0,xsb; u_int32_t hx; @@ -142,18 +142,17 @@ __ieee754_exp(double x) /* default IEEE double exp */ /* x is now in primary range */ t = x*x; + if(k >= -1021) + INSERT_WORDS(twopk,0x3ff00000+(k<<20), 0); + else + INSERT_WORDS(twopk,0x3ff00000+((k+1000)<<20), 0); c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); if(k==0) return one-((x*c)/(c-2.0)-x); else y = one-((lo-(x*c)/(2.0-c))-hi); if(k >= -1021) { - u_int32_t hy; - GET_HIGH_WORD(hy,y); - SET_HIGH_WORD(y,hy+(k<<20)); /* add k to y's exponent */ - return y; + if (k==1024) return y*2.0*0x1p1023; + return y*twopk; } else { - u_int32_t hy; - GET_HIGH_WORD(hy,y); - SET_HIGH_WORD(y,hy+((k+1000)<<20)); /* add k to y's exponent */ - return y*twom1000; + return y*twopk*twom1000; } } diff --git a/lib/msun/src/e_expf.c b/lib/msun/src/e_expf.c index 949fa54..45eda01 100644 --- a/lib/msun/src/e_expf.c +++ b/lib/msun/src/e_expf.c @@ -43,7 +43,7 @@ static volatile float twom100 = 7.8886090522e-31; /* 2**-100=0x0d800000 */ float __ieee754_expf(float x) /* default IEEE double exp */ { - float y,hi=0.0,lo=0.0,c,t; + float y,hi=0.0,lo=0.0,c,t,twopk; int32_t k=0,xsb; u_int32_t hx; @@ -80,18 +80,17 @@ __ieee754_expf(float x) /* default IEEE double exp */ /* x is now in primary range */ t = x*x; + if(k >= -125) + SET_FLOAT_WORD(twopk,0x3f800000+(k<<23)); + else + SET_FLOAT_WORD(twopk,0x3f800000+((k+100)<<23)); c = x - t*(P1+t*P2); 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); - SET_FLOAT_WORD(y,hy+(k<<23)); /* add k to y's exponent */ - return y; + if(k==128) return y*2.0F*0x1p127F; + return y*twopk; } else { - u_int32_t hy; - GET_FLOAT_WORD(hy,y); - SET_FLOAT_WORD(y,hy+((k+100)<<23)); /* add k to y's exponent */ - return y*twom100; + return y*twopk*twom100; } } diff --git a/lib/msun/src/s_exp2.c b/lib/msun/src/s_exp2.c index 56554f1..63b8997 100644 --- a/lib/msun/src/s_exp2.c +++ b/lib/msun/src/s_exp2.c @@ -340,7 +340,7 @@ static const double tbl[TBLSIZE * 2] = { double exp2(double x) { - double r, t, z; + double r, t, twopk, twopkp1000, z; uint32_t hx, hr, ix, lx, i0; int k; @@ -375,19 +375,19 @@ exp2(double x) /* Compute r = exp2(y) = exp2t[i0] * p(z - eps[i]). */ t = tbl[i0]; /* exp2t[i0] */ z -= tbl[i0 + 1]; /* eps[i0] */ + if (k >= -1021 << 20) + INSERT_WORDS(twopk, 0x3ff00000 + k, 0); + else + INSERT_WORDS(twopkp1000, 0x3ff00000 + k + (1000 << 20), 0); r = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * P5)))); /* Scale by 2**(k>>20). */ if(k >= -1021 << 20) { - if (k != 0) { - GET_HIGH_WORD(hr, r); - SET_HIGH_WORD(r, hr + k); - } - return (r); + if (k == 1024 << 20) + return (r * 2.0 * 0x1p1023); + return (r * twopk); } else { - GET_HIGH_WORD(hr, r); - SET_HIGH_WORD(r, hr + (k + (1000 << 20))); - return (r * twom1000); + return (r * twopkp1000 * twom1000); } } diff --git a/lib/msun/src/s_exp2f.c b/lib/msun/src/s_exp2f.c index ffe58a2..d0f6b2f 100644 --- a/lib/msun/src/s_exp2f.c +++ b/lib/msun/src/s_exp2f.c @@ -92,7 +92,7 @@ static const double exp2ft[TBLSIZE] = { float exp2f(float x) { - double tv; + double tv, twopk; float t, z; uint32_t hx, htv, ix, i0; int32_t k; @@ -129,9 +129,6 @@ exp2f(float x) tv = tv + tv * (z * (P1 + z * (P2 + z * (P3 + z * P4)))); /* Scale by 2**(k>>20). */ - if (k != 0) { - GET_HIGH_WORD(htv, tv); - SET_HIGH_WORD(tv, htv + k); - } - return (tv); + INSERT_WORDS(twopk, 0x3ff00000 + k, 0); + return (tv * twopk); } |