1/* log1pl.c
2 *
3 * Relative error logarithm
4 * Natural logarithm of 1+x, 128-bit long double precision
5 *
6 *
7 *
8 * SYNOPSIS:
9 *
10 * long double x, y, log1pl();
11 *
12 * y = log1pl( x );
13 *
14 *
15 *
16 * DESCRIPTION:
17 *
18 * Returns the base e (2.718...) logarithm of 1+x.
19 *
20 * The argument 1+x is separated into its exponent and fractional
21 * parts. If the exponent is between -1 and +1, the logarithm
22 * of the fraction is approximated by
23 *
24 * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
25 *
26 * Otherwise, setting z = 2(w-1)/(w+1),
27 *
28 * log(w) = z + z^3 P(z)/Q(z).
29 *
30 *
31 *
32 * ACCURACY:
33 *
34 * Relative error:
35 * arithmetic domain # trials peak rms
36 * IEEE -1, 8 100000 1.9e-34 4.3e-35
37 */
38
39/* Copyright 2001 by Stephen L. Moshier
40
41 This library is free software; you can redistribute it and/or
42 modify it under the terms of the GNU Lesser General Public
43 License as published by the Free Software Foundation; either
44 version 2.1 of the License, or (at your option) any later version.
45
46 This library is distributed in the hope that it will be useful,
47 but WITHOUT ANY WARRANTY; without even the implied warranty of
48 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
49 Lesser General Public License for more details.
50
51 You should have received a copy of the GNU Lesser General Public
52 License along with this library; if not, see
53 <http://www.gnu.org/licenses/>. */
54
55
56#include <float.h>
57#include <math.h>
58#include <math_private.h>
59
60/* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
61 * 1/sqrt(2) <= 1+x < sqrt(2)
62 * Theoretical peak relative error = 5.3e-37,
63 * relative peak error spread = 2.3e-14
64 */
65static const _Float128
66 P12 = L(1.538612243596254322971797716843006400388E-6),
67 P11 = L(4.998469661968096229986658302195402690910E-1),
68 P10 = L(2.321125933898420063925789532045674660756E1),
69 P9 = L(4.114517881637811823002128927449878962058E2),
70 P8 = L(3.824952356185897735160588078446136783779E3),
71 P7 = L(2.128857716871515081352991964243375186031E4),
72 P6 = L(7.594356839258970405033155585486712125861E4),
73 P5 = L(1.797628303815655343403735250238293741397E5),
74 P4 = L(2.854829159639697837788887080758954924001E5),
75 P3 = L(3.007007295140399532324943111654767187848E5),
76 P2 = L(2.014652742082537582487669938141683759923E5),
77 P1 = L(7.771154681358524243729929227226708890930E4),
78 P0 = L(1.313572404063446165910279910527789794488E4),
79 /* Q12 = 1.000000000000000000000000000000000000000E0L, */
80 Q11 = L(4.839208193348159620282142911143429644326E1),
81 Q10 = L(9.104928120962988414618126155557301584078E2),
82 Q9 = L(9.147150349299596453976674231612674085381E3),
83 Q8 = L(5.605842085972455027590989944010492125825E4),
84 Q7 = L(2.248234257620569139969141618556349415120E5),
85 Q6 = L(6.132189329546557743179177159925690841200E5),
86 Q5 = L(1.158019977462989115839826904108208787040E6),
87 Q4 = L(1.514882452993549494932585972882995548426E6),
88 Q3 = L(1.347518538384329112529391120390701166528E6),
89 Q2 = L(7.777690340007566932935753241556479363645E5),
90 Q1 = L(2.626900195321832660448791748036714883242E5),
91 Q0 = L(3.940717212190338497730839731583397586124E4);
92
93/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
94 * where z = 2(x-1)/(x+1)
95 * 1/sqrt(2) <= x < sqrt(2)
96 * Theoretical peak relative error = 1.1e-35,
97 * relative peak error spread 1.1e-9
98 */
99static const _Float128
100 R5 = L(-8.828896441624934385266096344596648080902E-1),
101 R4 = L(8.057002716646055371965756206836056074715E1),
102 R3 = L(-2.024301798136027039250415126250455056397E3),
103 R2 = L(2.048819892795278657810231591630928516206E4),
104 R1 = L(-8.977257995689735303686582344659576526998E4),
105 R0 = L(1.418134209872192732479751274970992665513E5),
106 /* S6 = 1.000000000000000000000000000000000000000E0L, */
107 S5 = L(-1.186359407982897997337150403816839480438E2),
108 S4 = L(3.998526750980007367835804959888064681098E3),
109 S3 = L(-5.748542087379434595104154610899551484314E4),
110 S2 = L(4.001557694070773974936904547424676279307E5),
111 S1 = L(-1.332535117259762928288745111081235577029E6),
112 S0 = L(1.701761051846631278975701529965589676574E6);
113
114/* C1 + C2 = ln 2 */
115static const _Float128 C1 = L(6.93145751953125E-1);
116static const _Float128 C2 = L(1.428606820309417232121458176568075500134E-6);
117
118static const _Float128 sqrth = L(0.7071067811865475244008443621048490392848);
119/* ln (2^16384 * (1 - 2^-113)) */
120static const _Float128 zero = 0;
121
122_Float128
123__log1pl (_Float128 xm1)
124{
125 _Float128 x, y, z, r, s;
126 ieee854_long_double_shape_type u;
127 int32_t hx;
128 int e;
129
130 /* Test for NaN or infinity input. */
131 u.value = xm1;
132 hx = u.parts32.w0;
133 if ((hx & 0x7fffffff) >= 0x7fff0000)
134 return xm1 + fabsl (xm1);
135
136 /* log1p(+- 0) = +- 0. */
137 if (((hx & 0x7fffffff) == 0)
138 && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
139 return xm1;
140
141 if ((hx & 0x7fffffff) < 0x3f8e0000)
142 {
143 math_check_force_underflow (xm1);
144 if ((int) xm1 == 0)
145 return xm1;
146 }
147
148 if (xm1 >= L(0x1p113))
149 x = xm1;
150 else
151 x = xm1 + 1;
152
153 /* log1p(-1) = -inf */
154 if (x <= 0)
155 {
156 if (x == 0)
157 return (-1 / zero); /* log1p(-1) = -inf */
158 else
159 return (zero / (x - x));
160 }
161
162 /* Separate mantissa from exponent. */
163
164 /* Use frexp used so that denormal numbers will be handled properly. */
165 x = __frexpl (x, &e);
166
167 /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
168 where z = 2(x-1)/x+1). */
169 if ((e > 2) || (e < -2))
170 {
171 if (x < sqrth)
172 { /* 2( 2x-1 )/( 2x+1 ) */
173 e -= 1;
174 z = x - L(0.5);
175 y = L(0.5) * z + L(0.5);
176 }
177 else
178 { /* 2 (x-1)/(x+1) */
179 z = x - L(0.5);
180 z -= L(0.5);
181 y = L(0.5) * x + L(0.5);
182 }
183 x = z / y;
184 z = x * x;
185 r = ((((R5 * z
186 + R4) * z
187 + R3) * z
188 + R2) * z
189 + R1) * z
190 + R0;
191 s = (((((z
192 + S5) * z
193 + S4) * z
194 + S3) * z
195 + S2) * z
196 + S1) * z
197 + S0;
198 z = x * (z * r / s);
199 z = z + e * C2;
200 z = z + x;
201 z = z + e * C1;
202 return (z);
203 }
204
205
206 /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
207
208 if (x < sqrth)
209 {
210 e -= 1;
211 if (e != 0)
212 x = 2 * x - 1; /* 2x - 1 */
213 else
214 x = xm1;
215 }
216 else
217 {
218 if (e != 0)
219 x = x - 1;
220 else
221 x = xm1;
222 }
223 z = x * x;
224 r = (((((((((((P12 * x
225 + P11) * x
226 + P10) * x
227 + P9) * x
228 + P8) * x
229 + P7) * x
230 + P6) * x
231 + P5) * x
232 + P4) * x
233 + P3) * x
234 + P2) * x
235 + P1) * x
236 + P0;
237 s = (((((((((((x
238 + Q11) * x
239 + Q10) * x
240 + Q9) * x
241 + Q8) * x
242 + Q7) * x
243 + Q6) * x
244 + Q5) * x
245 + Q4) * x
246 + Q3) * x
247 + Q2) * x
248 + Q1) * x
249 + Q0;
250 y = x * (z * r / s);
251 y = y + e * C2;
252 z = y - L(0.5) * z;
253 z = z + x;
254 z = z + e * C1;
255 return (z);
256}
257