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 199 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
210 } | 210 } |
211 | 211 |
212 | 212 |
213 var InitTrigonometricFunctions; | 213 var InitTrigonometricFunctions; |
214 | 214 |
215 | 215 |
216 // Define constants and interpolation functions. | 216 // Define constants and interpolation functions. |
217 // Also define the initialization function that populates the lookup table | 217 // Also define the initialization function that populates the lookup table |
218 // and then wires up the function definitions. | 218 // and then wires up the function definitions. |
219 function SetupTrigonometricFunctions() { | 219 function SetupTrigonometricFunctions() { |
220 // TODO(yangguo): The following table size has been chosen to satisfy | 220 var samples = 1800; // Table size. Do not change arbitrarily. |
221 // Sunspider's brittle result verification. Reconsider relevance. | 221 var inverse_pi_half = 0.636619772367581343; // 2 / pi |
222 var samples = 4489; | 222 var inverse_pi_half_s_26 = 9.48637384723993156e-9; // 2 / pi / (2^26) |
223 var pi = 3.1415926535897932; | 223 var s_26 = 1 << 26; |
224 var pi_half = pi / 2; | 224 var two_step_threshold = 1 << 27; |
225 var inverse_pi_half = 2 / pi; | 225 var index_convert = 1145.915590261646418; // samples / (pi / 2) |
226 var two_pi = 2 * pi; | 226 // pi / 2 rounded up |
227 var four_pi = 4 * pi; | 227 var pi_half = 1.570796326794896780; // 0x192d4454fb21f93f |
228 var interval = pi_half / samples; | 228 // We use two parts for pi/2 to emulate a higher precision. |
229 var inverse_interval = samples / pi_half; | 229 // pi_half_1 only has 26 significant bits for mantissa. |
| 230 // Note that pi_half > pi_half_1 + pi_half_2 |
| 231 var pi_half_1 = 1.570796325802803040; // 0x00000054fb21f93f |
| 232 var pi_half_2 = 9.920935796805404252e-10; // 0x3326a611460b113e |
230 var table_sin; | 233 var table_sin; |
231 var table_cos_interval; | 234 var table_cos_interval; |
232 | 235 |
233 // This implements sine using the following algorithm. | 236 // This implements sine using the following algorithm. |
234 // 1) Multiplication takes care of to-number conversion. | 237 // 1) Multiplication takes care of to-number conversion. |
235 // 2) Reduce x to the first quadrant [0, pi/2]. | 238 // 2) Reduce x to the first quadrant [0, pi/2]. |
236 // Conveniently enough, in case of +/-Infinity, we get NaN. | 239 // Conveniently enough, in case of +/-Infinity, we get NaN. |
| 240 // Note that we try to use only 26 instead of 52 significant bits for |
| 241 // mantissa to avoid rounding errors when multiplying. For very large |
| 242 // input we therefore have additional steps. |
237 // 3) Replace x by (pi/2-x) if x was in the 2nd or 4th quadrant. | 243 // 3) Replace x by (pi/2-x) if x was in the 2nd or 4th quadrant. |
238 // 4) Do a table lookup for the closest samples to the left and right of x. | 244 // 4) Do a table lookup for the closest samples to the left and right of x. |
239 // 5) Find the derivatives at those sampling points by table lookup: | 245 // 5) Find the derivatives at those sampling points by table lookup: |
240 // dsin(x)/dx = cos(x) = sin(pi/2-x) for x in [0, pi/2]. | 246 // dsin(x)/dx = cos(x) = sin(pi/2-x) for x in [0, pi/2]. |
241 // 6) Use cubic spline interpolation to approximate sin(x). | 247 // 6) Use cubic spline interpolation to approximate sin(x). |
242 // 7) Negate the result if x was in the 3rd or 4th quadrant. | 248 // 7) Negate the result if x was in the 3rd or 4th quadrant. |
243 // 8) Get rid of -0 by adding 0. | 249 // 8) Get rid of -0 by adding 0. |
244 var Interpolation = function(x) { | 250 var Interpolation = function(x, phase) { |
245 var double_index = x * inverse_interval; | 251 if (x < 0 || x > pi_half) { |
| 252 var multiple; |
| 253 while (x < -two_step_threshold || x > two_step_threshold) { |
| 254 // Let's assume the loop does not terminate. |
| 255 // All numbers x in each loop forms a set S. |
| 256 // abs(x) > 2^27 for all x in S. |
| 257 // abs(multiple) != 0 since (2^27 * inverse_pi_half_s26) > 1 |
| 258 // multiple is rounded down in 2^26 steps, so the rounding error is at |
| 259 // most max(ulp, 2^26), so for x > 2^27, we subtract at most 3/2 x |
| 260 // and at least 1/2 x. We end up with x' so that abs(x') <= abs(x)/2 |
| 261 // Note that since the difference is at least x/2, it cannot be simply |
| 262 // rounded off. |
| 263 // Such a set S cannot exist. |
| 264 multiple = MathFloor(x * inverse_pi_half_s_26) * s_26; |
| 265 x = x - multiple * pi_half_1 - multiple * pi_half_2; |
| 266 } |
| 267 multiple = MathFloor(x * inverse_pi_half); |
| 268 x = x - multiple * pi_half_1 - multiple * pi_half_2; |
| 269 phase += multiple; |
| 270 } |
| 271 var double_index = x * index_convert; |
| 272 if (phase & 1) double_index = samples - double_index; |
246 var index = double_index | 0; | 273 var index = double_index | 0; |
247 var t1 = double_index - index; | 274 var t1 = double_index - index; |
248 var t2 = 1 - t1; | 275 var t2 = 1 - t1; |
249 var y1 = table_sin[index]; | 276 var y1 = table_sin[index]; |
250 var y2 = table_sin[index + 1]; | 277 var y2 = table_sin[index + 1]; |
251 var dy = y2 - y1; | 278 var dy = y2 - y1; |
252 return (t2 * y1 + t1 * y2 + | 279 return (t2 * y1 + t1 * y2 + |
253 t1 * t2 * ((table_cos_interval[index] - dy) * t2 + | 280 t1 * t2 * ((table_cos_interval[index] - dy) * t2 + |
254 (dy - table_cos_interval[index + 1]) * t1)); | 281 (dy - table_cos_interval[index + 1]) * t1)) |
| 282 * (1 - (phase & 2)) + 0; |
255 } | 283 } |
256 | 284 |
257 var MathSinInterpolation = function(x) { | 285 var MathSinInterpolation = function(x) { |
258 // This is to make Sunspider's result verification happy. | 286 x = x * 1; // Convert to number and deal with -0. |
259 if (x > four_pi) x -= four_pi; | 287 if (%_IsMinusZero(x)) return x; |
260 var multiple = MathFloor(x * inverse_pi_half); | 288 return Interpolation(x, 0); |
261 if (%_IsMinusZero(multiple)) return multiple; | |
262 x = (multiple & 1) * pi_half + | |
263 (1 - ((multiple & 1) << 1)) * (x - multiple * pi_half); | |
264 return Interpolation(x) * (1 - (multiple & 2)) + 0; | |
265 } | 289 } |
266 | 290 |
267 // Cosine is sine with a phase offset of pi/2. | 291 // Cosine is sine with a phase offset. |
268 var MathCosInterpolation = function(x) { | 292 var MathCosInterpolation = function(x) { |
269 var multiple = MathFloor(x * inverse_pi_half); | 293 x = MathAbs(x); // Convert to number and get rid of -0. |
270 var phase = multiple + 1; | 294 return Interpolation(x, 1); |
271 x = (phase & 1) * pi_half + | |
272 (1 - ((phase & 1) << 1)) * (x - multiple * pi_half); | |
273 return Interpolation(x) * (1 - (phase & 2)) + 0; | |
274 }; | 295 }; |
275 | 296 |
276 %SetInlineBuiltinFlag(Interpolation); | 297 %SetInlineBuiltinFlag(Interpolation); |
277 %SetInlineBuiltinFlag(MathSinInterpolation); | 298 %SetInlineBuiltinFlag(MathSinInterpolation); |
278 %SetInlineBuiltinFlag(MathCosInterpolation); | 299 %SetInlineBuiltinFlag(MathCosInterpolation); |
279 | 300 |
280 InitTrigonometricFunctions = function() { | 301 InitTrigonometricFunctions = function() { |
281 table_sin = new global.Float64Array(samples + 2); | 302 table_sin = new global.Float64Array(samples + 2); |
282 table_cos_interval = new global.Float64Array(samples + 2); | 303 table_cos_interval = new global.Float64Array(samples + 2); |
283 %PopulateTrigonometricTable(table_sin, table_cos_interval, samples); | 304 %PopulateTrigonometricTable(table_sin, table_cos_interval, samples); |
(...skipping 77 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
361 "min", MathMin, | 382 "min", MathMin, |
362 "imul", MathImul | 383 "imul", MathImul |
363 )); | 384 )); |
364 | 385 |
365 %SetInlineBuiltinFlag(MathSin); | 386 %SetInlineBuiltinFlag(MathSin); |
366 %SetInlineBuiltinFlag(MathCos); | 387 %SetInlineBuiltinFlag(MathCos); |
367 %SetInlineBuiltinFlag(MathTan); | 388 %SetInlineBuiltinFlag(MathTan); |
368 } | 389 } |
369 | 390 |
370 SetUpMath(); | 391 SetUpMath(); |
OLD | NEW |