2 * Base 2 logarithm, 128-bit long double precision
8 * long double x, y, log2l();
16 * Returns the base 2 logarithm of x.
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
22 * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
24 * Otherwise, setting z = 2(x-1)/x+1),
26 * log(x) = z + z^3 P(z)/Q(z).
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
37 * In the tests over the interval exp(+-10000), the logarithms
38 * of the random arguments were uniformly distributed over
44 Cephes Math Library Release 2.2: January, 1991
45 Copyright 1984, 1991 by Stephen L. Moshier
46 Adapted for glibc November, 2001
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.
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.
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/>.
63 #include <math_private.h>
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
70 static const long double P
[13] =
72 1.313572404063446165910279910527789794488E4L
,
73 7.771154681358524243729929227226708890930E4L
,
74 2.014652742082537582487669938141683759923E5L
,
75 3.007007295140399532324943111654767187848E5L
,
76 2.854829159639697837788887080758954924001E5L
,
77 1.797628303815655343403735250238293741397E5L
,
78 7.594356839258970405033155585486712125861E4L
,
79 2.128857716871515081352991964243375186031E4L
,
80 3.824952356185897735160588078446136783779E3L
,
81 4.114517881637811823002128927449878962058E2L
,
82 2.321125933898420063925789532045674660756E1L
,
83 4.998469661968096229986658302195402690910E-1L,
84 1.538612243596254322971797716843006400388E-6L
86 static const long double Q
[12] =
88 3.940717212190338497730839731583397586124E4L
,
89 2.626900195321832660448791748036714883242E5L
,
90 7.777690340007566932935753241556479363645E5L
,
91 1.347518538384329112529391120390701166528E6L
,
92 1.514882452993549494932585972882995548426E6L
,
93 1.158019977462989115839826904108208787040E6L
,
94 6.132189329546557743179177159925690841200E5L
,
95 2.248234257620569139969141618556349415120E5L
,
96 5.605842085972455027590989944010492125825E4L
,
97 9.147150349299596453976674231612674085381E3L
,
98 9.104928120962988414618126155557301584078E2L
,
99 4.839208193348159620282142911143429644326E1L
100 /* 1.000000000000000000000000000000000000000E0L, */
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
109 static const long double R
[6] =
111 1.418134209872192732479751274970992665513E5L
,
112 -8.977257995689735303686582344659576526998E4L
,
113 2.048819892795278657810231591630928516206E4L
,
114 -2.024301798136027039250415126250455056397E3L
,
115 8.057002716646055371965756206836056074715E1L
,
116 -8.828896441624934385266096344596648080902E-1L
118 static const long double S
[6] =
120 1.701761051846631278975701529965589676574E6L
,
121 -1.332535117259762928288745111081235577029E6L
,
122 4.001557694070773974936904547424676279307E5L
,
123 -5.748542087379434595104154610899551484314E4L
,
124 3.998526750980007367835804959888064681098E3L
,
125 -1.186359407982897997337150403816839480438E2L
126 /* 1.000000000000000000000000000000000000000E0L, */
129 static const long double
131 LOG2EA
= 4.4269504088896340735992468100189213742664595E-1L,
133 SQRTH
= 7.071067811865475244008443621048490392848359E-1L;
136 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
139 neval (long double x
, const long double *p
, int n
)
154 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
157 deval (long double x
, const long double *p
, int n
)
182 /* Test for domain */
183 GET_LDOUBLE_WORDS64 (hx
, lx
, x
);
184 if (((hx
& 0x7fffffffffffffffLL
) | lx
) == 0)
185 return (-1.0L / (x
- x
));
187 return (x
- x
) / (x
- x
);
188 if (hx
>= 0x7fff000000000000LL
)
194 /* separate mantissa from exponent */
196 /* Note, frexp is used so that denormal numbers
197 * will be handled properly.
199 x
= __frexpl (x
, &e
);
202 /* logarithm using log(x) = z + z**3 P(z)/Q(z),
203 * where z = 2(x-1)/x+1)
205 if ((e
> 2) || (e
< -2))
208 { /* 2( 2x-1 )/( 2x+1 ) */
214 { /* 2 (x-1)/(x+1) */
221 y
= x
* (z
* neval (z
, R
, 5) / deval (z
, S
, 5));
226 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
231 x
= 2.0 * x
- 1.0L; /* 2x - 1 */
238 y
= x
* (z
* neval (x
, P
, 12) / deval (x
, Q
, 11));
243 /* Multiply log of fraction by log2(e)
244 * and base 2 exponent by 1
253 strong_alias (__ieee754_log2l
, __log2l_finite
)