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 <https://www.gnu.org/licenses/>.
61 */
62
63 #include <math.h>
64 #include <math_private.h>
65 #include <libm-alias-finite.h>
66
67 /* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
68 * 1/sqrt(2) <= x < sqrt(2)
69 * Theoretical peak relative error = 5.3e-37,
70 * relative peak error spread = 2.3e-14
71 */
72 static const _Float128 P[13] =
73 {
74 L(1.313572404063446165910279910527789794488E4),
75 L(7.771154681358524243729929227226708890930E4),
76 L(2.014652742082537582487669938141683759923E5),
77 L(3.007007295140399532324943111654767187848E5),
78 L(2.854829159639697837788887080758954924001E5),
79 L(1.797628303815655343403735250238293741397E5),
80 L(7.594356839258970405033155585486712125861E4),
81 L(2.128857716871515081352991964243375186031E4),
82 L(3.824952356185897735160588078446136783779E3),
83 L(4.114517881637811823002128927449878962058E2),
84 L(2.321125933898420063925789532045674660756E1),
85 L(4.998469661968096229986658302195402690910E-1),
86 L(1.538612243596254322971797716843006400388E-6)
87 };
88 static const _Float128 Q[12] =
89 {
90 L(3.940717212190338497730839731583397586124E4),
91 L(2.626900195321832660448791748036714883242E5),
92 L(7.777690340007566932935753241556479363645E5),
93 L(1.347518538384329112529391120390701166528E6),
94 L(1.514882452993549494932585972882995548426E6),
95 L(1.158019977462989115839826904108208787040E6),
96 L(6.132189329546557743179177159925690841200E5),
97 L(2.248234257620569139969141618556349415120E5),
98 L(5.605842085972455027590989944010492125825E4),
99 L(9.147150349299596453976674231612674085381E3),
100 L(9.104928120962988414618126155557301584078E2),
101 L(4.839208193348159620282142911143429644326E1)
102 /* 1.000000000000000000000000000000000000000E0L, */
103 };
104
105 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
106 * where z = 2(x-1)/(x+1)
107 * 1/sqrt(2) <= x < sqrt(2)
108 * Theoretical peak relative error = 1.1e-35,
109 * relative peak error spread 1.1e-9
110 */
111 static const _Float128 R[6] =
112 {
113 L(1.418134209872192732479751274970992665513E5),
114 L(-8.977257995689735303686582344659576526998E4),
115 L(2.048819892795278657810231591630928516206E4),
116 L(-2.024301798136027039250415126250455056397E3),
117 L(8.057002716646055371965756206836056074715E1),
118 L(-8.828896441624934385266096344596648080902E-1)
119 };
120 static const _Float128 S[6] =
121 {
122 L(1.701761051846631278975701529965589676574E6),
123 L(-1.332535117259762928288745111081235577029E6),
124 L(4.001557694070773974936904547424676279307E5),
125 L(-5.748542087379434595104154610899551484314E4),
126 L(3.998526750980007367835804959888064681098E3),
127 L(-1.186359407982897997337150403816839480438E2)
128 /* 1.000000000000000000000000000000000000000E0L, */
129 };
130
131 static const _Float128
132 /* log10(2) */
133 L102A = L(0.3125),
134 L102B = L(-1.14700043360188047862611052755069732318101185E-2),
135 /* log10(e) */
136 L10EA = L(0.5),
137 L10EB = L(-6.570551809674817234887108108339491770560299E-2),
138 /* sqrt(2)/2 */
139 SQRTH = L(7.071067811865475244008443621048490392848359E-1);
140
141
142
143 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
144
145 static _Float128
neval(_Float128 x,const _Float128 * p,int n)146 neval (_Float128 x, const _Float128 *p, int n)
147 {
148 _Float128 y;
149
150 p += n;
151 y = *p--;
152 do
153 {
154 y = y * x + *p--;
155 }
156 while (--n > 0);
157 return y;
158 }
159
160
161 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
162
163 static _Float128
deval(_Float128 x,const _Float128 * p,int n)164 deval (_Float128 x, const _Float128 *p, int n)
165 {
166 _Float128 y;
167
168 p += n;
169 y = x + *p--;
170 do
171 {
172 y = y * x + *p--;
173 }
174 while (--n > 0);
175 return y;
176 }
177
178
179
180 _Float128
__ieee754_log10l(_Float128 x)181 __ieee754_log10l (_Float128 x)
182 {
183 _Float128 z;
184 _Float128 y;
185 int e;
186 int64_t hx, lx;
187
188 /* Test for domain */
189 GET_LDOUBLE_WORDS64 (hx, lx, x);
190 if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
191 return (-1 / fabsl (x)); /* log10l(+-0)=-inf */
192 if (hx < 0)
193 return (x - x) / (x - x);
194 if (hx >= 0x7fff000000000000LL)
195 return (x + x);
196
197 if (x == 1)
198 return 0;
199
200 /* separate mantissa from exponent */
201
202 /* Note, frexp is used so that denormal numbers
203 * will be handled properly.
204 */
205 x = __frexpl (x, &e);
206
207
208 /* logarithm using log(x) = z + z**3 P(z)/Q(z),
209 * where z = 2(x-1)/x+1)
210 */
211 if ((e > 2) || (e < -2))
212 {
213 if (x < SQRTH)
214 { /* 2( 2x-1 )/( 2x+1 ) */
215 e -= 1;
216 z = x - L(0.5);
217 y = L(0.5) * z + L(0.5);
218 }
219 else
220 { /* 2 (x-1)/(x+1) */
221 z = x - L(0.5);
222 z -= L(0.5);
223 y = L(0.5) * x + L(0.5);
224 }
225 x = z / y;
226 z = x * x;
227 y = x * (z * neval (z, R, 5) / deval (z, S, 5));
228 goto done;
229 }
230
231
232 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
233
234 if (x < SQRTH)
235 {
236 e -= 1;
237 x = 2.0 * x - 1; /* 2x - 1 */
238 }
239 else
240 {
241 x = x - 1;
242 }
243 z = x * x;
244 y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
245 y = y - 0.5 * z;
246
247 done:
248
249 /* Multiply log of fraction by log10(e)
250 * and base 2 exponent by log10(2).
251 */
252 z = y * L10EB;
253 z += x * L10EB;
254 z += e * L102B;
255 z += y * L10EA;
256 z += x * L10EA;
257 z += e * L102A;
258 return (z);
259 }
260 libm_alias_finite (__ieee754_log10l, __log10l)
261