1/* log10l.c
2 *
3 * Common logarithm, 128-bit long double precision
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * long double x, y, log10l();
10 *
11 * y = log10l( x );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Returns the base 10 logarithm of x.
18 *
19 * The argument is separated into its exponent and fractional
20 * parts. If the exponent is between -1 and +1, the logarithm
21 * of the fraction is approximated by
22 *
23 * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
24 *
25 * Otherwise, setting z = 2(x-1)/x+1),
26 *
27 * log(x) = z + z^3 P(z)/Q(z).
28 *
29 *
30 *
31 * ACCURACY:
32 *
33 * Relative error:
34 * arithmetic domain # trials peak rms
35 * IEEE 0.5, 2.0 30000 2.3e-34 4.9e-35
36 * IEEE exp(+-10000) 30000 1.0e-34 4.1e-35
37 *
38 * In the tests over the interval exp(+-10000), the logarithms
39 * of the random arguments were uniformly distributed over
40 * [-10000, +10000].
41 *
42 */
43
44/*
45 Cephes Math Library Release 2.2: January, 1991
46 Copyright 1984, 1991 by Stephen L. Moshier
47 Adapted for glibc November, 2001
48
49 This library is free software; you can redistribute it and/or
50 modify it under the terms of the GNU Lesser General Public
51 License as published by the Free Software Foundation; either
52 version 2.1 of the License, or (at your option) any later version.
53
54 This library is distributed in the hope that it will be useful,
55 but WITHOUT ANY WARRANTY; without even the implied warranty of
56 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
57 Lesser General Public License for more details.
58
59 You should have received a copy of the GNU Lesser General Public
60 License along with this library; if not, see <http://www.gnu.org/licenses/>.
61 */
62
63#include <math.h>
64#include <math_private.h>
65
66/* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
67 * 1/sqrt(2) <= x < sqrt(2)
68 * Theoretical peak relative error = 5.3e-37,
69 * relative peak error spread = 2.3e-14
70 */
71static const _Float128 P[13] =
72{
73 L(1.313572404063446165910279910527789794488E4),
74 L(7.771154681358524243729929227226708890930E4),
75 L(2.014652742082537582487669938141683759923E5),
76 L(3.007007295140399532324943111654767187848E5),
77 L(2.854829159639697837788887080758954924001E5),
78 L(1.797628303815655343403735250238293741397E5),
79 L(7.594356839258970405033155585486712125861E4),
80 L(2.128857716871515081352991964243375186031E4),
81 L(3.824952356185897735160588078446136783779E3),
82 L(4.114517881637811823002128927449878962058E2),
83 L(2.321125933898420063925789532045674660756E1),
84 L(4.998469661968096229986658302195402690910E-1),
85 L(1.538612243596254322971797716843006400388E-6)
86};
87static const _Float128 Q[12] =
88{
89 L(3.940717212190338497730839731583397586124E4),
90 L(2.626900195321832660448791748036714883242E5),
91 L(7.777690340007566932935753241556479363645E5),
92 L(1.347518538384329112529391120390701166528E6),
93 L(1.514882452993549494932585972882995548426E6),
94 L(1.158019977462989115839826904108208787040E6),
95 L(6.132189329546557743179177159925690841200E5),
96 L(2.248234257620569139969141618556349415120E5),
97 L(5.605842085972455027590989944010492125825E4),
98 L(9.147150349299596453976674231612674085381E3),
99 L(9.104928120962988414618126155557301584078E2),
100 L(4.839208193348159620282142911143429644326E1)
101/* 1.000000000000000000000000000000000000000E0L, */
102};
103
104/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
105 * where z = 2(x-1)/(x+1)
106 * 1/sqrt(2) <= x < sqrt(2)
107 * Theoretical peak relative error = 1.1e-35,
108 * relative peak error spread 1.1e-9
109 */
110static const _Float128 R[6] =
111{
112 L(1.418134209872192732479751274970992665513E5),
113 L(-8.977257995689735303686582344659576526998E4),
114 L(2.048819892795278657810231591630928516206E4),
115 L(-2.024301798136027039250415126250455056397E3),
116 L(8.057002716646055371965756206836056074715E1),
117 L(-8.828896441624934385266096344596648080902E-1)
118};
119static const _Float128 S[6] =
120{
121 L(1.701761051846631278975701529965589676574E6),
122 L(-1.332535117259762928288745111081235577029E6),
123 L(4.001557694070773974936904547424676279307E5),
124 L(-5.748542087379434595104154610899551484314E4),
125 L(3.998526750980007367835804959888064681098E3),
126 L(-1.186359407982897997337150403816839480438E2)
127/* 1.000000000000000000000000000000000000000E0L, */
128};
129
130static const _Float128
131/* log10(2) */
132L102A = L(0.3125),
133L102B = L(-1.14700043360188047862611052755069732318101185E-2),
134/* log10(e) */
135L10EA = L(0.5),
136L10EB = L(-6.570551809674817234887108108339491770560299E-2),
137/* sqrt(2)/2 */
138SQRTH = L(7.071067811865475244008443621048490392848359E-1);
139
140
141
142/* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
143
144static _Float128
145neval (_Float128 x, const _Float128 *p, int n)
146{
147 _Float128 y;
148
149 p += n;
150 y = *p--;
151 do
152 {
153 y = y * x + *p--;
154 }
155 while (--n > 0);
156 return y;
157}
158
159
160/* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
161
162static _Float128
163deval (_Float128 x, const _Float128 *p, int n)
164{
165 _Float128 y;
166
167 p += n;
168 y = x + *p--;
169 do
170 {
171 y = y * x + *p--;
172 }
173 while (--n > 0);
174 return y;
175}
176
177
178
179_Float128
180__ieee754_log10l (_Float128 x)
181{
182 _Float128 z;
183 _Float128 y;
184 int e;
185 int64_t hx, lx;
186
187/* Test for domain */
188 GET_LDOUBLE_WORDS64 (hx, lx, x);
189 if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
190 return (-1 / fabsl (x)); /* log10l(+-0)=-inf */
191 if (hx < 0)
192 return (x - x) / (x - x);
193 if (hx >= 0x7fff000000000000LL)
194 return (x + x);
195
196 if (x == 1)
197 return 0;
198
199/* separate mantissa from exponent */
200
201/* Note, frexp is used so that denormal numbers
202 * will be handled properly.
203 */
204 x = __frexpl (x, &e);
205
206
207/* logarithm using log(x) = z + z**3 P(z)/Q(z),
208 * where z = 2(x-1)/x+1)
209 */
210 if ((e > 2) || (e < -2))
211 {
212 if (x < SQRTH)
213 { /* 2( 2x-1 )/( 2x+1 ) */
214 e -= 1;
215 z = x - L(0.5);
216 y = L(0.5) * z + L(0.5);
217 }
218 else
219 { /* 2 (x-1)/(x+1) */
220 z = x - L(0.5);
221 z -= L(0.5);
222 y = L(0.5) * x + L(0.5);
223 }
224 x = z / y;
225 z = x * x;
226 y = x * (z * neval (z, R, 5) / deval (z, S, 5));
227 goto done;
228 }
229
230
231/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
232
233 if (x < SQRTH)
234 {
235 e -= 1;
236 x = 2.0 * x - 1; /* 2x - 1 */
237 }
238 else
239 {
240 x = x - 1;
241 }
242 z = x * x;
243 y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
244 y = y - 0.5 * z;
245
246done:
247
248 /* Multiply log of fraction by log10(e)
249 * and base 2 exponent by log10(2).
250 */
251 z = y * L10EB;
252 z += x * L10EB;
253 z += e * L102B;
254 z += y * L10EA;
255 z += x * L10EA;
256 z += e * L102A;
257 return (z);
258}
259strong_alias (__ieee754_log10l, __log10l_finite)
260