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 <https://www.gnu.org/licenses/>. */
54
55
56#include <float.h>
57#include <math.h>
58#include <math_private.h>
59#include <math-underflow.h>
60
61/* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
62 * 1/sqrt(2) <= 1+x < sqrt(2)
63 * Theoretical peak relative error = 5.3e-37,
64 * relative peak error spread = 2.3e-14
65 */
66static const _Float128
67 P12 = L(1.538612243596254322971797716843006400388E-6),
68 P11 = L(4.998469661968096229986658302195402690910E-1),
69 P10 = L(2.321125933898420063925789532045674660756E1),
70 P9 = L(4.114517881637811823002128927449878962058E2),
71 P8 = L(3.824952356185897735160588078446136783779E3),
72 P7 = L(2.128857716871515081352991964243375186031E4),
73 P6 = L(7.594356839258970405033155585486712125861E4),
74 P5 = L(1.797628303815655343403735250238293741397E5),
75 P4 = L(2.854829159639697837788887080758954924001E5),
76 P3 = L(3.007007295140399532324943111654767187848E5),
77 P2 = L(2.014652742082537582487669938141683759923E5),
78 P1 = L(7.771154681358524243729929227226708890930E4),
79 P0 = L(1.313572404063446165910279910527789794488E4),
80 /* Q12 = 1.000000000000000000000000000000000000000E0L, */
81 Q11 = L(4.839208193348159620282142911143429644326E1),
82 Q10 = L(9.104928120962988414618126155557301584078E2),
83 Q9 = L(9.147150349299596453976674231612674085381E3),
84 Q8 = L(5.605842085972455027590989944010492125825E4),
85 Q7 = L(2.248234257620569139969141618556349415120E5),
86 Q6 = L(6.132189329546557743179177159925690841200E5),
87 Q5 = L(1.158019977462989115839826904108208787040E6),
88 Q4 = L(1.514882452993549494932585972882995548426E6),
89 Q3 = L(1.347518538384329112529391120390701166528E6),
90 Q2 = L(7.777690340007566932935753241556479363645E5),
91 Q1 = L(2.626900195321832660448791748036714883242E5),
92 Q0 = L(3.940717212190338497730839731583397586124E4);
93
94/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
95 * where z = 2(x-1)/(x+1)
96 * 1/sqrt(2) <= x < sqrt(2)
97 * Theoretical peak relative error = 1.1e-35,
98 * relative peak error spread 1.1e-9
99 */
100static const _Float128
101 R5 = L(-8.828896441624934385266096344596648080902E-1),
102 R4 = L(8.057002716646055371965756206836056074715E1),
103 R3 = L(-2.024301798136027039250415126250455056397E3),
104 R2 = L(2.048819892795278657810231591630928516206E4),
105 R1 = L(-8.977257995689735303686582344659576526998E4),
106 R0 = L(1.418134209872192732479751274970992665513E5),
107 /* S6 = 1.000000000000000000000000000000000000000E0L, */
108 S5 = L(-1.186359407982897997337150403816839480438E2),
109 S4 = L(3.998526750980007367835804959888064681098E3),
110 S3 = L(-5.748542087379434595104154610899551484314E4),
111 S2 = L(4.001557694070773974936904547424676279307E5),
112 S1 = L(-1.332535117259762928288745111081235577029E6),
113 S0 = L(1.701761051846631278975701529965589676574E6);
114
115/* C1 + C2 = ln 2 */
116static const _Float128 C1 = L(6.93145751953125E-1);
117static const _Float128 C2 = L(1.428606820309417232121458176568075500134E-6);
118
119static const _Float128 sqrth = L(0.7071067811865475244008443621048490392848);
120/* ln (2^16384 * (1 - 2^-113)) */
121static const _Float128 zero = 0;
122
123_Float128
124__log1pl (_Float128 xm1)
125{
126 _Float128 x, y, z, r, s;
127 ieee854_long_double_shape_type u;
128 int32_t hx;
129 int e;
130
131 /* Test for NaN or infinity input. */
132 u.value = xm1;
133 hx = u.parts32.w0;
134 if ((hx & 0x7fffffff) >= 0x7fff0000)
135 return xm1 + fabsl (xm1);
136
137 /* log1p(+- 0) = +- 0. */
138 if (((hx & 0x7fffffff) == 0)
139 && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
140 return xm1;
141
142 if ((hx & 0x7fffffff) < 0x3f8e0000)
143 {
144 math_check_force_underflow (xm1);
145 if ((int) xm1 == 0)
146 return xm1;
147 }
148
149 if (xm1 >= L(0x1p113))
150 x = xm1;
151 else
152 x = xm1 + 1;
153
154 /* log1p(-1) = -inf */
155 if (x <= 0)
156 {
157 if (x == 0)
158 return (-1 / zero); /* log1p(-1) = -inf */
159 else
160 return (zero / (x - x));
161 }
162
163 /* Separate mantissa from exponent. */
164
165 /* Use frexp used so that denormal numbers will be handled properly. */
166 x = __frexpl (x, &e);
167
168 /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
169 where z = 2(x-1)/x+1). */
170 if ((e > 2) || (e < -2))
171 {
172 if (x < sqrth)
173 { /* 2( 2x-1 )/( 2x+1 ) */
174 e -= 1;
175 z = x - L(0.5);
176 y = L(0.5) * z + L(0.5);
177 }
178 else
179 { /* 2 (x-1)/(x+1) */
180 z = x - L(0.5);
181 z -= L(0.5);
182 y = L(0.5) * x + L(0.5);
183 }
184 x = z / y;
185 z = x * x;
186 r = ((((R5 * z
187 + R4) * z
188 + R3) * z
189 + R2) * z
190 + R1) * z
191 + R0;
192 s = (((((z
193 + S5) * z
194 + S4) * z
195 + S3) * z
196 + S2) * z
197 + S1) * z
198 + S0;
199 z = x * (z * r / s);
200 z = z + e * C2;
201 z = z + x;
202 z = z + e * C1;
203 return (z);
204 }
205
206
207 /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
208
209 if (x < sqrth)
210 {
211 e -= 1;
212 if (e != 0)
213 x = 2 * x - 1; /* 2x - 1 */
214 else
215 x = xm1;
216 }
217 else
218 {
219 if (e != 0)
220 x = x - 1;
221 else
222 x = xm1;
223 }
224 z = x * x;
225 r = (((((((((((P12 * x
226 + P11) * x
227 + P10) * x
228 + P9) * x
229 + P8) * x
230 + P7) * x
231 + P6) * x
232 + P5) * x
233 + P4) * x
234 + P3) * x
235 + P2) * x
236 + P1) * x
237 + P0;
238 s = (((((((((((x
239 + Q11) * x
240 + Q10) * x
241 + Q9) * x
242 + Q8) * x
243 + Q7) * x
244 + Q6) * x
245 + Q5) * x
246 + Q4) * x
247 + Q3) * x
248 + Q2) * x
249 + Q1) * x
250 + Q0;
251 y = x * (z * r / s);
252 y = y + e * C2;
253 z = y - L(0.5) * z;
254 z = z + x;
255 z = z + e * C1;
256 return (z);
257}
258