From 39c910fb061114e6aa5c3bf2c94b1d7262d62221 Mon Sep 17 00:00:00 2001 From: Szabolcs Nagy Date: Wed, 4 Sep 2013 07:51:11 +0000 Subject: [PATCH] math: fix underflow in exp*.c and long double handling in exp2l * don't care about inexact flag * use double_t and float_t (faster, smaller, more precise on x86) * exp: underflow when result is zero or subnormal and not -inf * exp2: underflow when result is zero or subnormal and not exact * expm1: underflow when result is zero or subnormal * expl: don't underflow on -inf * exp2: fix incorrect comment * expm1: simplify special case handling and overflow properly * expm1: cleanup final scaling and fix negative left shift ub (twopk) --- src/math/exp.c | 21 +++++++------- src/math/exp2.c | 40 +++++++++++++------------- src/math/exp2f.c | 52 ++++++++++++++++------------------ src/math/exp2l.c | 44 +++++++++++++---------------- src/math/expf.c | 25 ++++++++-------- src/math/expl.c | 2 +- src/math/expm1.c | 72 ++++++++++++++++++----------------------------- src/math/expm1f.c | 65 +++++++++++++++++------------------------- 8 files changed, 139 insertions(+), 182 deletions(-) diff --git a/src/math/exp.c b/src/math/exp.c index 5ec8f8a7..2b182acd 100644 --- a/src/math/exp.c +++ b/src/math/exp.c @@ -80,7 +80,7 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ double exp(double x) { - double hi, lo, c, xx; + double_t hi, lo, c, xx, y; int k, sign; uint32_t hx; @@ -89,20 +89,19 @@ double exp(double x) hx &= 0x7fffffff; /* high word of |x| */ /* special cases */ - if (hx >= 0x40862e42) { /* if |x| >= 709.78... */ + if (hx >= 0x4086232b) { /* if |x| >= 708.39... */ if (isnan(x)) return x; - if (hx == 0x7ff00000 && sign) /* -inf */ - return 0; if (x > 709.782712893383973096) { /* overflow if x!=inf */ STRICT_ASSIGN(double, x, 0x1p1023 * x); return x; } - if (x < -745.13321910194110842) { - /* underflow */ - STRICT_ASSIGN(double, x, 0x1p-1000 * 0x1p-1000); - return x; + if (x < -708.39641853226410622) { + /* underflow if x!=-inf */ + FORCE_EVAL((float)(-0x1p-149/x)); + if (x < -745.13321910194110842) + return 0; } } @@ -128,8 +127,8 @@ double exp(double x) /* x is now in primary range */ xx = x*x; c = x - xx*(P1+xx*(P2+xx*(P3+xx*(P4+xx*P5)))); - x = 1 + (x*c/(2-c) - lo + hi); + y = 1 + (x*c/(2-c) - lo + hi); if (k == 0) - return x; - return scalbn(x, k); + return y; + return scalbn(y, k); } 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] */ diff --git a/src/math/exp2f.c b/src/math/exp2f.c index ea50db4a..91738f04 100644 --- a/src/math/exp2f.c +++ b/src/math/exp2f.c @@ -63,7 +63,7 @@ static const double exp2ft[TBLSIZE] = { * Method: (equally-spaced 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 exp2f(x) = 2**k * exp2(y). * * Reduce y: @@ -83,46 +83,42 @@ static const double exp2ft[TBLSIZE] = { */ float exp2f(float x) { - double tv, twopk, u, z; - float t; - uint32_t hx, ix, i0, k; + double_t t, r, z; + union {float f; uint32_t i;} u = {x}; + union {double f; uint64_t i;} uk; + uint32_t ix, i0, k; /* Filter out exceptional cases. */ - GET_FLOAT_WORD(hx, x); - ix = hx & 0x7fffffff; - if (ix >= 0x43000000) { /* |x| >= 128 */ - if (ix >= 0x7f800000) { - if (hx == 0xff800000) /* -inf */ - return 0; - return x; - } - if (x >= 128) { + ix = u.i & 0x7fffffff; + if (ix > 0x42fc0000) { /* |x| > 126 */ + if (u.i >= 0x43000000 && u.i < 0x80000000) { /* x >= 128 */ STRICT_ASSIGN(float, x, x * 0x1p127f); return x; } - if (x <= -150) { - STRICT_ASSIGN(float, x, 0x1p-100f*0x1p-100f); - return x; + if (u.i >= 0x80000000) { /* x < -126 */ + if (u.i >= 0xc3160000 || (u.i & 0x0000ffff)) + FORCE_EVAL(-0x1p-149f/x); + if (u.i >= 0xc3160000) /* x <= -150 */ + return 0; } } else if (ix <= 0x33000000) { /* |x| <= 0x1p-25 */ return 1.0f + x; } /* Reduce x, computing z, i0, and k. */ - STRICT_ASSIGN(float, t, x + redux); - GET_FLOAT_WORD(i0, t); + u.f = x + redux; + i0 = u.i; i0 += TBLSIZE / 2; - k = (i0 / TBLSIZE) << 20; + k = i0 / TBLSIZE; + uk.i = (uint64_t)(0x3ff + k)<<52; i0 &= TBLSIZE - 1; - t -= redux; - z = x - t; - INSERT_WORDS(twopk, 0x3ff00000 + k, 0); - + u.f -= redux; + z = x - u.f; /* Compute r = exp2(y) = exp2ft[i0] * p(z). */ - tv = exp2ft[i0]; - u = tv * z; - tv = tv + u * (P1 + z * P2) + u * (z * z) * (P3 + z * P4); + r = exp2ft[i0]; + t = r * z; + r = r + t * (P1 + z * P2) + t * (z * z) * (P3 + z * P4); - /* Scale by 2**(k>>20). */ - return tv * twopk; + /* Scale by 2**k */ + return r * uk.f; } diff --git a/src/math/exp2l.c b/src/math/exp2l.c index 145c20fe..8fc4037c 100644 --- a/src/math/exp2l.c +++ b/src/math/exp2l.c @@ -33,13 +33,9 @@ long double exp2l(long double x) return exp2(x); } #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 - #define TBLBITS 7 #define TBLSIZE (1 << TBLBITS) -#define BIAS (LDBL_MAX_EXP - 1) -#define EXPMASK (BIAS + LDBL_MAX_EXP) - static const double redux = 0x1.8p63 / TBLSIZE, P1 = 0x1.62e42fefa39efp-1, @@ -203,29 +199,29 @@ static const double tbl[TBLSIZE * 2] = { */ long double exp2l(long double x) { - union IEEEl2bits u, v; + union ldshape u = {x}; + int e = u.i.se & 0x7fff; long double r, z; - uint32_t hx, ix, i0; + uint32_t i0; 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 == EXPMASK) { - if (u.xbits.man == 1ULL << 63 && hx == 0xffff) /* -inf */ + if (e >= 0x3fff + 13) { /* |x| >= 8192 or x is NaN */ + if (u.i.se >= 0x3fff + 14 && u.i.se >> 15 == 0) + /* overflow */ + return x * 0x1p16383L; + if (e == 0x7fff) /* -inf or -nan */ + return -1/x; + if (x < -16382) { + if (x <= -16446 || x - 0x1p63 + 0x1p63 != x) + /* underflow */ + FORCE_EVAL((float)(-0x1p-149/x)); + if (x <= -16446) return 0; - return x; - } - if (x >= 16384) { - x *= 0x1p16383L; - return x; } - if (x <= -16446) - return 0x1p-10000L*0x1p-10000L; - } else if (ix < BIAS - 64) /* |x| < 0x1p-64 */ + } else if (e < 0x3fff - 64) { return 1 + x; + } /* * Reduce x, computing z, i0, and k. The low bits of x + redux @@ -238,13 +234,13 @@ long double exp2l(long double x) * We split this into k = 0xabc and i0 = 0x12 (adjusted to * index into the table), then we compute z = 0x0.003456p0. */ - u.e = x + redux; - i0 = u.bits.manl + TBLSIZE / 2; + u.f = x + redux; + i0 = u.i.m + TBLSIZE / 2; k.u = i0 / TBLSIZE * TBLSIZE; k.i /= TBLSIZE; i0 %= TBLSIZE; - u.e -= redux; - z = x - u.e; + u.f -= redux; + z = x - u.f; /* Compute r = exp2l(y) = exp2lt[i0] * p(z). */ long double t_hi = tbl[2*i0]; diff --git a/src/math/expf.c b/src/math/expf.c index 8aefc917..5572bbf6 100644 --- a/src/math/expf.c +++ b/src/math/expf.c @@ -29,7 +29,7 @@ P2 = -2.7667332906e-3f; /* -0xb55215.0p-32 */ float expf(float x) { - float hi, lo, c, xx; + float_t hi, lo, c, xx, y; int k, sign; uint32_t hx; @@ -38,20 +38,17 @@ float expf(float x) hx &= 0x7fffffff; /* high word of |x| */ /* special cases */ - if (hx >= 0x42b17218) { /* if |x| >= 88.722839f or NaN */ - if (hx > 0x7f800000) /* NaN */ - return x; - if (!sign) { - /* overflow if x!=inf */ + if (hx >= 0x42aeac50) { /* if |x| >= -87.33655f or NaN */ + if (hx >= 0x42b17218 && !sign) { /* x >= 88.722839f */ + /* overflow */ STRICT_ASSIGN(float, x, x * 0x1p127f); return x; } - if (hx == 0x7f800000) /* -inf */ - return 0; - if (hx >= 0x42cff1b5) { /* x <= -103.972084f */ + if (sign) { /* underflow */ - STRICT_ASSIGN(float, x, 0x1p-100f*0x1p-100f); - return x; + FORCE_EVAL(-0x1p-149f/x); + if (hx >= 0x42cff1b5) /* x <= -103.972084f */ + return 0; } } @@ -77,8 +74,8 @@ float expf(float x) /* x is now in primary range */ xx = x*x; c = x - xx*(P1+xx*P2); - x = 1 + (x*c/(2-c) - lo + hi); + y = 1 + (x*c/(2-c) - lo + hi); if (k == 0) - return x; - return scalbnf(x, k); + return y; + return scalbnf(y, k); } diff --git a/src/math/expl.c b/src/math/expl.c index 50a04297..b62980fa 100644 --- a/src/math/expl.c +++ b/src/math/expl.c @@ -100,7 +100,7 @@ long double expl(long double x) if (x > 11356.5234062941439488L) /* x > ln(2^16384 - 0.5) */ return x * 0x1p16383L; if (x < -11399.4985314888605581L) /* x < ln(2^-16446) */ - return 0x1p-10000L * 0x1p-10000L; + return -0x1p-16445L/x; /* Express e**x = e**f 2**k * = e**(f + k ln(2)) diff --git a/src/math/expm1.c b/src/math/expm1.c index f8f32c46..a7eb2c0b 100644 --- a/src/math/expm1.c +++ b/src/math/expm1.c @@ -31,7 +31,7 @@ * R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r) * = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r)) * = 1 - r^2/60 + r^4/2520 - r^6/100800 + ... - * We use a special Reme algorithm on [0,0.347] to generate + * We use a special Remez algorithm on [0,0.347] to generate * a polynomial of degree 5 in r*r to approximate R1. The * maximum error of this polynomial approximation is bounded * by 2**-61. In other words, @@ -107,8 +107,6 @@ #include "libm.h" static const double -huge = 1.0e+300, -tiny = 1.0e-300, o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ ln2_hi = 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ ln2_lo = 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */ @@ -122,39 +120,27 @@ Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */ double expm1(double x) { - double y,hi,lo,c,t,e,hxs,hfx,r1,twopk; - int32_t k,xsb; - uint32_t hx; - - GET_HIGH_WORD(hx, x); - xsb = hx&0x80000000; /* sign bit of x */ - hx &= 0x7fffffff; /* high word of |x| */ + double_t y,hi,lo,c,t,e,hxs,hfx,r1,twopk; + union {double f; uint64_t i;} u = {x}; + uint32_t hx = u.i>>32 & 0x7fffffff; + int k, sign = u.i>>63; /* filter out huge and non-finite argument */ if (hx >= 0x4043687A) { /* if |x|>=56*ln2 */ - if (hx >= 0x40862E42) { /* if |x|>=709.78... */ - if (hx >= 0x7ff00000) { - uint32_t low; - - GET_LOW_WORD(low, x); - if (((hx&0xfffff)|low) != 0) /* NaN */ - return x+x; - return xsb==0 ? x : -1.0; /* exp(+-inf)={inf,-1} */ - } - if(x > o_threshold) - return huge*huge; /* overflow */ - } - if (xsb != 0) { /* x < -56*ln2, return -1.0 with inexact */ - /* raise inexact */ - if(x+tiny<0.0) - return tiny-1.0; /* return -1 */ + if (isnan(x)) + return x; + if (sign) + return -1; + if (x > o_threshold) { + x *= 0x1p1023; + return x; } } /* argument reduction */ if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ - if (xsb == 0) { + if (!sign) { hi = x - ln2_hi; lo = ln2_lo; k = 1; @@ -164,7 +150,7 @@ double expm1(double x) k = -1; } } else { - k = invln2*x + (xsb==0 ? 0.5 : -0.5); + k = invln2*x + (sign ? -0.5 : 0.5); t = k; hi = x - t*ln2_hi; /* t*ln2_hi is exact here */ lo = t*ln2_lo; @@ -172,9 +158,9 @@ double expm1(double x) STRICT_ASSIGN(double, x, hi - lo); c = (hi-x)-lo; } else if (hx < 0x3c900000) { /* |x| < 2**-54, return x */ - /* raise inexact flags when x != 0 */ - t = huge+x; - return x - (t-(huge+x)); + if (hx < 0x00100000) + FORCE_EVAL((float)x); + return x; } else k = 0; @@ -186,9 +172,9 @@ double expm1(double x) e = hxs*((r1-t)/(6.0 - x*t)); if (k == 0) /* c is 0 */ return x - (x*e-hxs); - INSERT_WORDS(twopk, 0x3ff00000+(k<<20), 0); /* 2^k */ e = x*(e-c) - c; e -= hxs; + /* exp(x) ~ 2^k (x_reduced - e + 1) */ if (k == -1) return 0.5*(x-e) - 0.5; if (k == 1) { @@ -196,24 +182,20 @@ double expm1(double x) return -2.0*(e-(x+0.5)); return 1.0+2.0*(x-e); } - if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */ - y = 1.0 - (e-x); + u.i = (uint64_t)(0x3ff + k)<<52; /* 2^k */ + twopk = u.f; + if (k < 0 || k > 56) { /* suffice to return exp(x)-1 */ + y = x - e + 1.0; if (k == 1024) y = y*2.0*0x1p1023; else y = y*twopk; return y - 1.0; } - t = 1.0; - if (k < 20) { - SET_HIGH_WORD(t, 0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */ - y = t-(e-x); - y = y*twopk; - } else { - SET_HIGH_WORD(t, ((0x3ff-k)<<20)); /* 2^-k */ - y = x-(e+t); - y += 1.0; - y = y*twopk; - } + u.i = (uint64_t)(0x3ff - k)<<52; /* 2^-k */ + if (k < 20) + y = (x-e+(1-u.f))*twopk; + else + y = (x-(e+u.f)+1)*twopk; return y; } diff --git a/src/math/expm1f.c b/src/math/expm1f.c index d9568466..698ab45f 100644 --- a/src/math/expm1f.c +++ b/src/math/expm1f.c @@ -16,8 +16,6 @@ #include "libm.h" static const float -huge = 1.0e+30, -tiny = 1.0e-30, o_threshold = 8.8721679688e+01, /* 0x42b17180 */ ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ @@ -32,35 +30,27 @@ Q2 = 1.5807170421e-3; /* 0xcf3010.0p-33 */ float expm1f(float x) { - float y,hi,lo,c,t,e,hxs,hfx,r1,twopk; - int32_t k,xsb; - uint32_t hx; - - GET_FLOAT_WORD(hx, x); - xsb = hx&0x80000000; /* sign bit of x */ - hx &= 0x7fffffff; /* high word of |x| */ + float_t y,hi,lo,c,t,e,hxs,hfx,r1,twopk; + union {float f; uint32_t i;} u = {x}; + uint32_t hx = u.i & 0x7fffffff; + int k, sign = u.i >> 31; /* filter out huge and non-finite argument */ if (hx >= 0x4195b844) { /* if |x|>=27*ln2 */ - if (hx >= 0x42b17218) { /* if |x|>=88.721... */ - if (hx > 0x7f800000) /* NaN */ - return x+x; - if (hx == 0x7f800000) /* exp(+-inf)={inf,-1} */ - return xsb==0 ? x : -1.0; - if (x > o_threshold) - return huge*huge; /* overflow */ - } - if (xsb != 0) { /* x < -27*ln2 */ - /* raise inexact */ - if (x+tiny < 0.0f) - return tiny-1.0f; /* return -1 */ + if (hx > 0x7f800000) /* NaN */ + return x; + if (sign) + return -1; + if (x > o_threshold) { + x *= 0x1p127f; + return x; } } /* argument reduction */ if (hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ if (hx < 0x3F851592) { /* and |x| < 1.5 ln2 */ - if (xsb == 0) { + if (!sign) { hi = x - ln2_hi; lo = ln2_lo; k = 1; @@ -70,7 +60,7 @@ float expm1f(float x) k = -1; } } else { - k = invln2*x + (xsb==0 ? 0.5f : -0.5f); + k = invln2*x + (sign ? -0.5f : 0.5f); t = k; hi = x - t*ln2_hi; /* t*ln2_hi is exact here */ lo = t*ln2_lo; @@ -78,8 +68,9 @@ float expm1f(float x) STRICT_ASSIGN(float, x, hi - lo); c = (hi-x)-lo; } else if (hx < 0x33000000) { /* when |x|<2**-25, return x */ - t = huge+x; /* return x with inexact flags when x!=0 */ - return x - (t-(huge+x)); + if (hx < 0x00800000) + FORCE_EVAL(x*x); + return x; } else k = 0; @@ -91,9 +82,9 @@ float expm1f(float x) e = hxs*((r1-t)/(6.0f - x*t)); if (k == 0) /* c is 0 */ return x - (x*e-hxs); - SET_FLOAT_WORD(twopk, 0x3f800000+(k<<23)); /* 2^k */ e = x*(e-c) - c; e -= hxs; + /* exp(x) ~ 2^k (x_reduced - e + 1) */ if (k == -1) return 0.5f*(x-e) - 0.5f; if (k == 1) { @@ -101,24 +92,20 @@ float expm1f(float x) return -2.0f*(e-(x+0.5f)); return 1.0f + 2.0f*(x-e); } - if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */ - y = 1.0f - (e - x); + u.i = (0x7f+k)<<23; /* 2^k */ + twopk = u.f; + if (k < 0 || k > 56) { /* suffice to return exp(x)-1 */ + y = x - e + 1.0f; if (k == 128) y = y*2.0f*0x1p127f; else y = y*twopk; return y - 1.0f; } - t = 1.0f; - if (k < 23) { - SET_FLOAT_WORD(t, 0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */ - y = t - (e - x); - y = y*twopk; - } else { - SET_FLOAT_WORD(t, (0x7f-k)<<23); /* 2^-k */ - y = x - (e + t); - y += 1.0f; - y = y*twopk; - } + u.i = (0x7f-k)<<23; /* 2^-k */ + if (k < 23) + y = (x-e+(1-u.f))*twopk; + else + y = (x-(e+u.f)+1)*twopk; return y; } -- 2.25.1