OLD | NEW |
1 // Copyright 2012 the V8 project authors. All rights reserved. | 1 // Copyright 2012 the V8 project authors. All rights reserved. |
2 // Redistribution and use in source and binary forms, with or without | 2 // Redistribution and use in source and binary forms, with or without |
3 // modification, are permitted provided that the following conditions are | 3 // modification, are permitted provided that the following conditions are |
4 // met: | 4 // met: |
5 // | 5 // |
6 // * Redistributions of source code must retain the above copyright | 6 // * Redistributions of source code must retain the above copyright |
7 // notice, this list of conditions and the following disclaimer. | 7 // notice, this list of conditions and the following disclaimer. |
8 // * Redistributions in binary form must reproduce the above | 8 // * Redistributions in binary form must reproduce the above |
9 // copyright notice, this list of conditions and the following | 9 // copyright notice, this list of conditions and the following |
10 // disclaimer in the documentation and/or other materials provided | 10 // disclaimer in the documentation and/or other materials provided |
(...skipping 61 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
72 return %Math_atan2(TO_NUMBER_INLINE(y), TO_NUMBER_INLINE(x)); | 72 return %Math_atan2(TO_NUMBER_INLINE(y), TO_NUMBER_INLINE(x)); |
73 } | 73 } |
74 | 74 |
75 // ECMA 262 - 15.8.2.6 | 75 // ECMA 262 - 15.8.2.6 |
76 function MathCeil(x) { | 76 function MathCeil(x) { |
77 return %Math_ceil(TO_NUMBER_INLINE(x)); | 77 return %Math_ceil(TO_NUMBER_INLINE(x)); |
78 } | 78 } |
79 | 79 |
80 // ECMA 262 - 15.8.2.7 | 80 // ECMA 262 - 15.8.2.7 |
81 function MathCos(x) { | 81 function MathCos(x) { |
82 return MathCosImpl(x); | 82 x = MathAbs(x); // Convert to number and get rid of -0. |
| 83 return TrigonometricInterpolation(x, 1); |
83 } | 84 } |
84 | 85 |
85 // ECMA 262 - 15.8.2.8 | 86 // ECMA 262 - 15.8.2.8 |
86 function MathExp(x) { | 87 function MathExp(x) { |
87 return %Math_exp(TO_NUMBER_INLINE(x)); | 88 return %Math_exp(TO_NUMBER_INLINE(x)); |
88 } | 89 } |
89 | 90 |
90 // ECMA 262 - 15.8.2.9 | 91 // ECMA 262 - 15.8.2.9 |
91 function MathFloor(x) { | 92 function MathFloor(x) { |
92 x = TO_NUMBER_INLINE(x); | 93 x = TO_NUMBER_INLINE(x); |
(...skipping 86 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
179 return (x < 0 ? (x + 0x100000000) : x) * 2.3283064365386962890625e-10; | 180 return (x < 0 ? (x + 0x100000000) : x) * 2.3283064365386962890625e-10; |
180 } | 181 } |
181 | 182 |
182 // ECMA 262 - 15.8.2.15 | 183 // ECMA 262 - 15.8.2.15 |
183 function MathRound(x) { | 184 function MathRound(x) { |
184 return %RoundNumber(TO_NUMBER_INLINE(x)); | 185 return %RoundNumber(TO_NUMBER_INLINE(x)); |
185 } | 186 } |
186 | 187 |
187 // ECMA 262 - 15.8.2.16 | 188 // ECMA 262 - 15.8.2.16 |
188 function MathSin(x) { | 189 function MathSin(x) { |
189 return MathSinImpl(x); | 190 x = x * 1; // Convert to number and deal with -0. |
| 191 if (%_IsMinusZero(x)) return x; |
| 192 return TrigonometricInterpolation(x, 0); |
190 } | 193 } |
191 | 194 |
192 // ECMA 262 - 15.8.2.17 | 195 // ECMA 262 - 15.8.2.17 |
193 function MathSqrt(x) { | 196 function MathSqrt(x) { |
194 return %_MathSqrt(TO_NUMBER_INLINE(x)); | 197 return %_MathSqrt(TO_NUMBER_INLINE(x)); |
195 } | 198 } |
196 | 199 |
197 // ECMA 262 - 15.8.2.18 | 200 // ECMA 262 - 15.8.2.18 |
198 function MathTan(x) { | 201 function MathTan(x) { |
199 return MathSinImpl(x) / MathCosImpl(x); | 202 return MathSin(x) / MathCos(x); |
200 } | 203 } |
201 | 204 |
202 // Non-standard extension. | 205 // Non-standard extension. |
203 function MathImul(x, y) { | 206 function MathImul(x, y) { |
204 return %NumberImul(TO_NUMBER_INLINE(x), TO_NUMBER_INLINE(y)); | 207 return %NumberImul(TO_NUMBER_INLINE(x), TO_NUMBER_INLINE(y)); |
205 } | 208 } |
206 | 209 |
207 | 210 |
208 var MathSinImpl = function(x) { | 211 var kInversePiHalf = 0.636619772367581343; // 2 / pi |
209 InitTrigonometricFunctions(); | 212 var kInversePiHalfS26 = 9.48637384723993156e-9; // 2 / pi / (2^26) |
210 return MathSinImpl(x); | 213 var kS26 = 1 << 26; |
| 214 var kTwoStepThreshold = 1 << 27; |
| 215 // pi / 2 rounded up |
| 216 var kPiHalf = 1.570796326794896780; // 0x192d4454fb21f93f |
| 217 // We use two parts for pi/2 to emulate a higher precision. |
| 218 // pi_half_1 only has 26 significant bits for mantissa. |
| 219 // Note that pi_half > pi_half_1 + pi_half_2 |
| 220 var kPiHalf1 = 1.570796325802803040; // 0x00000054fb21f93f |
| 221 var kPiHalf2 = 9.920935796805404252e-10; // 0x3326a611460b113e |
| 222 |
| 223 var kSamples; // Initialized to a number during genesis. |
| 224 var kIndexConvert; // Initialized to a kSamples / (pi/2) during genesis. |
| 225 var kSinTable; // Initialized to a Float64Array during genesis. |
| 226 var kCosXIntervalTable; // Initialized to a Float64Array during genesis. |
| 227 |
| 228 // This implements sine using the following algorithm. |
| 229 // 1) Multiplication takes care of to-number conversion. |
| 230 // 2) Reduce x to the first quadrant [0, pi/2]. |
| 231 // Conveniently enough, in case of +/-Infinity, we get NaN. |
| 232 // Note that we try to use only 26 instead of 52 significant bits for |
| 233 // mantissa to avoid rounding errors when multiplying. For very large |
| 234 // input we therefore have additional steps. |
| 235 // 3) Replace x by (pi/2-x) if x was in the 2nd or 4th quadrant. |
| 236 // 4) Do a table lookup for the closest samples to the left and right of x. |
| 237 // 5) Find the derivatives at those sampling points by table lookup: |
| 238 // dsin(x)/dx = cos(x) = sin(pi/2-x) for x in [0, pi/2]. |
| 239 // 6) Use cubic spline interpolation to approximate sin(x). |
| 240 // 7) Negate the result if x was in the 3rd or 4th quadrant. |
| 241 // 8) Get rid of -0 by adding 0. |
| 242 function TrigonometricInterpolation(x, phase) { |
| 243 if (x < 0 || x > kPiHalf) { |
| 244 var multiple; |
| 245 while (x < -kTwoStepThreshold || x > kTwoStepThreshold) { |
| 246 // Let's assume this loop does not terminate. |
| 247 // All numbers x in each loop forms a set S. |
| 248 // (1) abs(x) > 2^27 for all x in S. |
| 249 // (2) abs(multiple) != 0 since (2^27 * inverse_pi_half_s26) > 1 |
| 250 // (3) multiple is rounded down in 2^26 steps, so the rounding error is |
| 251 // at most max(ulp, 2^26). |
| 252 // (4) so for x > 2^27, we subtract at most (1+pi/4)x and at least |
| 253 // (1-pi/4)x |
| 254 // (5) The subtraction results in x' so that abs(x') <= abs(x)*pi/4. |
| 255 // Note that this difference cannot be simply rounded off. |
| 256 // Set S cannot exist since (5) violates (1). Loop must terminate. |
| 257 multiple = MathFloor(x * kInversePiHalfS26) * kS26; |
| 258 x = x - multiple * kPiHalf1 - multiple * kPiHalf2; |
| 259 } |
| 260 multiple = MathFloor(x * kInversePiHalf); |
| 261 x = x - multiple * kPiHalf1 - multiple * kPiHalf2; |
| 262 phase += multiple; |
| 263 } |
| 264 var double_index = x * kIndexConvert; |
| 265 if (phase & 1) double_index = kSamples - double_index; |
| 266 var index = double_index | 0; |
| 267 var t1 = double_index - index; |
| 268 var t2 = 1 - t1; |
| 269 var y1 = kSinTable[index]; |
| 270 var y2 = kSinTable[index + 1]; |
| 271 var dy = y2 - y1; |
| 272 return (t2 * y1 + t1 * y2 + |
| 273 t1 * t2 * ((kCosXIntervalTable[index] - dy) * t2 + |
| 274 (dy - kCosXIntervalTable[index + 1]) * t1)) |
| 275 * (1 - (phase & 2)) + 0; |
211 } | 276 } |
212 | 277 |
213 | |
214 var MathCosImpl = function(x) { | |
215 InitTrigonometricFunctions(); | |
216 return MathCosImpl(x); | |
217 } | |
218 | |
219 | |
220 var InitTrigonometricFunctions; | |
221 | |
222 | |
223 // Define constants and interpolation functions. | |
224 // Also define the initialization function that populates the lookup table | |
225 // and then wires up the function definitions. | |
226 function SetupTrigonometricFunctions() { | |
227 var samples = 1800; // Table size. Do not change arbitrarily. | |
228 var inverse_pi_half = 0.636619772367581343; // 2 / pi | |
229 var inverse_pi_half_s_26 = 9.48637384723993156e-9; // 2 / pi / (2^26) | |
230 var s_26 = 1 << 26; | |
231 var two_step_threshold = 1 << 27; | |
232 var index_convert = 1145.915590261646418; // samples / (pi / 2) | |
233 // pi / 2 rounded up | |
234 var pi_half = 1.570796326794896780; // 0x192d4454fb21f93f | |
235 // We use two parts for pi/2 to emulate a higher precision. | |
236 // pi_half_1 only has 26 significant bits for mantissa. | |
237 // Note that pi_half > pi_half_1 + pi_half_2 | |
238 var pi_half_1 = 1.570796325802803040; // 0x00000054fb21f93f | |
239 var pi_half_2 = 9.920935796805404252e-10; // 0x3326a611460b113e | |
240 var table_sin; | |
241 var table_cos_interval; | |
242 | |
243 // This implements sine using the following algorithm. | |
244 // 1) Multiplication takes care of to-number conversion. | |
245 // 2) Reduce x to the first quadrant [0, pi/2]. | |
246 // Conveniently enough, in case of +/-Infinity, we get NaN. | |
247 // Note that we try to use only 26 instead of 52 significant bits for | |
248 // mantissa to avoid rounding errors when multiplying. For very large | |
249 // input we therefore have additional steps. | |
250 // 3) Replace x by (pi/2-x) if x was in the 2nd or 4th quadrant. | |
251 // 4) Do a table lookup for the closest samples to the left and right of x. | |
252 // 5) Find the derivatives at those sampling points by table lookup: | |
253 // dsin(x)/dx = cos(x) = sin(pi/2-x) for x in [0, pi/2]. | |
254 // 6) Use cubic spline interpolation to approximate sin(x). | |
255 // 7) Negate the result if x was in the 3rd or 4th quadrant. | |
256 // 8) Get rid of -0 by adding 0. | |
257 var Interpolation = function(x, phase) { | |
258 if (x < 0 || x > pi_half) { | |
259 var multiple; | |
260 while (x < -two_step_threshold || x > two_step_threshold) { | |
261 // Let's assume this loop does not terminate. | |
262 // All numbers x in each loop forms a set S. | |
263 // (1) abs(x) > 2^27 for all x in S. | |
264 // (2) abs(multiple) != 0 since (2^27 * inverse_pi_half_s26) > 1 | |
265 // (3) multiple is rounded down in 2^26 steps, so the rounding error is | |
266 // at most max(ulp, 2^26). | |
267 // (4) so for x > 2^27, we subtract at most (1+pi/4)x and at least | |
268 // (1-pi/4)x | |
269 // (5) The subtraction results in x' so that abs(x') <= abs(x)*pi/4. | |
270 // Note that this difference cannot be simply rounded off. | |
271 // Set S cannot exist since (5) violates (1). Loop must terminate. | |
272 multiple = MathFloor(x * inverse_pi_half_s_26) * s_26; | |
273 x = x - multiple * pi_half_1 - multiple * pi_half_2; | |
274 } | |
275 multiple = MathFloor(x * inverse_pi_half); | |
276 x = x - multiple * pi_half_1 - multiple * pi_half_2; | |
277 phase += multiple; | |
278 } | |
279 var double_index = x * index_convert; | |
280 if (phase & 1) double_index = samples - double_index; | |
281 var index = double_index | 0; | |
282 var t1 = double_index - index; | |
283 var t2 = 1 - t1; | |
284 var y1 = table_sin[index]; | |
285 var y2 = table_sin[index + 1]; | |
286 var dy = y2 - y1; | |
287 return (t2 * y1 + t1 * y2 + | |
288 t1 * t2 * ((table_cos_interval[index] - dy) * t2 + | |
289 (dy - table_cos_interval[index + 1]) * t1)) | |
290 * (1 - (phase & 2)) + 0; | |
291 } | |
292 | |
293 var MathSinInterpolation = function(x) { | |
294 x = x * 1; // Convert to number and deal with -0. | |
295 if (%_IsMinusZero(x)) return x; | |
296 return Interpolation(x, 0); | |
297 } | |
298 | |
299 // Cosine is sine with a phase offset. | |
300 var MathCosInterpolation = function(x) { | |
301 x = MathAbs(x); // Convert to number and get rid of -0. | |
302 return Interpolation(x, 1); | |
303 }; | |
304 | |
305 %SetInlineBuiltinFlag(Interpolation); | |
306 %SetInlineBuiltinFlag(MathSinInterpolation); | |
307 %SetInlineBuiltinFlag(MathCosInterpolation); | |
308 | |
309 InitTrigonometricFunctions = function() { | |
310 table_sin = new global.Float64Array(samples + 2); | |
311 table_cos_interval = new global.Float64Array(samples + 2); | |
312 %PopulateTrigonometricTable(table_sin, table_cos_interval, samples); | |
313 MathSinImpl = MathSinInterpolation; | |
314 MathCosImpl = MathCosInterpolation; | |
315 } | |
316 } | |
317 | |
318 SetupTrigonometricFunctions(); | |
319 | |
320 | |
321 // ------------------------------------------------------------------- | 278 // ------------------------------------------------------------------- |
322 | 279 |
323 function SetUpMath() { | 280 function SetUpMath() { |
324 %CheckIsBootstrapping(); | 281 %CheckIsBootstrapping(); |
325 | 282 |
326 %SetPrototype($Math, $Object.prototype); | 283 %SetPrototype($Math, $Object.prototype); |
327 %SetProperty(global, "Math", $Math, DONT_ENUM); | 284 %SetProperty(global, "Math", $Math, DONT_ENUM); |
328 %FunctionSetInstanceClassName(MathConstructor, 'Math'); | 285 %FunctionSetInstanceClassName(MathConstructor, 'Math'); |
329 | 286 |
330 // Set up math constants. | 287 // Set up math constants. |
(...skipping 56 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
387 "atan2", MathAtan2, | 344 "atan2", MathAtan2, |
388 "pow", MathPow, | 345 "pow", MathPow, |
389 "max", MathMax, | 346 "max", MathMax, |
390 "min", MathMin, | 347 "min", MathMin, |
391 "imul", MathImul | 348 "imul", MathImul |
392 )); | 349 )); |
393 | 350 |
394 %SetInlineBuiltinFlag(MathSin); | 351 %SetInlineBuiltinFlag(MathSin); |
395 %SetInlineBuiltinFlag(MathCos); | 352 %SetInlineBuiltinFlag(MathCos); |
396 %SetInlineBuiltinFlag(MathTan); | 353 %SetInlineBuiltinFlag(MathTan); |
| 354 %SetInlineBuiltinFlag(TrigonometricInterpolation); |
397 } | 355 } |
398 | 356 |
399 SetUpMath(); | 357 SetUpMath(); |
OLD | NEW |