1 | /* @(#)s_log1p.c 5.1 93/09/24 */ |
2 | /* |
3 | * ==================================================== |
4 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
5 | * |
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 |
9 | * is preserved. |
10 | * ==================================================== |
11 | */ |
12 | /* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/25, |
13 | for performance improvement on pipelined processors. |
14 | */ |
15 | |
16 | /* double log1p(double x) |
17 | * |
18 | * Method : |
19 | * 1. Argument Reduction: find k and f such that |
20 | * 1+x = 2^k * (1+f), |
21 | * where sqrt(2)/2 < 1+f < sqrt(2) . |
22 | * |
23 | * Note. If k=0, then f=x is exact. However, if k!=0, then f |
24 | * may not be representable exactly. In that case, a correction |
25 | * term is need. Let u=1+x rounded. Let c = (1+x)-u, then |
26 | * log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u), |
27 | * and add back the correction term c/u. |
28 | * (Note: when x > 2**53, one can simply return log(x)) |
29 | * |
30 | * 2. Approximation of log1p(f). |
31 | * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) |
32 | * = 2s + 2/3 s**3 + 2/5 s**5 + ....., |
33 | * = 2s + s*R |
34 | * We use a special Reme algorithm on [0,0.1716] to generate |
35 | * a polynomial of degree 14 to approximate R The maximum error |
36 | * of this polynomial approximation is bounded by 2**-58.45. In |
37 | * other words, |
38 | * 2 4 6 8 10 12 14 |
39 | * R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s +Lp6*s +Lp7*s |
40 | * (the values of Lp1 to Lp7 are listed in the program) |
41 | * and |
42 | * | 2 14 | -58.45 |
43 | * | Lp1*s +...+Lp7*s - R(z) | <= 2 |
44 | * | | |
45 | * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. |
46 | * In order to guarantee error in log below 1ulp, we compute log |
47 | * by |
48 | * log1p(f) = f - (hfsq - s*(hfsq+R)). |
49 | * |
50 | * 3. Finally, log1p(x) = k*ln2 + log1p(f). |
51 | * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) |
52 | * Here ln2 is split into two floating point number: |
53 | * ln2_hi + ln2_lo, |
54 | * where n*ln2_hi is always exact for |n| < 2000. |
55 | * |
56 | * Special cases: |
57 | * log1p(x) is NaN with signal if x < -1 (including -INF) ; |
58 | * log1p(+INF) is +INF; log1p(-1) is -INF with signal; |
59 | * log1p(NaN) is that NaN with no signal. |
60 | * |
61 | * Accuracy: |
62 | * according to an error analysis, the error is always less than |
63 | * 1 ulp (unit in the last place). |
64 | * |
65 | * Constants: |
66 | * The hexadecimal values are the intended ones for the following |
67 | * constants. The decimal values may be used, provided that the |
68 | * compiler will convert from decimal to binary accurately enough |
69 | * to produce the hexadecimal values shown. |
70 | * |
71 | * Note: Assuming log() return accurate answer, the following |
72 | * algorithm can be used to compute log1p(x) to within a few ULP: |
73 | * |
74 | * u = 1+x; |
75 | * if(u==1.0) return x ; else |
76 | * return log(u)*(x/(u-1.0)); |
77 | * |
78 | * See HP-15C Advanced Functions Handbook, p.193. |
79 | */ |
80 | |
81 | #include <float.h> |
82 | #include <math.h> |
83 | #include <math_private.h> |
84 | |
85 | static const double |
86 | ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ |
87 | ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ |
88 | two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ |
89 | Lp[] = { 0.0, 6.666666666666735130e-01, /* 3FE55555 55555593 */ |
90 | 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ |
91 | 2.857142874366239149e-01, /* 3FD24924 94229359 */ |
92 | 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ |
93 | 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ |
94 | 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ |
95 | 1.479819860511658591e-01 }; /* 3FC2F112 DF3E5244 */ |
96 | |
97 | static const double zero = 0.0; |
98 | |
99 | double |
100 | __log1p (double x) |
101 | { |
102 | double hfsq, f, c, s, z, R, u, z2, z4, z6, R1, R2, R3, R4; |
103 | int32_t k, hx, hu, ax; |
104 | |
105 | GET_HIGH_WORD (hx, x); |
106 | ax = hx & 0x7fffffff; |
107 | |
108 | k = 1; |
109 | if (hx < 0x3FDA827A) /* x < 0.41422 */ |
110 | { |
111 | if (__glibc_unlikely (ax >= 0x3ff00000)) /* x <= -1.0 */ |
112 | { |
113 | if (x == -1.0) |
114 | return -two54 / zero; /* log1p(-1)=-inf */ |
115 | else |
116 | return (x - x) / (x - x); /* log1p(x<-1)=NaN */ |
117 | } |
118 | if (__glibc_unlikely (ax < 0x3e200000)) /* |x| < 2**-29 */ |
119 | { |
120 | math_force_eval (two54 + x); /* raise inexact */ |
121 | if (ax < 0x3c900000) /* |x| < 2**-54 */ |
122 | { |
123 | math_check_force_underflow (x); |
124 | return x; |
125 | } |
126 | else |
127 | return x - x * x * 0.5; |
128 | } |
129 | if (hx > 0 || hx <= ((int32_t) 0xbfd2bec3)) |
130 | { |
131 | k = 0; f = x; hu = 1; |
132 | } /* -0.2929<x<0.41422 */ |
133 | } |
134 | else if (__glibc_unlikely (hx >= 0x7ff00000)) |
135 | return x + x; |
136 | if (k != 0) |
137 | { |
138 | if (hx < 0x43400000) |
139 | { |
140 | u = 1.0 + x; |
141 | GET_HIGH_WORD (hu, u); |
142 | k = (hu >> 20) - 1023; |
143 | c = (k > 0) ? 1.0 - (u - x) : x - (u - 1.0); /* correction term */ |
144 | c /= u; |
145 | } |
146 | else |
147 | { |
148 | u = x; |
149 | GET_HIGH_WORD (hu, u); |
150 | k = (hu >> 20) - 1023; |
151 | c = 0; |
152 | } |
153 | hu &= 0x000fffff; |
154 | if (hu < 0x6a09e) |
155 | { |
156 | SET_HIGH_WORD (u, hu | 0x3ff00000); /* normalize u */ |
157 | } |
158 | else |
159 | { |
160 | k += 1; |
161 | SET_HIGH_WORD (u, hu | 0x3fe00000); /* normalize u/2 */ |
162 | hu = (0x00100000 - hu) >> 2; |
163 | } |
164 | f = u - 1.0; |
165 | } |
166 | hfsq = 0.5 * f * f; |
167 | if (hu == 0) /* |f| < 2**-20 */ |
168 | { |
169 | if (f == zero) |
170 | { |
171 | if (k == 0) |
172 | return zero; |
173 | else |
174 | { |
175 | c += k * ln2_lo; return k * ln2_hi + c; |
176 | } |
177 | } |
178 | R = hfsq * (1.0 - 0.66666666666666666 * f); |
179 | if (k == 0) |
180 | return f - R; |
181 | else |
182 | return k * ln2_hi - ((R - (k * ln2_lo + c)) - f); |
183 | } |
184 | s = f / (2.0 + f); |
185 | z = s * s; |
186 | R1 = z * Lp[1]; z2 = z * z; |
187 | R2 = Lp[2] + z * Lp[3]; z4 = z2 * z2; |
188 | R3 = Lp[4] + z * Lp[5]; z6 = z4 * z2; |
189 | R4 = Lp[6] + z * Lp[7]; |
190 | R = R1 + z2 * R2 + z4 * R3 + z6 * R4; |
191 | if (k == 0) |
192 | return f - (hfsq - s * (hfsq + R)); |
193 | else |
194 | return k * ln2_hi - ((hfsq - (s * (hfsq + R) + (k * ln2_lo + c))) - f); |
195 | } |
196 | |