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