From 159c7655d06f04aa56a57385d633699c4c63d72c Mon Sep 17 00:00:00 2001 From: Szabolcs Nagy Date: Sat, 17 Nov 2012 23:39:39 +0100 Subject: [PATCH] math: cleanup exp2.c exp2f.c and exp2l.c * old code relied on sign extension on right shift * exp2l ld64 wrapper was wrong * use scalbn instead of bithacks --- src/math/exp2.c | 53 ++++++++++++++++++------------------------- src/math/exp2f.c | 30 ++++++++++++------------ src/math/exp2l.c | 59 ++++++++++++++++-------------------------------- 3 files changed, 56 insertions(+), 86 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); } diff --git a/src/math/exp2f.c b/src/math/exp2f.c index 55c22eac..279f32de 100644 --- a/src/math/exp2f.c +++ b/src/math/exp2f.c @@ -27,19 +27,15 @@ #include "libm.h" -#define TBLBITS 4 -#define TBLSIZE (1 << TBLBITS) +#define TBLSIZE 16 static const float -huge = 0x1p100f, redux = 0x1.8p23f / TBLSIZE, P1 = 0x1.62e430p-1f, P2 = 0x1.ebfbe0p-3f, P3 = 0x1.c6b348p-5f, P4 = 0x1.3b2c9cp-7f; -static const volatile float twom100 = 0x1p-100f; - static const double exp2ft[TBLSIZE] = { 0x1.6a09e667f3bcdp-1, 0x1.7a11473eb0187p-1, @@ -89,23 +85,25 @@ float exp2f(float x) { double tv, twopk, u, z; float t; - uint32_t hx, ix, i0; - int32_t k; + uint32_t hx, ix, i0, k; /* Filter out exceptional cases. */ GET_FLOAT_WORD(hx, x); ix = hx & 0x7fffffff; if (ix >= 0x43000000) { /* |x| >= 128 */ if (ix >= 0x7f800000) { - if ((ix & 0x7fffff) != 0 || (hx & 0x80000000) == 0) - return x + x; /* x is NaN or +Inf */ - else - return 0.0; /* x is -Inf */ + if (hx == 0xff800000) /* -inf */ + return 0; + return x; + } + if (x >= 128) { + STRICT_ASSIGN(float, x, x * 0x1p127); + return x; + } + if (x <= -150) { + STRICT_ASSIGN(float, x, 0x1p-100*0x1p-100); + return x; } - if (x >= 0x1.0p7f) - return huge * huge; /* overflow */ - if (x <= -0x1.2cp7f) - return twom100 * twom100; /* underflow */ } else if (ix <= 0x33000000) { /* |x| <= 0x1p-25 */ return 1.0f + x; } @@ -114,7 +112,7 @@ float exp2f(float x) STRICT_ASSIGN(float, t, x + redux); GET_FLOAT_WORD(i0, t); i0 += TBLSIZE / 2; - k = (i0 >> TBLBITS) << 20; + k = (i0 / TBLSIZE) << 20; i0 &= TBLSIZE - 1; t -= redux; z = x - t; diff --git a/src/math/exp2l.c b/src/math/exp2l.c index b587308f..145c20fe 100644 --- a/src/math/exp2l.c +++ b/src/math/exp2l.c @@ -30,7 +30,7 @@ #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 long double exp2l(long double x) { - return exp2l(x); + return exp2(x); } #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 @@ -40,10 +40,6 @@ long double exp2l(long double x) #define BIAS (LDBL_MAX_EXP - 1) #define EXPMASK (BIAS + LDBL_MAX_EXP) -static const long double huge = 0x1p10000L; -/* XXX Prevent gcc from erroneously constant folding this. */ -static const volatile long double twom10000 = 0x1p-10000L; - static const double redux = 0x1.8p63 / TBLSIZE, P1 = 0x1.62e42fefa39efp-1, @@ -208,27 +204,28 @@ static const double tbl[TBLSIZE * 2] = { long double exp2l(long double x) { union IEEEl2bits u, v; - long double r, twopk, twopkp10000, z; + long double r, z; uint32_t hx, ix, i0; - int k; + union {uint32_t u; int32_t i;} k; /* Filter out exceptional cases. */ u.e = x; hx = u.xbits.expsign; ix = hx & EXPMASK; if (ix >= BIAS + 14) { /* |x| >= 16384 or x is NaN */ - if (ix == BIAS + LDBL_MAX_EXP) { - if (u.xbits.man != 1ULL << 63 || (hx & 0x8000) == 0) - return x + x; /* x is +Inf or NaN */ - return 0.0; /* x is -Inf */ + if (ix == EXPMASK) { + if (u.xbits.man == 1ULL << 63 && hx == 0xffff) /* -inf */ + return 0; + return x; + } + if (x >= 16384) { + x *= 0x1p16383L; + return x; } - if (x >= 16384) - return huge * huge; /* overflow */ if (x <= -16446) - return twom10000 * twom10000; /* underflow */ - } else if (ix <= BIAS - 66) { /* |x| < 0x1p-66 */ - return 1.0 + x; - } + return 0x1p-10000L*0x1p-10000L; + } else if (ix < BIAS - 64) /* |x| < 0x1p-64 */ + return 1 + x; /* * Reduce x, computing z, i0, and k. The low bits of x + redux @@ -240,38 +237,22 @@ long double exp2l(long double x) * Then the low-order word of x + redux is 0x000abc12, * We split this into k = 0xabc and i0 = 0x12 (adjusted to * index into the table), then we compute z = 0x0.003456p0. - * - * XXX If the exponent is negative, the computation of k depends on - * '>>' doing sign extension. */ u.e = x + redux; i0 = u.bits.manl + TBLSIZE / 2; - k = (int)i0 >> TBLBITS; - i0 = (i0 & (TBLSIZE - 1)) << 1; + k.u = i0 / TBLSIZE * TBLSIZE; + k.i /= TBLSIZE; + i0 %= TBLSIZE; 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]; - long double t_lo = tbl[i0 + 1]; + long double t_hi = tbl[2*i0]; + long double t_lo = tbl[2*i0 + 1]; /* XXX This gives > 1 ulp errors outside of FE_TONEAREST mode */ r = t_lo + (t_hi + t_lo) * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * (P5 + z * P6))))) + t_hi; - /* Scale by 2**k. */ - if (k >= LDBL_MIN_EXP) { - if (k == LDBL_MAX_EXP) - return r * 2.0 * 0x1p16383L; - return r * twopk; - } - return r * twopkp10000 * twom10000; + return scalbnl(r, k.i); } #endif -- 2.20.1