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