1 /* Copyright JS Foundation and other contributors, http://js.foundation
2 *
3 * Licensed under the Apache License, Version 2.0 (the "License");
4 * you may not use this file except in compliance with the License.
5 * You may obtain a copy of the License at
6 *
7 * http://www.apache.org/licenses/LICENSE-2.0
8 *
9 * Unless required by applicable law or agreed to in writing, software
10 * distributed under the License is distributed on an "AS IS" BASIS
11 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 * See the License for the specific language governing permissions and
13 * limitations under the License.
14 *
15 * This file is based on work under the following copyright and permission
16 * notice:
17 *
18 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
19 *
20 * Developed at SunSoft, a Sun Microsystems, Inc. business.
21 * Permission to use, copy, modify, and distribute this
22 * software is freely granted, provided that this notice
23 * is preserved.
24 *
25 * @(#)e_log2.c 1.3 95/01/18
26 */
27
28 #include "jerry-libm-internal.h"
29
30 /* log2(x)
31 * Return the base 2 logarithm of x. See e_log.c and k_log.h for most
32 * comments.
33 *
34 * This reduces x to {k, 1+f} exactly as in e_log.c, then calls the kernel,
35 * then does the combining and scaling steps
36 * log2(x) = (f - 0.5*f*f + k_log1p(f)) / ln2 + k
37 * in not-quite-routine extra precision.
38 */
39
40 #define zero 0.0
41 #define two54 1.80143985094819840000e+16 /* 0x43500000, 0x00000000 */
42 #define ivln2hi 1.44269504072144627571e+00 /* 0x3FF71547, 0x65200000 */
43 #define ivln2lo 1.67517131648865118353e-10 /* 0x3DE705FC, 0x2EEFA200 */
44 #define Lg1 6.666666666666735130e-01 /* 0x3FE55555, 0x55555593 */
45 #define Lg2 3.999999999940941908e-01 /* 0x3FD99999, 0x9997FA04 */
46 #define Lg3 2.857142874366239149e-01 /* 0x3FD24924, 0x94229359 */
47 #define Lg4 2.222219843214978396e-01 /* 0x3FCC71C5, 0x1D8E78AF */
48 #define Lg5 1.818357216161805012e-01 /* 0x3FC74664, 0x96CB03DE */
49 #define Lg6 1.531383769920937332e-01 /* 0x3FC39A09, 0xD078C69F */
50 #define Lg7 1.479819860511658591e-01 /* 0x3FC2F112, 0xDF3E5244 */
51
52 double
log2(double x)53 log2 (double x)
54 {
55 double f, hfsq, hi, lo, r, val_hi, val_lo, w, y;
56 int i, k, hx;
57 unsigned int lx;
58 double_accessor temp;
59
60 hx = __HI (x); /* high word of x */
61 lx = __LO (x); /* low word of x */
62
63 k = 0;
64 if (hx < 0x00100000)
65 { /* x < 2**-1022 */
66 if (((hx & 0x7fffffff) | lx) == 0)
67 {
68 return -two54 / zero; /* log(+-0)=-inf */
69 }
70 if (hx < 0)
71 {
72 return (x - x) / zero; /* log(-#) = NaN */
73 }
74 k -= 54;
75 x *= two54; /* subnormal number, scale up x */
76 hx = __HI (x); /* high word of x */
77 }
78 if (hx >= 0x7ff00000)
79 {
80 return x + x;
81 }
82 if (hx == 0x3ff00000 && lx == 0)
83 {
84 return zero; /* log(1) = +0 */
85 }
86 k += (hx >> 20) - 1023;
87 hx &= 0x000fffff;
88 i = (hx + 0x95f64) & 0x100000;
89 temp.dbl = x;
90 temp.as_int.hi = hx | (i ^ 0x3ff00000); /* normalize x or x/2 */
91 k += (i >> 20);
92 y = (double) k;
93 f = temp.dbl - 1.0;
94 hfsq = 0.5 * f * f;
95 double s, z, R, t1, t2;
96
97 s = f / (2.0 + f);
98 z = s * s;
99 w = z * z;
100 t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
101 t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
102 R = t2 + t1;
103 r = s * (hfsq + R);
104 /*
105 * f-hfsq must (for args near 1) be evaluated in extra precision
106 * to avoid a large cancellation when x is near sqrt(2) or 1/sqrt(2).
107 * This is fairly efficient since f-hfsq only depends on f, so can
108 * be evaluated in parallel with R. Not combining hfsq with R also
109 * keeps R small (though not as small as a true `lo' term would be),
110 * so that extra precision is not needed for terms involving R.
111 *
112 * Compiler bugs involving extra precision used to break Dekker's
113 * theorem for spitting f-hfsq as hi+lo, unless double_t was used
114 * or the multi-precision calculations were avoided when double_t
115 * has extra precision. These problems are now automatically
116 * avoided as a side effect of the optimization of combining the
117 * Dekker splitting step with the clear-low-bits step.
118 *
119 * y must (for args near sqrt(2) and 1/sqrt(2)) be added in extra
120 * precision to avoid a very large cancellation when x is very near
121 * these values. Unlike the above cancellations, this problem is
122 * specific to base 2. It is strange that adding +-1 is so much
123 * harder than adding +-ln2 or +-log10_2.
124 *
125 * This uses Dekker's theorem to normalize y+val_hi, so the
126 * compiler bugs are back in some configurations, sigh. And I
127 * don't want to used double_t to avoid them, since that gives a
128 * pessimization and the support for avoiding the pessimization
129 * is not yet available.
130 *
131 * The multi-precision calculations for the multiplications are
132 * routine.
133 */
134 hi = f - hfsq;
135 temp.dbl = hi;
136 temp.as_int.lo = 0;
137
138 lo = (f - hi) - hfsq + r;
139 val_hi = hi * ivln2hi;
140 val_lo = (lo + hi) * ivln2lo + lo * ivln2hi;
141
142 /* spadd(val_hi, val_lo, y), except for not using double_t: */
143 w = y + val_hi;
144 val_lo += (y - w) + val_hi;
145 val_hi = w;
146
147 return val_lo + val_hi;
148 } /* log2 */
149
150 #undef zero
151 #undef two54
152 #undef ivln2hi
153 #undef ivln2lo
154 #undef Lg1
155 #undef Lg2
156 #undef Lg3
157 #undef Lg4
158 #undef Lg5
159 #undef Lg6
160 #undef Lg7
161