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