OLD | NEW |
| (Empty) |
1 // The following is adapted from fdlibm (http://www.netlib.org/fdlibm). | |
2 // | |
3 // ==================================================== | |
4 // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | |
5 // | |
6 // Developed at SunSoft, a Sun Microsystems, Inc. business. | |
7 // Permission to use, copy, modify, and distribute this | |
8 // software is freely granted, provided that this notice | |
9 // is preserved. | |
10 // ==================================================== | |
11 // | |
12 // The original source code covered by the above license above has been | |
13 // modified significantly by Google Inc. | |
14 // Copyright 2014 the V8 project authors. All rights reserved. | |
15 | |
16 #include "src/v8.h" | |
17 | |
18 #include "src/double.h" | |
19 #include "third_party/fdlibm/fdlibm.h" | |
20 | |
21 | |
22 namespace v8 { | |
23 namespace fdlibm { | |
24 | |
25 #ifdef _MSC_VER | |
26 inline double scalbn(double x, int y) { return _scalb(x, y); } | |
27 #endif // _MSC_VER | |
28 | |
29 const double MathConstants::constants[] = { | |
30 6.36619772367581382433e-01, // invpio2 0 | |
31 1.57079632673412561417e+00, // pio2_1 1 | |
32 6.07710050650619224932e-11, // pio2_1t 2 | |
33 6.07710050630396597660e-11, // pio2_2 3 | |
34 2.02226624879595063154e-21, // pio2_2t 4 | |
35 2.02226624871116645580e-21, // pio2_3 5 | |
36 8.47842766036889956997e-32, // pio2_3t 6 | |
37 -1.66666666666666324348e-01, // S1 7 coefficients for sin | |
38 8.33333333332248946124e-03, // 8 | |
39 -1.98412698298579493134e-04, // 9 | |
40 2.75573137070700676789e-06, // 10 | |
41 -2.50507602534068634195e-08, // 11 | |
42 1.58969099521155010221e-10, // S6 12 | |
43 4.16666666666666019037e-02, // C1 13 coefficients for cos | |
44 -1.38888888888741095749e-03, // 14 | |
45 2.48015872894767294178e-05, // 15 | |
46 -2.75573143513906633035e-07, // 16 | |
47 2.08757232129817482790e-09, // 17 | |
48 -1.13596475577881948265e-11, // C6 18 | |
49 3.33333333333334091986e-01, // T0 19 coefficients for tan | |
50 1.33333333333201242699e-01, // 20 | |
51 5.39682539762260521377e-02, // 21 | |
52 2.18694882948595424599e-02, // 22 | |
53 8.86323982359930005737e-03, // 23 | |
54 3.59207910759131235356e-03, // 24 | |
55 1.45620945432529025516e-03, // 25 | |
56 5.88041240820264096874e-04, // 26 | |
57 2.46463134818469906812e-04, // 27 | |
58 7.81794442939557092300e-05, // 28 | |
59 7.14072491382608190305e-05, // 29 | |
60 -1.85586374855275456654e-05, // 30 | |
61 2.59073051863633712884e-05, // T12 31 | |
62 7.85398163397448278999e-01, // pio4 32 | |
63 3.06161699786838301793e-17, // pio4lo 33 | |
64 6.93147180369123816490e-01, // ln2_hi 34 | |
65 1.90821492927058770002e-10, // ln2_lo 35 | |
66 1.80143985094819840000e+16, // 2^54 36 | |
67 6.666666666666666666e-01, // 2/3 37 | |
68 6.666666666666735130e-01, // LP1 38 coefficients for log1p | |
69 3.999999999940941908e-01, // 39 | |
70 2.857142874366239149e-01, // 40 | |
71 2.222219843214978396e-01, // 41 | |
72 1.818357216161805012e-01, // 42 | |
73 1.531383769920937332e-01, // 43 | |
74 1.479819860511658591e-01, // LP7 44 | |
75 7.09782712893383973096e+02, // 45 overflow threshold for expm1 | |
76 1.44269504088896338700e+00, // 1/ln2 46 | |
77 -3.33333333333331316428e-02, // Q1 47 coefficients for expm1 | |
78 1.58730158725481460165e-03, // 48 | |
79 -7.93650757867487942473e-05, // 49 | |
80 4.00821782732936239552e-06, // 50 | |
81 -2.01099218183624371326e-07, // Q5 51 | |
82 710.4758600739439 // 52 overflow threshold sinh, cosh | |
83 }; | |
84 | |
85 | |
86 // Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi | |
87 static const int two_over_pi[] = { | |
88 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C, | |
89 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 0x424DD2, 0xE00649, | |
90 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 0xA73EE8, 0x8235F5, 0x2EBB44, | |
91 0x84E99C, 0x7026B4, 0x5F7E41, 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, | |
92 0xBDF928, 0x3B1FF8, 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, | |
93 0x367ECF, 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, | |
94 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 0x560330, | |
95 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 0x91615E, 0xE61B08, | |
96 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 0x4D7327, 0x310606, 0x1556CA, | |
97 0x73A8C9, 0x60E27B, 0xC08C6B}; | |
98 | |
99 static const double zero = 0.0; | |
100 static const double two24 = 1.6777216e+07; | |
101 static const double one = 1.0; | |
102 static const double twon24 = 5.9604644775390625e-08; | |
103 | |
104 static const double PIo2[] = { | |
105 1.57079625129699707031e+00, // 0x3FF921FB, 0x40000000 | |
106 7.54978941586159635335e-08, // 0x3E74442D, 0x00000000 | |
107 5.39030252995776476554e-15, // 0x3CF84698, 0x80000000 | |
108 3.28200341580791294123e-22, // 0x3B78CC51, 0x60000000 | |
109 1.27065575308067607349e-29, // 0x39F01B83, 0x80000000 | |
110 1.22933308981111328932e-36, // 0x387A2520, 0x40000000 | |
111 2.73370053816464559624e-44, // 0x36E38222, 0x80000000 | |
112 2.16741683877804819444e-51 // 0x3569F31D, 0x00000000 | |
113 }; | |
114 | |
115 | |
116 int __kernel_rem_pio2(double* x, double* y, int e0, int nx) { | |
117 static const int32_t jk = 3; | |
118 double fw; | |
119 int32_t jx = nx - 1; | |
120 int32_t jv = (e0 - 3) / 24; | |
121 if (jv < 0) jv = 0; | |
122 int32_t q0 = e0 - 24 * (jv + 1); | |
123 int32_t m = jx + jk; | |
124 | |
125 double f[10]; | |
126 for (int i = 0, j = jv - jx; i <= m; i++, j++) { | |
127 f[i] = (j < 0) ? zero : static_cast<double>(two_over_pi[j]); | |
128 } | |
129 | |
130 double q[10]; | |
131 for (int i = 0; i <= jk; i++) { | |
132 fw = 0.0; | |
133 for (int j = 0; j <= jx; j++) fw += x[j] * f[jx + i - j]; | |
134 q[i] = fw; | |
135 } | |
136 | |
137 int32_t jz = jk; | |
138 | |
139 recompute: | |
140 | |
141 int32_t iq[10]; | |
142 double z = q[jz]; | |
143 for (int i = 0, j = jz; j > 0; i++, j--) { | |
144 fw = static_cast<double>(static_cast<int32_t>(twon24 * z)); | |
145 iq[i] = static_cast<int32_t>(z - two24 * fw); | |
146 z = q[j - 1] + fw; | |
147 } | |
148 | |
149 z = scalbn(z, q0); | |
150 z -= 8.0 * std::floor(z * 0.125); | |
151 int32_t n = static_cast<int32_t>(z); | |
152 z -= static_cast<double>(n); | |
153 int32_t ih = 0; | |
154 if (q0 > 0) { | |
155 int32_t i = (iq[jz - 1] >> (24 - q0)); | |
156 n += i; | |
157 iq[jz - 1] -= i << (24 - q0); | |
158 ih = iq[jz - 1] >> (23 - q0); | |
159 } else if (q0 == 0) { | |
160 ih = iq[jz - 1] >> 23; | |
161 } else if (z >= 0.5) { | |
162 ih = 2; | |
163 } | |
164 | |
165 if (ih > 0) { | |
166 n += 1; | |
167 int32_t carry = 0; | |
168 for (int i = 0; i < jz; i++) { | |
169 int32_t j = iq[i]; | |
170 if (carry == 0) { | |
171 if (j != 0) { | |
172 carry = 1; | |
173 iq[i] = 0x1000000 - j; | |
174 } | |
175 } else { | |
176 iq[i] = 0xffffff - j; | |
177 } | |
178 } | |
179 if (q0 == 1) { | |
180 iq[jz - 1] &= 0x7fffff; | |
181 } else if (q0 == 2) { | |
182 iq[jz - 1] &= 0x3fffff; | |
183 } | |
184 if (ih == 2) { | |
185 z = one - z; | |
186 if (carry != 0) z -= scalbn(one, q0); | |
187 } | |
188 } | |
189 | |
190 if (z == zero) { | |
191 int32_t j = 0; | |
192 for (int i = jz - 1; i >= jk; i--) j |= iq[i]; | |
193 if (j == 0) { | |
194 int32_t k = 1; | |
195 while (iq[jk - k] == 0) k++; | |
196 for (int i = jz + 1; i <= jz + k; i++) { | |
197 f[jx + i] = static_cast<double>(two_over_pi[jv + i]); | |
198 for (j = 0, fw = 0.0; j <= jx; j++) fw += x[j] * f[jx + i - j]; | |
199 q[i] = fw; | |
200 } | |
201 jz += k; | |
202 goto recompute; | |
203 } | |
204 } | |
205 | |
206 if (z == 0.0) { | |
207 jz -= 1; | |
208 q0 -= 24; | |
209 while (iq[jz] == 0) { | |
210 jz--; | |
211 q0 -= 24; | |
212 } | |
213 } else { | |
214 z = scalbn(z, -q0); | |
215 if (z >= two24) { | |
216 fw = static_cast<double>(static_cast<int32_t>(twon24 * z)); | |
217 iq[jz] = static_cast<int32_t>(z - two24 * fw); | |
218 jz += 1; | |
219 q0 += 24; | |
220 iq[jz] = static_cast<int32_t>(fw); | |
221 } else { | |
222 iq[jz] = static_cast<int32_t>(z); | |
223 } | |
224 } | |
225 | |
226 fw = scalbn(one, q0); | |
227 for (int i = jz; i >= 0; i--) { | |
228 q[i] = fw * static_cast<double>(iq[i]); | |
229 fw *= twon24; | |
230 } | |
231 | |
232 double fq[10]; | |
233 for (int i = jz; i >= 0; i--) { | |
234 fw = 0.0; | |
235 for (int k = 0; k <= jk && k <= jz - i; k++) fw += PIo2[k] * q[i + k]; | |
236 fq[jz - i] = fw; | |
237 } | |
238 | |
239 fw = 0.0; | |
240 for (int i = jz; i >= 0; i--) fw += fq[i]; | |
241 y[0] = (ih == 0) ? fw : -fw; | |
242 fw = fq[0] - fw; | |
243 for (int i = 1; i <= jz; i++) fw += fq[i]; | |
244 y[1] = (ih == 0) ? fw : -fw; | |
245 return n & 7; | |
246 } | |
247 | |
248 | |
249 int rempio2(double x, double* y) { | |
250 int32_t hx = static_cast<int32_t>(internal::double_to_uint64(x) >> 32); | |
251 int32_t ix = hx & 0x7fffffff; | |
252 | |
253 if (ix >= 0x7ff00000) { | |
254 *y = base::OS::nan_value(); | |
255 return 0; | |
256 } | |
257 | |
258 int32_t e0 = (ix >> 20) - 1046; | |
259 uint64_t zi = internal::double_to_uint64(x) & 0xFFFFFFFFu; | |
260 zi |= static_cast<uint64_t>(ix - (e0 << 20)) << 32; | |
261 double z = internal::uint64_to_double(zi); | |
262 | |
263 double tx[3]; | |
264 for (int i = 0; i < 2; i++) { | |
265 tx[i] = static_cast<double>(static_cast<int32_t>(z)); | |
266 z = (z - tx[i]) * two24; | |
267 } | |
268 tx[2] = z; | |
269 | |
270 int nx = 3; | |
271 while (tx[nx - 1] == zero) nx--; | |
272 int n = __kernel_rem_pio2(tx, y, e0, nx); | |
273 if (hx < 0) { | |
274 y[0] = -y[0]; | |
275 y[1] = -y[1]; | |
276 return -n; | |
277 } | |
278 return n; | |
279 } | |
280 } | |
281 } // namespace v8::internal | |
OLD | NEW |