diff options
Diffstat (limited to 'lib/libm/common_source/expm1.c')
| -rw-r--r-- | lib/libm/common_source/expm1.c | 32 |
1 files changed, 16 insertions, 16 deletions
diff --git a/lib/libm/common_source/expm1.c b/lib/libm/common_source/expm1.c index 760d2be..8a84f14 100644 --- a/lib/libm/common_source/expm1.c +++ b/lib/libm/common_source/expm1.c @@ -38,36 +38,36 @@ static char sccsid[] = "@(#)expm1.c 8.1 (Berkeley) 6/4/93"; /* EXPM1(X) * RETURN THE EXPONENTIAL OF X MINUS ONE * DOUBLE PRECISION (IEEE 53 BITS, VAX D FORMAT 56 BITS) - * CODED IN C BY K.C. NG, 1/19/85; + * CODED IN C BY K.C. NG, 1/19/85; * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/21/85, 4/16/85. * * Required system supported functions: - * scalb(x,n) - * copysign(x,y) + * scalb(x,n) + * copysign(x,y) * finite(x) * * Kernel function: * exp__E(x,c) * * Method: - * 1. Argument Reduction: given the input x, find r and integer k such + * 1. Argument Reduction: given the input x, find r and integer k such * that - * x = k*ln2 + r, |r| <= 0.5*ln2 . + * x = k*ln2 + r, |r| <= 0.5*ln2 . * r will be represented as r := z+c for better accuracy. * - * 2. Compute EXPM1(r)=exp(r)-1 by + * 2. Compute EXPM1(r)=exp(r)-1 by * * EXPM1(r=z+c) := z + exp__E(z,c) * * 3. EXPM1(x) = 2^k * ( EXPM1(r) + 1-2^-k ). * - * Remarks: + * Remarks: * 1. When k=1 and z < -0.25, we use the following formula for * better accuracy: * EXPM1(x) = 2 * ( (z+0.5) + exp__E(z,c) ) * 2. To avoid rounding error in 1-2^-k where k is large, we use * EXPM1(x) = 2^k * { [z+(exp__E(z,c)-2^-k )] + 1 } - * when k>56. + * when k>56. * * Special cases: * EXPM1(INF) is INF, EXPM1(NaN) is NaN; @@ -108,7 +108,7 @@ ic(invln2, 1.4426950408889633870E0, 0, 1.71547652B82FE) double expm1(x) double x; { - const static double one=1.0, half=1.0/2.0; + const static double one=1.0, half=1.0/2.0; double z,hi,lo,c; int k; #if defined(vax)||defined(tahoe) @@ -126,13 +126,13 @@ double x; /* argument reduction : x - k*ln2 */ k= invln2 *x+copysign(0.5,x); /* k=NINT(x/ln2) */ - hi=x-k*ln2hi ; + hi=x-k*ln2hi ; z=hi-(lo=k*ln2lo); c=(hi-z)-lo; if(k==0) return(z+__exp__E(z,c)); if(k==1) - if(z< -0.25) + if(z< -0.25) {x=z+half;x +=__exp__E(z,c); return(x+x);} else {z+=__exp__E(z,c); x=half+z; return(x+x);} @@ -143,17 +143,17 @@ double x; { x=one-scalb(one,-k); z += __exp__E(z,c);} else if(k<100) { x = __exp__E(z,c)-scalb(one,-k); x+=z; z=one;} - else + else { x = __exp__E(z,c)+z; z=one;} - return (scalb(x+z,k)); + return (scalb(x+z,k)); } } /* end of x > lnunfl */ - else + else /* expm1(-big#) rounded to -1 (inexact) */ - if(finite(x)) + if(finite(x)) { ln2hi+ln2lo; return(-one);} /* expm1(-INF) is -1 */ @@ -161,7 +161,7 @@ double x; } /* end of x < lnhuge */ - else + else /* expm1(INF) is INF, expm1(+big#) overflows to INF */ return( finite(x) ? scalb(one,5000) : x); } |
