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 long double P[13] =
73 {
74 1.313572404063446165910279910527789794488E4L,
75 7.771154681358524243729929227226708890930E4L,
76 2.014652742082537582487669938141683759923E5L,
77 3.007007295140399532324943111654767187848E5L,
78 2.854829159639697837788887080758954924001E5L,
79 1.797628303815655343403735250238293741397E5L,
80 7.594356839258970405033155585486712125861E4L,
81 2.128857716871515081352991964243375186031E4L,
82 3.824952356185897735160588078446136783779E3L,
83 4.114517881637811823002128927449878962058E2L,
84 2.321125933898420063925789532045674660756E1L,
85 4.998469661968096229986658302195402690910E-1L,
86 1.538612243596254322971797716843006400388E-6L
87 };
88 static const long double Q[12] =
89 {
90 3.940717212190338497730839731583397586124E4L,
91 2.626900195321832660448791748036714883242E5L,
92 7.777690340007566932935753241556479363645E5L,
93 1.347518538384329112529391120390701166528E6L,
94 1.514882452993549494932585972882995548426E6L,
95 1.158019977462989115839826904108208787040E6L,
96 6.132189329546557743179177159925690841200E5L,
97 2.248234257620569139969141618556349415120E5L,
98 5.605842085972455027590989944010492125825E4L,
99 9.147150349299596453976674231612674085381E3L,
100 9.104928120962988414618126155557301584078E2L,
101 4.839208193348159620282142911143429644326E1L
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 long double R[6] =
112 {
113 1.418134209872192732479751274970992665513E5L,
114 -8.977257995689735303686582344659576526998E4L,
115 2.048819892795278657810231591630928516206E4L,
116 -2.024301798136027039250415126250455056397E3L,
117 8.057002716646055371965756206836056074715E1L,
118 -8.828896441624934385266096344596648080902E-1L
119 };
120 static const long double S[6] =
121 {
122 1.701761051846631278975701529965589676574E6L,
123 -1.332535117259762928288745111081235577029E6L,
124 4.001557694070773974936904547424676279307E5L,
125 -5.748542087379434595104154610899551484314E4L,
126 3.998526750980007367835804959888064681098E3L,
127 -1.186359407982897997337150403816839480438E2L
128 /* 1.000000000000000000000000000000000000000E0L, */
129 };
130
131 static const long double
132 /* log10(2) */
133 L102A = 0.3125L,
134 L102B = -1.14700043360188047862611052755069732318101185E-2L,
135 /* log10(e) */
136 L10EA = 0.5L,
137 L10EB = -6.570551809674817234887108108339491770560299E-2L,
138 /* sqrt(2)/2 */
139 SQRTH = 7.071067811865475244008443621048490392848359E-1L;
140
141
142
143 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
144
145 static long double
146 neval (long double x, const long double *p, int n)
147 {
148 long double 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 long double
164 deval (long double x, const long double *p, int n)
165 {
166 long double 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 long double
181 __ieee754_log10l (long double x)
182 {
183 long double z;
184 long double y;
185 int e;
186 int64_t hx;
187 double xhi;
188
189 /* Test for domain */
190 xhi = ldbl_high (x);
191 EXTRACT_WORDS64 (hx, xhi);
192 if ((hx & 0x7fffffffffffffffLL) == 0)
193 return (-1.0L / fabsl (x)); /* log10l(+-0)=-inf */
194 if (hx < 0)
195 return (x - x) / (x - x);
196 if (hx >= 0x7ff0000000000000LL)
197 return (x + x);
198
199 if (x == 1.0L)
200 return 0.0L;
201
202 /* separate mantissa from exponent */
203
204 /* Note, frexp is used so that denormal numbers
205 * will be handled properly.
206 */
207 x = __frexpl (x, &e);
208
209
210 /* logarithm using log(x) = z + z**3 P(z)/Q(z),
211 * where z = 2(x-1)/x+1)
212 */
213 if ((e > 2) || (e < -2))
214 {
215 if (x < SQRTH)
216 { /* 2( 2x-1 )/( 2x+1 ) */
217 e -= 1;
218 z = x - 0.5L;
219 y = 0.5L * z + 0.5L;
220 }
221 else
222 { /* 2 (x-1)/(x+1) */
223 z = x - 0.5L;
224 z -= 0.5L;
225 y = 0.5L * x + 0.5L;
226 }
227 x = z / y;
228 z = x * x;
229 y = x * (z * neval (z, R, 5) / deval (z, S, 5));
230 goto done;
231 }
232
233
234 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
235
236 if (x < SQRTH)
237 {
238 e -= 1;
239 x = 2.0 * x - 1.0L; /* 2x - 1 */
240 }
241 else
242 {
243 x = x - 1.0L;
244 }
245 z = x * x;
246 y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
247 y = y - 0.5 * z;
248
249 done:
250
251 /* Multiply log of fraction by log10(e)
252 * and base 2 exponent by log10(2).
253 */
254 z = y * L10EB;
255 z += x * L10EB;
256 z += e * L102B;
257 z += y * L10EA;
258 z += x * L10EA;
259 z += e * L102A;
260 return (z);
261 }
262 libm_alias_finite (__ieee754_log10l, __log10l)