1 /* @(#)s_log1p.c 5.1 93/09/24 */
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
10 * ====================================================
13 #if defined(LIBM_SCCS) && !defined(lint)
14 static char rcsid
[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $";
17 /* double log1p(double x)
20 * 1. Argument Reduction: find k and f such that
22 * where sqrt(2)/2 < 1+f < sqrt(2) .
24 * Note. If k=0, then f=x is exact. However, if k!=0, then f
25 * may not be representable exactly. In that case, a correction
26 * term is need. Let u=1+x rounded. Let c = (1+x)-u, then
27 * log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u),
28 * and add back the correction term c/u.
29 * (Note: when x > 2**53, one can simply return log(x))
31 * 2. Approximation of log1p(f).
32 * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
33 * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
35 * We use a special Reme algorithm on [0,0.1716] to generate
36 * a polynomial of degree 14 to approximate R The maximum error
37 * of this polynomial approximation is bounded by 2**-58.45. In
40 * R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s +Lp6*s +Lp7*s
41 * (the values of Lp1 to Lp7 are listed in the program)
44 * | Lp1*s +...+Lp7*s - R(z) | <= 2
46 * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
47 * In order to guarantee error in log below 1ulp, we compute log
49 * log1p(f) = f - (hfsq - s*(hfsq+R)).
51 * 3. Finally, log1p(x) = k*ln2 + log1p(f).
52 * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
53 * Here ln2 is split into two floating point number:
55 * where n*ln2_hi is always exact for |n| < 2000.
58 * log1p(x) is NaN with signal if x < -1 (including -INF) ;
59 * log1p(+INF) is +INF; log1p(-1) is -INF with signal;
60 * log1p(NaN) is that NaN with no signal.
63 * according to an error analysis, the error is always less than
64 * 1 ulp (unit in the last place).
67 * The hexadecimal values are the intended ones for the following
68 * constants. The decimal values may be used, provided that the
69 * compiler will convert from decimal to binary accurately enough
70 * to produce the hexadecimal values shown.
72 * Note: Assuming log() return accurate answer, the following
73 * algorithm can be used to compute log1p(x) to within a few ULP:
76 * if(u==1.0) return x ; else
77 * return log(u)*(x/(u-1.0));
79 * See HP-15C Advanced Functions Handbook, p.193.
83 #include "math_private.h"
90 ln2_hi
= 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
91 ln2_lo
= 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
92 two54
= 1.80143985094819840000e+16, /* 43500000 00000000 */
93 Lp1
= 6.666666666666735130e-01, /* 3FE55555 55555593 */
94 Lp2
= 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
95 Lp3
= 2.857142874366239149e-01, /* 3FD24924 94229359 */
96 Lp4
= 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
97 Lp5
= 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
98 Lp6
= 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
99 Lp7
= 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
102 static const double zero
= 0.0;
104 static double zero
= 0.0;
108 double __log1p(double x
)
114 double hfsq
,f
,c
,s
,z
,R
,u
;
121 if (hx
< 0x3FDA827A) { /* x < 0.41422 */
122 if(ax
>=0x3ff00000) { /* x <= -1.0 */
123 if(x
==-1.0) return -two54
/zero
; /* log1p(-1)=+inf */
124 else return (x
-x
)/(x
-x
); /* log1p(x<-1)=NaN */
126 if(ax
<0x3e200000) { /* |x| < 2**-29 */
127 if(two54
+x
>zero
/* raise inexact */
128 &&ax
<0x3c900000) /* |x| < 2**-54 */
133 if(hx
>0||hx
<=((int32_t)0xbfd2bec3)) {
134 k
=0;f
=x
;hu
=1;} /* -0.2929<x<0.41422 */
136 if (hx
>= 0x7ff00000) return x
+x
;
142 c
= (k
>0)? 1.0-(u
-x
):x
-(u
-1.0);/* correction term */
152 SET_HIGH_WORD(u
,hu
|0x3ff00000); /* normalize u */
155 SET_HIGH_WORD(u
,hu
|0x3fe00000); /* normalize u/2 */
156 hu
= (0x00100000-hu
)>>2;
161 if(hu
==0) { /* |f| < 2**-20 */
162 if(f
==zero
) if(k
==0) return zero
;
163 else {c
+= k
*ln2_lo
; return k
*ln2_hi
+c
;}
164 R
= hfsq
*(1.0-0.66666666666666666*f
);
165 if(k
==0) return f
-R
; else
166 return k
*ln2_hi
-((R
-(k
*ln2_lo
+c
))-f
);
170 R
= z
*(Lp1
+z
*(Lp2
+z
*(Lp3
+z
*(Lp4
+z
*(Lp5
+z
*(Lp6
+z
*Lp7
))))));
171 if(k
==0) return f
-(hfsq
-s
*(hfsq
+R
)); else
172 return k
*ln2_hi
-((hfsq
-(s
*(hfsq
+R
)+(k
*ln2_lo
+c
)))-f
);
174 weak_alias (__log1p
, log1p
)
175 #ifdef NO_LONG_DOUBLE
176 strong_alias (__log1p
, __log1pl
)
177 weak_alias (__log1p
, log1pl
)