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 */
63 static 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 */
97 static 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 */
113 static const __float128 C1 = 6.93145751953125E-1Q;
114 static const __float128 C2 = 1.428606820309417232121458176568075500134E-6Q;
115
116 static const __float128 sqrth = 0.7071067811865475244008443621048490392848Q;
117 static const __float128 zero = 0.0Q;
118
119
120 __float128
log1pq(__float128 xm1)121 log1pq (__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 >= 0x7fff0000)
132 return 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 if ((int) xm1 == 0)
142 return xm1;
143 }
144
145 x = xm1 + 1.0Q;
146
147 /* log1p(-1) = -inf */
148 if (x <= 0.0Q)
149 {
150 if (x == 0.0Q)
151 return (-1.0Q / (x - x));
152 else
153 return (zero / (x - x));
154 }
155
156 /* Separate mantissa from exponent. */
157
158 /* Use frexp used so that denormal numbers will be handled properly. */
159 x = frexpq (x, &e);
160
161 /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
162 where z = 2(x-1)/x+1). */
163 if ((e > 2) || (e < -2))
164 {
165 if (x < sqrth)
166 { /* 2( 2x-1 )/( 2x+1 ) */
167 e -= 1;
168 z = x - 0.5Q;
169 y = 0.5Q * z + 0.5Q;
170 }
171 else
172 { /* 2 (x-1)/(x+1) */
173 z = x - 0.5Q;
174 z -= 0.5Q;
175 y = 0.5Q * x + 0.5Q;
176 }
177 x = z / y;
178 z = x * x;
179 r = ((((R5 * z
180 + R4) * z
181 + R3) * z
182 + R2) * z
183 + R1) * z
184 + R0;
185 s = (((((z
186 + S5) * z
187 + S4) * z
188 + S3) * z
189 + S2) * z
190 + S1) * z
191 + S0;
192 z = x * (z * r / s);
193 z = z + e * C2;
194 z = z + x;
195 z = z + e * C1;
196 return (z);
197 }
198
199
200 /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
201
202 if (x < sqrth)
203 {
204 e -= 1;
205 if (e != 0)
206 x = 2.0Q * x - 1.0Q; /* 2x - 1 */
207 else
208 x = xm1;
209 }
210 else
211 {
212 if (e != 0)
213 x = x - 1.0Q;
214 else
215 x = xm1;
216 }
217 z = x * x;
218 r = (((((((((((P12 * x
219 + P11) * x
220 + P10) * x
221 + P9) * x
222 + P8) * x
223 + P7) * x
224 + P6) * x
225 + P5) * x
226 + P4) * x
227 + P3) * x
228 + P2) * x
229 + P1) * x
230 + P0;
231 s = (((((((((((x
232 + Q11) * x
233 + Q10) * x
234 + Q9) * x
235 + Q8) * x
236 + Q7) * x
237 + Q6) * x
238 + Q5) * x
239 + Q4) * x
240 + Q3) * x
241 + Q2) * x
242 + Q1) * x
243 + Q0;
244 y = x * (z * r / s);
245 y = y + e * C2;
246 z = y - 0.5Q * z;
247 z = z + x;
248 z = z + e * C1;
249 return (z);
250 }
251