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