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 * ====================================================
12 /* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/25,
13 for performance improvement on pipelined processors.
16 #if defined(LIBM_SCCS) && !defined(lint)
17 static char rcsid
[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $";
20 /* double log1p(double x)
23 * 1. Argument Reduction: find k and f such that
25 * where sqrt(2)/2 < 1+f < sqrt(2) .
27 * Note. If k=0, then f=x is exact. However, if k!=0, then f
28 * may not be representable exactly. In that case, a correction
29 * term is need. Let u=1+x rounded. Let c = (1+x)-u, then
30 * log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u),
31 * and add back the correction term c/u.
32 * (Note: when x > 2**53, one can simply return log(x))
34 * 2. Approximation of log1p(f).
35 * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
36 * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
38 * We use a special Reme algorithm on [0,0.1716] to generate
39 * a polynomial of degree 14 to approximate R The maximum error
40 * of this polynomial approximation is bounded by 2**-58.45. In
43 * R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s +Lp6*s +Lp7*s
44 * (the values of Lp1 to Lp7 are listed in the program)
47 * | Lp1*s +...+Lp7*s - R(z) | <= 2
49 * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
50 * In order to guarantee error in log below 1ulp, we compute log
52 * log1p(f) = f - (hfsq - s*(hfsq+R)).
54 * 3. Finally, log1p(x) = k*ln2 + log1p(f).
55 * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
56 * Here ln2 is split into two floating point number:
58 * where n*ln2_hi is always exact for |n| < 2000.
61 * log1p(x) is NaN with signal if x < -1 (including -INF) ;
62 * log1p(+INF) is +INF; log1p(-1) is -INF with signal;
63 * log1p(NaN) is that NaN with no signal.
66 * according to an error analysis, the error is always less than
67 * 1 ulp (unit in the last place).
70 * The hexadecimal values are the intended ones for the following
71 * constants. The decimal values may be used, provided that the
72 * compiler will convert from decimal to binary accurately enough
73 * to produce the hexadecimal values shown.
75 * Note: Assuming log() return accurate answer, the following
76 * algorithm can be used to compute log1p(x) to within a few ULP:
79 * if(u==1.0) return x ; else
80 * return log(u)*(x/(u-1.0));
82 * See HP-15C Advanced Functions Handbook, p.193.
86 #include "math_private.h"
93 ln2_hi
= 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
94 ln2_lo
= 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
95 two54
= 1.80143985094819840000e+16, /* 43500000 00000000 */
96 Lp
[] = {0.0, 6.666666666666735130e-01, /* 3FE55555 55555593 */
97 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
98 2.857142874366239149e-01, /* 3FD24924 94229359 */
99 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
100 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
101 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
102 1.479819860511658591e-01}; /* 3FC2F112 DF3E5244 */
105 static const double zero
= 0.0;
107 static double zero
= 0.0;
111 double __log1p(double x
)
117 double hfsq
,f
,c
,s
,z
,R
,u
,z2
,z4
,z6
,R1
,R2
,R3
,R4
;
124 if (hx
< 0x3FDA827A) { /* x < 0.41422 */
125 if(ax
>=0x3ff00000) { /* x <= -1.0 */
126 if(x
==-1.0) return -two54
/(x
-x
);/* log1p(-1)=+inf */
127 else return (x
-x
)/(x
-x
); /* log1p(x<-1)=NaN */
129 if(ax
<0x3e200000) { /* |x| < 2**-29 */
130 if(two54
+x
>zero
/* raise inexact */
131 &&ax
<0x3c900000) /* |x| < 2**-54 */
136 if(hx
>0||hx
<=((int32_t)0xbfd2bec3)) {
137 k
=0;f
=x
;hu
=1;} /* -0.2929<x<0.41422 */
139 if (hx
>= 0x7ff00000) return x
+x
;
145 c
= (k
>0)? 1.0-(u
-x
):x
-(u
-1.0);/* correction term */
155 SET_HIGH_WORD(u
,hu
|0x3ff00000); /* normalize u */
158 SET_HIGH_WORD(u
,hu
|0x3fe00000); /* normalize u/2 */
159 hu
= (0x00100000-hu
)>>2;
164 if(hu
==0) { /* |f| < 2**-20 */
166 if(k
==0) return zero
;
167 else {c
+= k
*ln2_lo
; return k
*ln2_hi
+c
;}
169 R
= hfsq
*(1.0-0.66666666666666666*f
);
170 if(k
==0) return f
-R
; else
171 return k
*ln2_hi
-((R
-(k
*ln2_lo
+c
))-f
);
175 #ifdef DO_NOT_USE_THIS
176 R
= z
*(Lp1
+z
*(Lp2
+z
*(Lp3
+z
*(Lp4
+z
*(Lp5
+z
*(Lp6
+z
*Lp7
))))));
178 R1
= z
*Lp
[1]; z2
=z
*z
;
179 R2
= Lp
[2]+z
*Lp
[3]; z4
=z2
*z2
;
180 R3
= Lp
[4]+z
*Lp
[5]; z6
=z4
*z2
;
182 R
= R1
+ z2
*R2
+ z4
*R3
+ z6
*R4
;
184 if(k
==0) return f
-(hfsq
-s
*(hfsq
+R
)); else
185 return k
*ln2_hi
-((hfsq
-(s
*(hfsq
+R
)+(k
*ln2_lo
+c
)))-f
);
187 weak_alias (__log1p
, log1p
)
188 #ifdef NO_LONG_DOUBLE
189 strong_alias (__log1p
, __log1pl
)
190 weak_alias (__log1p
, log1pl
)