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