summaryrefslogtreecommitdiff log msg author committer range
path: root/src/math/expm1.c
diff options
 context: 12345678910152025303540 space: includeignore mode: unifiedssdiffstat only
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.cindex 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; }