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