diff options
Diffstat (limited to 'src/math/exp2.c')
| -rw-r--r-- | src/math/exp2.c | 53 |
1 files changed, 22 insertions, 31 deletions
diff --git a/src/math/exp2.c b/src/math/exp2.c index 08c21a66..8e252280 100644 --- a/src/math/exp2.c +++ b/src/math/exp2.c @@ -27,11 +27,9 @@ #include "libm.h" -#define TBLBITS 8 -#define TBLSIZE (1 << TBLBITS) +#define TBLSIZE 256 static const double -huge = 0x1p1000, redux = 0x1.8p52 / TBLSIZE, P1 = 0x1.62e42fefa39efp-1, P2 = 0x1.ebfbdff82c575p-3, @@ -39,8 +37,6 @@ P3 = 0x1.c6b08d704a0a6p-5, P4 = 0x1.3b2ab88f70400p-7, P5 = 0x1.5d88003875c74p-10; -static const volatile double twom1000 = 0x1p-1000; - static const double tbl[TBLSIZE * 2] = { /* exp2(z + eps) eps */ 0x1.6a09e667f3d5dp-1, 0x1.9880p-44, @@ -334,25 +330,28 @@ static const double tbl[TBLSIZE * 2] = { */ double exp2(double x) { - double r, t, twopk, twopkp1000, z; - uint32_t hx, ix, lx, i0; - int k; + double r, t, z; + uint32_t hx, ix, i0; + union {uint32_t u; int32_t i;} k; /* Filter out exceptional cases. */ GET_HIGH_WORD(hx, x); ix = hx & 0x7fffffff; if (ix >= 0x40900000) { /* |x| >= 1024 */ if (ix >= 0x7ff00000) { - GET_LOW_WORD(lx, x); - if (((ix & 0xfffff) | lx) != 0 || (hx & 0x80000000) == 0) - return x + x; /* x is NaN or +Inf */ - else - return 0.0; /* x is -Inf */ + GET_LOW_WORD(ix, x); + if (hx == 0xfff00000 && ix == 0) /* -inf */ + return 0; + return x; + } + if (x >= 1024) { + STRICT_ASSIGN(double, x, x * 0x1p1023); + return x; + } + if (x <= -1075) { + STRICT_ASSIGN(double, x, 0x1p-1000*0x1p-1000); + return x; } - if (x >= 0x1.0p10) - return huge * huge; /* overflow */ - if (x <= -0x1.0ccp10) - return twom1000 * twom1000; /* underflow */ } else if (ix < 0x3c900000) { /* |x| < 0x1p-54 */ return 1.0 + x; } @@ -361,24 +360,16 @@ double exp2(double x) STRICT_ASSIGN(double, t, x + redux); GET_LOW_WORD(i0, t); i0 += TBLSIZE / 2; - k = (i0 >> TBLBITS) << 20; - i0 = (i0 & (TBLSIZE - 1)) << 1; + k.u = i0 / TBLSIZE * TBLSIZE; + k.i /= TBLSIZE; + i0 %= TBLSIZE; t -= redux; z = x - t; /* 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); + t = tbl[2*i0]; /* exp2t[i0] */ + z -= tbl[2*i0 + 1]; /* eps[i0] */ r = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * P5)))); - /* Scale by 2**(k>>20). */ - if (k < -1021 << 20) - return r * twopkp1000 * twom1000; - if (k == 1024 << 20) - return r * 2.0 * 0x1p1023; - return r * twopk; + return scalbn(r, k.i); } |