]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/ieee754/dbl-64/wordsize-64/e_log2.c
6dc7b7d217266a6a3826738482e9ef652b10ac14
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
9 * ====================================================
13 * Return the logarithm to base 2 of x
16 * 1. Argument Reduction: find k and f such that
18 * where sqrt(2)/2 < 1+f < sqrt(2) .
20 * 2. Approximation of log(1+f).
21 * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
22 * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
24 * We use a special Reme algorithm on [0,0.1716] to generate
25 * a polynomial of degree 14 to approximate R The maximum error
26 * of this polynomial approximation is bounded by 2**-58.45. In
29 * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
30 * (the values of Lg1 to Lg7 are listed in the program)
33 * | Lg1*s +...+Lg7*s - R(z) | <= 2
35 * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
36 * In order to guarantee error in log below 1ulp, we compute log
38 * log(1+f) = f - s*(f - R) (if f is not too large)
39 * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
41 * 3. Finally, log(x) = k + log(1+f).
42 * = k+(f-(hfsq-(s*(hfsq+R))))
45 * log2(x) is NaN with signal if x < 0 (including -INF) ;
46 * log2(+INF) is +INF; log(0) is -INF with signal;
47 * log2(NaN) is that NaN with no signal.
50 * The hexadecimal values are the intended ones for the following
51 * constants. The decimal values may be used, provided that the
52 * compiler will convert from decimal to binary accurately enough
53 * to produce the hexadecimal values shown.
57 #include <math_private.h>
59 static const double ln2
= 0.69314718055994530942;
60 static const double two54
= 1.80143985094819840000e+16; /* 4350000000000000 */
61 static const double Lg1
= 6.666666666666735130e-01; /* 3FE5555555555593 */
62 static const double Lg2
= 3.999999999940941908e-01; /* 3FD999999997FA04 */
63 static const double Lg3
= 2.857142874366239149e-01; /* 3FD2492494229359 */
64 static const double Lg4
= 2.222219843214978396e-01; /* 3FCC71C51D8E78AF */
65 static const double Lg5
= 1.818357216161805012e-01; /* 3FC7466496CB03DE */
66 static const double Lg6
= 1.531383769920937332e-01; /* 3FC39A09D078C69F */
67 static const double Lg7
= 1.479819860511658591e-01; /* 3FC2F112DF3E5244 */
69 static const double zero
= 0.0;
72 __ieee754_log2 (double x
)
74 double hfsq
, f
, s
, z
, R
, w
, t1
, t2
, dk
;
78 EXTRACT_WORDS64 (hx
, x
);
81 if (hx
< INT64_C(0x0010000000000000))
83 if (__builtin_expect ((hx
& UINT64_C(0x7fffffffffffffff)) == 0, 0))
84 return -two54
/ (x
- x
); /* log(+-0)=-inf */
85 if (__builtin_expect (hx
< 0, 0))
86 return (x
- x
) / (x
- x
); /* log(-#) = NaN */
88 x
*= two54
; /* subnormal number, scale up x */
89 EXTRACT_WORDS64 (hx
, x
);
91 if (__builtin_expect (hx
>= UINT64_C(0x7ff0000000000000), 0))
93 k
+= (hx
>> 52) - 1023;
94 hx
&= UINT64_C(0x000fffffffffffff);
95 i
= (hx
+ UINT64_C(0x95f6400000000)) & UINT64_C(0x10000000000000);
96 /* normalize x or x/2 */
97 INSERT_WORDS64 (x
, hx
| (i
^ UINT64_C(0x3ff0000000000000)));
101 if ((UINT64_C(0x000fffffffffffff) & (2 + hx
)) < 3)
105 R
= f
* f
* (0.5 - 0.33333333333333333 * f
);
106 return dk
- (R
- f
) / ln2
;
110 i
= hx
- UINT64_C(0x6147a00000000);
112 j
= UINT64_C(0x6b85100000000) - hx
;
113 t1
= w
* (Lg2
+ w
* (Lg4
+ w
* Lg6
));
114 t2
= z
* (Lg1
+ w
* (Lg3
+ w
* (Lg5
+ w
* Lg7
)));
120 return dk
- ((hfsq
- (s
* (hfsq
+ R
))) - f
) / ln2
;
124 return dk
- ((s
* (f
- R
)) - f
) / ln2
;
128 strong_alias (__ieee754_log2
, __log2_finite
)