Print this page
*** 158,168 ****
#define Q4 xxx[10]
#define Q5 xxx[11]
double
expm1(double x) {
! double y, hi, lo, c, t, e, hxs, hfx, r1;
int k, xsb;
unsigned hx;
hx = ((unsigned *) &x)[HIWORD]; /* high word of x */
xsb = hx & 0x80000000; /* sign bit of x */
--- 158,168 ----
#define Q4 xxx[10]
#define Q5 xxx[11]
double
expm1(double x) {
! double y, hi, lo, c = 0.0L, t, e, hxs, hfx, r1;
int k, xsb;
unsigned hx;
hx = ((unsigned *) &x)[HIWORD]; /* high word of x */
xsb = hx & 0x80000000; /* sign bit of x */
*** 170,182 ****
y = x;
else
y = -x; /* y = |x| */
hx &= 0x7fffffff; /* high word of |x| */
! /* filter out huge and non-finite arugment */
! if (hx >= 0x4043687A) { /* if |x|>=56*ln2 */
! if (hx >= 0x40862E42) { /* if |x|>=709.78... */
if (hx >= 0x7ff00000) {
if (((hx & 0xfffff) | ((int *) &x)[LOWORD])
!= 0)
return x * x; /* + -> * for Cheetah */
else
--- 170,183 ----
y = x;
else
y = -x; /* y = |x| */
hx &= 0x7fffffff; /* high word of |x| */
! /* filter out huge and non-finite argument */
! /* for example exp(38)-1 is approximately 3.1855932e+16 */
! if (hx >= 0x4043687A) { /* if |x|>=56*ln2 (~38.8162...)*/
! if (hx >= 0x40862E42) { /* if |x|>=709.78... -> inf */
if (hx >= 0x7ff00000) {
if (((hx & 0xfffff) | ((int *) &x)[LOWORD])
!= 0)
return x * x; /* + -> * for Cheetah */
else
*** 192,246 ****
}
/* 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 = (int) (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;
}
x = hi - lo;
! c = (hi - x) - lo;
}
else if (hx < 0x3c900000) { /* when |x|<2**-54, return x */
t = huge + x; /* return x w/inexact 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 * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5))));
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 */
y = one - (e - x);
((int *) &y)[HIWORD] += k << 20;
return y - one;
}
--- 193,248 ----
}
/* argument reduction */
if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
! if (xsb == 0) { /* positive number */
hi = x - ln2_hi;
lo = ln2_lo;
k = 1;
}
! else { /* negative number */
hi = x + ln2_hi;
lo = -ln2_lo;
k = -1;
}
}
! else { /* |x| > 1.5 ln2 */
k = (int) (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;
}
x = hi - lo;
! c = (hi - x) - lo; /* still at |x| > 0.5 ln2 */
}
else if (hx < 0x3c900000) { /* when |x|<2**-54, return x */
t = huge + x; /* return x w/inexact when x != 0 */
return x - (t - (huge + x));
}
! else /* |x| <= 0.5 ln2 */
k = 0;
/* x is now in primary range */
hfx = 0.5 * x;
hxs = x * hfx;
r1 = one + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5))));
t = 3.0 - r1 * hfx;
e = hxs * ((r1 - t) / (6.0 - x * t));
! if (k == 0) /* |x| <= 0.5 ln2 */
! return x - (x * e - hxs);
! else { /* |x| > 0.5 ln2 */
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 */
y = one - (e - x);
((int *) &y)[HIWORD] += k << 20;
return y - one;
}