* sysdeps/ieee754/ldbl-128/e_j0l.c (__ieee751_j0l): Use __finitel.
[glibc.git] / sysdeps / ieee754 / ldbl-128 / s_log1pl.c
blobdefbe76f5e72e9060eb35bcae9ed4b5592c85ca3
1 /* log1pl.c
3 * Relative error logarithm
4 * Natural logarithm of 1+x, 128-bit long double precision
8 * SYNOPSIS:
10 * long double x, y, log1pl();
12 * y = log1pl( x );
16 * DESCRIPTION:
18 * Returns the base e (2.718...) logarithm of 1+x.
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
24 * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
26 * Otherwise, setting z = 2(w-1)/(w+1),
28 * log(w) = z + z^3 P(z)/Q(z).
32 * ACCURACY:
34 * Relative error:
35 * arithmetic domain # trials peak rms
36 * IEEE -1, 8 100000 1.9e-34 4.3e-35
39 /* Copyright 2001 by Stephen L. Moshier
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.
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.
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 */
56 #include "math.h"
57 #include "math_private.h"
59 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
60 * 1/sqrt(2) <= 1+x < sqrt(2)
61 * Theoretical peak relative error = 5.3e-37,
62 * relative peak error spread = 2.3e-14
64 static const long double
65 P12 = 1.538612243596254322971797716843006400388E-6L,
66 P11 = 4.998469661968096229986658302195402690910E-1L,
67 P10 = 2.321125933898420063925789532045674660756E1L,
68 P9 = 4.114517881637811823002128927449878962058E2L,
69 P8 = 3.824952356185897735160588078446136783779E3L,
70 P7 = 2.128857716871515081352991964243375186031E4L,
71 P6 = 7.594356839258970405033155585486712125861E4L,
72 P5 = 1.797628303815655343403735250238293741397E5L,
73 P4 = 2.854829159639697837788887080758954924001E5L,
74 P3 = 3.007007295140399532324943111654767187848E5L,
75 P2 = 2.014652742082537582487669938141683759923E5L,
76 P1 = 7.771154681358524243729929227226708890930E4L,
77 P0 = 1.313572404063446165910279910527789794488E4L,
78 /* Q12 = 1.000000000000000000000000000000000000000E0L, */
79 Q11 = 4.839208193348159620282142911143429644326E1L,
80 Q10 = 9.104928120962988414618126155557301584078E2L,
81 Q9 = 9.147150349299596453976674231612674085381E3L,
82 Q8 = 5.605842085972455027590989944010492125825E4L,
83 Q7 = 2.248234257620569139969141618556349415120E5L,
84 Q6 = 6.132189329546557743179177159925690841200E5L,
85 Q5 = 1.158019977462989115839826904108208787040E6L,
86 Q4 = 1.514882452993549494932585972882995548426E6L,
87 Q3 = 1.347518538384329112529391120390701166528E6L,
88 Q2 = 7.777690340007566932935753241556479363645E5L,
89 Q1 = 2.626900195321832660448791748036714883242E5L,
90 Q0 = 3.940717212190338497730839731583397586124E4L;
92 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
93 * where z = 2(x-1)/(x+1)
94 * 1/sqrt(2) <= x < sqrt(2)
95 * Theoretical peak relative error = 1.1e-35,
96 * relative peak error spread 1.1e-9
98 static const long double
99 R5 = -8.828896441624934385266096344596648080902E-1L,
100 R4 = 8.057002716646055371965756206836056074715E1L,
101 R3 = -2.024301798136027039250415126250455056397E3L,
102 R2 = 2.048819892795278657810231591630928516206E4L,
103 R1 = -8.977257995689735303686582344659576526998E4L,
104 R0 = 1.418134209872192732479751274970992665513E5L,
105 /* S6 = 1.000000000000000000000000000000000000000E0L, */
106 S5 = -1.186359407982897997337150403816839480438E2L,
107 S4 = 3.998526750980007367835804959888064681098E3L,
108 S3 = -5.748542087379434595104154610899551484314E4L,
109 S2 = 4.001557694070773974936904547424676279307E5L,
110 S1 = -1.332535117259762928288745111081235577029E6L,
111 S0 = 1.701761051846631278975701529965589676574E6L;
113 /* C1 + C2 = ln 2 */
114 static const long double C1 = 6.93145751953125E-1L;
115 static const long double C2 = 1.428606820309417232121458176568075500134E-6L;
117 static const long double sqrth = 0.7071067811865475244008443621048490392848L;
118 /* ln (2^16384 * (1 - 2^-113)) */
119 static const long double maxlog = 1.1356523406294143949491931077970764891253E4L;
120 static const long double big = 2e4932L;
121 static const long double zero = 0.0L;
123 long double
124 __log1pl (long double xm1)
126 long double x, y, z, r, s;
127 ieee854_long_double_shape_type u;
128 int32_t hx;
129 int e;
131 /* Test for NaN or infinity input. */
132 u.value = xm1;
133 hx = u.parts32.w0;
134 if (hx >= 0x7fff0000)
135 return xm1;
137 /* log1p(+- 0) = +- 0. */
138 if (((hx & 0x7fffffff) == 0)
139 && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
140 return xm1;
142 x = xm1 + 1.0L;
144 /* log1p(-1) = -inf */
145 if (x <= 0.0L)
147 if (x == 0.0L)
148 return (-1.0L / (x - x));
149 else
150 return (zero / (x - x));
153 /* Separate mantissa from exponent. */
155 /* Use frexp used so that denormal numbers will be handled properly. */
156 x = __frexpl (x, &e);
158 /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
159 where z = 2(x-1)/x+1). */
160 if ((e > 2) || (e < -2))
162 if (x < sqrth)
163 { /* 2( 2x-1 )/( 2x+1 ) */
164 e -= 1;
165 z = x - 0.5L;
166 y = 0.5L * z + 0.5L;
168 else
169 { /* 2 (x-1)/(x+1) */
170 z = x - 0.5L;
171 z -= 0.5L;
172 y = 0.5L * x + 0.5L;
174 x = z / y;
175 z = x * x;
176 r = ((((R5 * z
177 + R4) * z
178 + R3) * z
179 + R2) * z
180 + R1) * z
181 + R0;
182 s = (((((z
183 + S5) * z
184 + S4) * z
185 + S3) * z
186 + S2) * z
187 + S1) * z
188 + S0;
189 z = x * (z * r / s);
190 z = z + e * C2;
191 z = z + x;
192 z = z + e * C1;
193 return (z);
197 /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
199 if (x < sqrth)
201 e -= 1;
202 if (e != 0)
203 x = 2.0L * x - 1.0L; /* 2x - 1 */
204 else
205 x = xm1;
207 else
209 if (e != 0)
210 x = x - 1.0L;
211 else
212 x = xm1;
214 z = x * x;
215 r = (((((((((((P12 * x
216 + P11) * x
217 + P10) * x
218 + P9) * x
219 + P8) * x
220 + P7) * x
221 + P6) * x
222 + P5) * x
223 + P4) * x
224 + P3) * x
225 + P2) * x
226 + P1) * x
227 + P0;
228 s = (((((((((((x
229 + Q11) * x
230 + Q10) * x
231 + Q9) * x
232 + Q8) * x
233 + Q7) * x
234 + Q6) * x
235 + Q5) * x
236 + Q4) * x
237 + Q3) * x
238 + Q2) * x
239 + Q1) * x
240 + Q0;
241 y = x * (z * r / s);
242 y = y + e * C2;
243 z = y - 0.5L * z;
244 z = z + x;
245 z = z + e * C1;
246 return (z);
249 weak_alias (__log1pl, log1pl)