diff options
Diffstat (limited to 'libc/sysdeps/ieee754/dbl-64/s_log1p.c')
| -rw-r--r-- | libc/sysdeps/ieee754/dbl-64/s_log1p.c | 52 |
1 files changed, 18 insertions, 34 deletions
diff --git a/libc/sysdeps/ieee754/dbl-64/s_log1p.c b/libc/sysdeps/ieee754/dbl-64/s_log1p.c index 0a9801a93..dc79a02bb 100644 --- a/libc/sysdeps/ieee754/dbl-64/s_log1p.c +++ b/libc/sysdeps/ieee754/dbl-64/s_log1p.c @@ -13,10 +13,6 @@ for performance improvement on pipelined processors. */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $"; -#endif - /* double log1p(double x) * * Method : @@ -34,14 +30,14 @@ static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $"; * 2. Approximation of log1p(f). * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) * = 2s + 2/3 s**3 + 2/5 s**5 + ....., - * = 2s + s*R + * = 2s + s*R * We use a special Reme algorithm on [0,0.1716] to generate - * a polynomial of degree 14 to approximate R The maximum error + * a polynomial of degree 14 to approximate R The maximum error * of this polynomial approximation is bounded by 2**-58.45. In * other words, - * 2 4 6 8 10 12 14 + * 2 4 6 8 10 12 14 * R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s +Lp6*s +Lp7*s - * (the values of Lp1 to Lp7 are listed in the program) + * (the values of Lp1 to Lp7 are listed in the program) * and * | 2 14 | -58.45 * | Lp1*s +...+Lp7*s - R(z) | <= 2 @@ -52,7 +48,7 @@ static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $"; * log1p(f) = f - (hfsq - s*(hfsq+R)). * * 3. Finally, log1p(x) = k*ln2 + log1p(f). - * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) + * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) * Here ln2 is split into two floating point number: * ln2_hi + ln2_lo, * where n*ln2_hi is always exact for |n| < 2000. @@ -73,7 +69,7 @@ static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $"; * to produce the hexadecimal values shown. * * Note: Assuming log() return accurate answer, the following - * algorithm can be used to compute log1p(x) to within a few ULP: + * algorithm can be used to compute log1p(x) to within a few ULP: * * u = 1+x; * if(u==1.0) return x ; else @@ -85,11 +81,7 @@ static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $"; #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const double -#else -static double -#endif ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ @@ -101,18 +93,10 @@ Lp[] = {0.0, 6.666666666666735130e-01, /* 3FE55555 55555593 */ 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ 1.479819860511658591e-01}; /* 3FC2F112 DF3E5244 */ -#ifdef __STDC__ static const double zero = 0.0; -#else -static double zero = 0.0; -#endif -#ifdef __STDC__ - double __log1p(double x) -#else - double __log1p(x) - double x; -#endif +double +__log1p(double x) { double hfsq,f,c,s,z,R,u,z2,z4,z6,R1,R2,R3,R4; int32_t k,hx,hu,ax; @@ -127,8 +111,8 @@ static double zero = 0.0; else return (x-x)/(x-x); /* log1p(x<-1)=NaN */ } if(ax<0x3e200000) { /* |x| < 2**-29 */ - if(two54+x>zero /* raise inexact */ - &&ax<0x3c900000) /* |x| < 2**-54 */ + math_force_eval(two54+x); /* raise inexact */ + if (ax<0x3c900000) /* |x| < 2**-54 */ return x; else return x - x*x*0.5; @@ -141,22 +125,22 @@ static double zero = 0.0; if(hx<0x43400000) { u = 1.0+x; GET_HIGH_WORD(hu,u); - k = (hu>>20)-1023; - c = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */ + k = (hu>>20)-1023; + c = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */ c /= u; } else { u = x; GET_HIGH_WORD(hu,u); - k = (hu>>20)-1023; + k = (hu>>20)-1023; c = 0; } hu &= 0x000fffff; if(hu<0x6a09e) { - SET_HIGH_WORD(u,hu|0x3ff00000); /* normalize u */ + SET_HIGH_WORD(u,hu|0x3ff00000); /* normalize u */ } else { - k += 1; + k += 1; SET_HIGH_WORD(u,hu|0x3fe00000); /* normalize u/2 */ - hu = (0x00100000-hu)>>2; + hu = (0x00100000-hu)>>2; } f = u-1.0; } @@ -168,9 +152,9 @@ static double zero = 0.0; } R = hfsq*(1.0-0.66666666666666666*f); if(k==0) return f-R; else - return k*ln2_hi-((R-(k*ln2_lo+c))-f); + return k*ln2_hi-((R-(k*ln2_lo+c))-f); } - s = f/(2.0+f); + s = f/(2.0+f); z = s*s; #ifdef DO_NOT_USE_THIS R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7)))))); |