diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/s_expm1.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_expm1.c | 223 |
1 files changed, 126 insertions, 97 deletions
diff --git a/sysdeps/ieee754/dbl-64/s_expm1.c b/sysdeps/ieee754/dbl-64/s_expm1.c index 87da7eb7e6..9c9bb34559 100644 --- a/sysdeps/ieee754/dbl-64/s_expm1.c +++ b/sysdeps/ieee754/dbl-64/s_expm1.c @@ -11,7 +11,7 @@ */ /* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/25, for performance improvement on pipelined processors. -*/ + */ /* expm1(x) * Returns exp(x)-1, the exponential of x minus 1. @@ -113,116 +113,145 @@ #include <math_private.h> #define one Q[0] 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 */ -invln2 = 1.44269504088896338700e+00,/* 0x3ff71547, 0x652b82fe */ - /* scaled coefficients related to expm1 */ -Q[] = {1.0, -3.33333333333331316428e-02, /* BFA11111 111110F4 */ - 1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */ - -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */ - 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */ - -2.01099218183624371326e-07}; /* BE8AFDB7 6E09C32D */ + 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 */ + invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */ +/* scaled coefficients related to expm1 */ + Q[] = { 1.0, -3.33333333333331316428e-02, /* BFA11111 111110F4 */ + 1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */ + -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */ + 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */ + -2.01099218183624371326e-07 }; /* BE8AFDB7 6E09C32D */ double -__expm1(double x) +__expm1 (double x) { - double y,hi,lo,c,t,e,hxs,hfx,r1,h2,h4,R1,R2,R3; - int32_t k,xsb; - u_int32_t hx; + double y, hi, lo, c, t, e, hxs, hfx, r1, h2, h4, R1, R2, R3; + int32_t k, xsb; + u_int32_t hx; - GET_HIGH_WORD(hx,x); - xsb = hx&0x80000000; /* sign bit of x */ - if(xsb==0) y=x; else y= -x; /* y = |x| */ - hx &= 0x7fffffff; /* high word of |x| */ + GET_HIGH_WORD (hx, x); + xsb = hx & 0x80000000; /* sign bit of x */ + if (xsb == 0) + y = x; + else + y = -x; /* y = |x| */ + hx &= 0x7fffffff; /* high word of |x| */ - /* filter out huge and non-finite argument */ - if(hx >= 0x4043687A) { /* if |x|>=56*ln2 */ - if(hx >= 0x40862E42) { /* if |x|>=709.78... */ - if(hx>=0x7ff00000) { - u_int32_t low; - GET_LOW_WORD(low,x); - if(((hx&0xfffff)|low)!=0) - return x+x; /* NaN */ - else return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */ - } - if(x > o_threshold) { - __set_errno (ERANGE); - return huge*huge; /* overflow */ - } + /* filter out huge and non-finite argument */ + if (hx >= 0x4043687A) /* if |x|>=56*ln2 */ + { + if (hx >= 0x40862E42) /* if |x|>=709.78... */ + { + if (hx >= 0x7ff00000) + { + u_int32_t low; + GET_LOW_WORD (low, x); + if (((hx & 0xfffff) | low) != 0) + return x + x; /* NaN */ + else + return (xsb == 0) ? x : -1.0; /* exp(+-inf)={inf,-1} */ } - if(xsb!=0) { /* x < -56*ln2, return -1.0 with inexact */ - math_force_eval(x+tiny); /* raise inexact */ - return tiny-one; /* return -1 */ + if (x > o_threshold) + { + __set_errno (ERANGE); + return huge * huge; /* overflow */ } } + if (xsb != 0) /* x < -56*ln2, return -1.0 with inexact */ + { + math_force_eval (x + tiny); /* raise inexact */ + return tiny - one; /* return -1 */ + } + } - /* argument reduction */ - if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ - if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ - if(xsb==0) - {hi = x - ln2_hi; lo = ln2_lo; k = 1;} - else - {hi = x + ln2_hi; lo = -ln2_lo; k = -1;} - } else { - k = invln2*x+((xsb==0)?0.5:-0.5); - t = k; - hi = x - t*ln2_hi; /* t*ln2_hi is exact here */ - lo = t*ln2_lo; + /* argument reduction */ + if (hx > 0x3fd62e42) /* if |x| > 0.5 ln2 */ + { + if (hx < 0x3FF0A2B2) /* and |x| < 1.5 ln2 */ + { + if (xsb == 0) + { + hi = x - ln2_hi; lo = ln2_lo; k = 1; + } + else + { + hi = x + ln2_hi; lo = -ln2_lo; k = -1; } - x = hi - lo; - c = (hi-x)-lo; } - else if(hx < 0x3c900000) { /* when |x|<2**-54, return x */ - t = huge+x; /* return x with inexact flags when x!=0 */ - return x - (t-(huge+x)); + else + { + k = invln2 * x + ((xsb == 0) ? 0.5 : -0.5); + t = k; + hi = x - t * ln2_hi; /* t*ln2_hi is exact here */ + lo = t * ln2_lo; } - else k = 0; + x = hi - lo; + c = (hi - x) - lo; + } + else if (hx < 0x3c900000) /* when |x|<2**-54, return x */ + { + t = huge + x; /* return x with inexact flags when x!=0 */ + return x - (t - (huge + x)); + } + else + k = 0; - /* x is now in primary range */ - hfx = 0.5*x; - hxs = x*hfx; - R1 = one+hxs*Q[1]; h2 = hxs*hxs; - R2 = Q[2]+hxs*Q[3]; h4 = h2*h2; - R3 = Q[4]+hxs*Q[5]; - r1 = R1 + h2*R2 + h4*R3; - t = 3.0-r1*hfx; - e = hxs*((r1-t)/(6.0 - x*t)); - if(k==0) return x - (x*e-hxs); /* c is 0 */ - else { - e = (x*(e-c)-c); - e -= hxs; - if(k== -1) return 0.5*(x-e)-0.5; - if(k==1) { - if(x < -0.25) return -2.0*(e-(x+0.5)); - else return one+2.0*(x-e); - } - if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */ - u_int32_t high; - y = one-(e-x); - GET_HIGH_WORD(high,y); - SET_HIGH_WORD(y,high+(k<<20)); /* add k to y's exponent */ - return y-one; - } - t = one; - if(k<20) { - u_int32_t high; - SET_HIGH_WORD(t,0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */ - y = t-(e-x); - GET_HIGH_WORD(high,y); - SET_HIGH_WORD(y,high+(k<<20)); /* add k to y's exponent */ - } else { - u_int32_t high; - SET_HIGH_WORD(t,((0x3ff-k)<<20)); /* 2^-k */ - y = x-(e+t); - y += one; - GET_HIGH_WORD(high,y); - SET_HIGH_WORD(y,high+(k<<20)); /* add k to y's exponent */ - } + /* x is now in primary range */ + hfx = 0.5 * x; + hxs = x * hfx; + R1 = one + hxs * Q[1]; h2 = hxs * hxs; + R2 = Q[2] + hxs * Q[3]; h4 = h2 * h2; + R3 = Q[4] + hxs * Q[5]; + r1 = R1 + h2 * R2 + h4 * R3; + t = 3.0 - r1 * hfx; + e = hxs * ((r1 - t) / (6.0 - x * t)); + if (k == 0) + return x - (x * e - hxs); /* c is 0 */ + else + { + e = (x * (e - c) - c); + e -= hxs; + if (k == -1) + return 0.5 * (x - e) - 0.5; + if (k == 1) + { + if (x < -0.25) + return -2.0 * (e - (x + 0.5)); + else + return one + 2.0 * (x - e); + } + if (k <= -2 || k > 56) /* suffice to return exp(x)-1 */ + { + u_int32_t high; + y = one - (e - x); + GET_HIGH_WORD (high, y); + SET_HIGH_WORD (y, high + (k << 20)); /* add k to y's exponent */ + return y - one; + } + t = one; + if (k < 20) + { + u_int32_t high; + SET_HIGH_WORD (t, 0x3ff00000 - (0x200000 >> k)); /* t=1-2^-k */ + y = t - (e - x); + GET_HIGH_WORD (high, y); + SET_HIGH_WORD (y, high + (k << 20)); /* add k to y's exponent */ + } + else + { + u_int32_t high; + SET_HIGH_WORD (t, ((0x3ff - k) << 20)); /* 2^-k */ + y = x - (e + t); + y += one; + GET_HIGH_WORD (high, y); + SET_HIGH_WORD (y, high + (k << 20)); /* add k to y's exponent */ } - return y; + } + return y; } weak_alias (__expm1, expm1) #ifdef NO_LONG_DOUBLE |