summaryrefslogtreecommitdiff
path: root/src/math/expm1.c
diff options
context:
space:
mode:
authorSzabolcs Nagy <nsz@port70.net>2013-09-04 07:51:11 +0000
committerSzabolcs Nagy <nsz@port70.net>2013-09-05 11:30:08 +0000
commit39c910fb061114e6aa5c3bf2c94b1d7262d62221 (patch)
tree5ca9b82746e8ac224f93dcda1b8d19f95b68eda2 /src/math/expm1.c
parentea9bb95a5b36c0a3d2ed8fb03808745b406c2633 (diff)
downloadmusl-39c910fb061114e6aa5c3bf2c94b1d7262d62221.tar.gz
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)
Diffstat (limited to 'src/math/expm1.c')
-rw-r--r--src/math/expm1.c72
1 files changed, 27 insertions, 45 deletions
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;
}