diff options
Diffstat (limited to 'src/math/exp2.c')
| -rw-r--r-- | src/math/exp2.c | 40 |
1 files changed, 20 insertions, 20 deletions
diff --git a/src/math/exp2.c b/src/math/exp2.c index 8e252280..2e078fb0 100644 --- a/src/math/exp2.c +++ b/src/math/exp2.c @@ -305,7 +305,7 @@ static const double tbl[TBLSIZE * 2] = { * Method: (accurate tables) * * Reduce x: - * x = 2**k + y, for integer k and |y| <= 1/2. + * x = k + y, for integer k and |y| <= 1/2. * Thus we have exp2(x) = 2**k * exp2(y). * * Reduce y: @@ -330,41 +330,41 @@ static const double tbl[TBLSIZE * 2] = { */ double exp2(double x) { - double r, t, z; - uint32_t hx, ix, i0; + double_t r, t, z; + uint32_t ix, i0; + union {double f; uint64_t i;} u = {x}; 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(ix, x); - if (hx == 0xfff00000 && ix == 0) /* -inf */ - return 0; - return x; - } - if (x >= 1024) { + ix = u.i>>32 & 0x7fffffff; + if (ix >= 0x408ff000) { /* |x| >= 1022 or nan */ + if (ix >= 0x40900000 && u.i>>63 == 0) { /* x >= 1024 or nan */ + /* overflow */ STRICT_ASSIGN(double, x, x * 0x1p1023); return x; } - if (x <= -1075) { - STRICT_ASSIGN(double, x, 0x1p-1000*0x1p-1000); - return x; + if (ix >= 0x7ff00000) /* -inf or -nan */ + return -1/x; + if (u.i>>63) { /* x <= -1022 */ + /* underflow */ + if (x <= -1075 || x - 0x1p52 + 0x1p52 != x) + FORCE_EVAL((float)(-0x1p-149/x)); + if (x <= -1075) + return 0; } } else if (ix < 0x3c900000) { /* |x| < 0x1p-54 */ return 1.0 + x; } /* Reduce x, computing z, i0, and k. */ - STRICT_ASSIGN(double, t, x + redux); - GET_LOW_WORD(i0, t); + u.f = x + redux; + i0 = u.i; i0 += TBLSIZE / 2; k.u = i0 / TBLSIZE * TBLSIZE; k.i /= TBLSIZE; i0 %= TBLSIZE; - t -= redux; - z = x - t; + u.f -= redux; + z = x - u.f; /* Compute r = exp2(y) = exp2t[i0] * p(z - eps[i]). */ t = tbl[2*i0]; /* exp2t[i0] */ |