@@ -79,7 +79,8 @@ function MathCeil(x) {
7979
8080// ECMA 262 - 15.8.2.7
8181function MathCos ( x ) {
82- return MathCosImpl ( x ) ;
82+ x = MathAbs ( x ) ; // Convert to number and get rid of -0.
83+ return TrigonometricInterpolation ( x , 1 ) ;
8384}
8485
8586// ECMA 262 - 15.8.2.8
@@ -179,7 +180,9 @@ function MathRound(x) {
179180
180181// ECMA 262 - 15.8.2.16
181182function MathSin ( x ) {
182- return MathSinImpl ( x ) ;
183+ x = x * 1 ; // Convert to number and deal with -0.
184+ if ( % _IsMinusZero ( x ) ) return x ;
185+ return TrigonometricInterpolation ( x , 0 ) ;
183186}
184187
185188// ECMA 262 - 15.8.2.17
@@ -189,7 +192,7 @@ function MathSqrt(x) {
189192
190193// ECMA 262 - 15.8.2.18
191194function MathTan ( x ) {
192- return MathSinImpl ( x ) / MathCosImpl ( x ) ;
195+ return MathSin ( x ) / MathCos ( x ) ;
193196}
194197
195198// Non-standard extension.
@@ -198,119 +201,73 @@ function MathImul(x, y) {
198201}
199202
200203
201- var MathSinImpl = function ( x ) {
202- InitTrigonometricFunctions ( ) ;
203- return MathSinImpl ( x ) ;
204- }
205-
206-
207- var MathCosImpl = function ( x ) {
208- InitTrigonometricFunctions ( ) ;
209- return MathCosImpl ( x ) ;
210- }
211-
212-
213- var InitTrigonometricFunctions ;
214-
215-
216- // Define constants and interpolation functions.
217- // Also define the initialization function that populates the lookup table
218- // and then wires up the function definitions.
219- function SetupTrigonometricFunctions ( ) {
220- var samples = 1800 ; // Table size. Do not change arbitrarily.
221- var inverse_pi_half = 0.636619772367581343 ; // 2 / pi
222- var inverse_pi_half_s_26 = 9.48637384723993156e-9 ; // 2 / pi / (2^26)
223- var s_26 = 1 << 26 ;
224- var two_step_threshold = 1 << 27 ;
225- var index_convert = 1145.915590261646418 ; // samples / (pi / 2)
226- // pi / 2 rounded up
227- var pi_half = 1.570796326794896780 ; // 0x192d4454fb21f93f
228- // We use two parts for pi/2 to emulate a higher precision.
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
233- var table_sin ;
234- var table_cos_interval ;
235-
236- // This implements sine using the following algorithm.
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 ;
204+ var kInversePiHalf = 0.636619772367581343 ; // 2 / pi
205+ var kInversePiHalfS26 = 9.48637384723993156e-9 ; // 2 / pi / (2^26)
206+ var kS26 = 1 << 26 ;
207+ var kTwoStepThreshold = 1 << 27 ;
208+ // pi / 2 rounded up
209+ var kPiHalf = 1.570796326794896780 ; // 0x192d4454fb21f93f
210+ // We use two parts for pi/2 to emulate a higher precision.
211+ // pi_half_1 only has 26 significant bits for mantissa.
212+ // Note that pi_half > pi_half_1 + pi_half_2
213+ var kPiHalf1 = 1.570796325802803040 ; // 0x00000054fb21f93f
214+ var kPiHalf2 = 9.920935796805404252e-10 ; // 0x3326a611460b113e
215+
216+ var kSamples ; // Initialized to a number during genesis.
217+ var kIndexConvert ; // Initialized to kSamples / (pi/2) during genesis.
218+ var kSinTable ; // Initialized to a Float64Array during genesis.
219+ var kCosXIntervalTable ; // Initialized to a Float64Array during genesis.
220+
221+ // This implements sine using the following algorithm.
222+ // 1) Multiplication takes care of to-number conversion.
223+ // 2) Reduce x to the first quadrant [0, pi/2].
224+ // Conveniently enough, in case of +/-Infinity, we get NaN.
225+ // Note that we try to use only 26 instead of 52 significant bits for
226+ // mantissa to avoid rounding errors when multiplying. For very large
227+ // input we therefore have additional steps.
228+ // 3) Replace x by (pi/2-x) if x was in the 2nd or 4th quadrant.
229+ // 4) Do a table lookup for the closest samples to the left and right of x.
230+ // 5) Find the derivatives at those sampling points by table lookup:
231+ // dsin(x)/dx = cos(x) = sin(pi/2-x) for x in [0, pi/2].
232+ // 6) Use cubic spline interpolation to approximate sin(x).
233+ // 7) Negate the result if x was in the 3rd or 4th quadrant.
234+ // 8) Get rid of -0 by adding 0.
235+ function TrigonometricInterpolation ( x , phase ) {
236+ if ( x < 0 || x > kPiHalf ) {
237+ var multiple ;
238+ while ( x < - kTwoStepThreshold || x > kTwoStepThreshold ) {
239+ // Let's assume this loop does not terminate.
240+ // All numbers x in each loop forms a set S.
241+ // (1) abs(x) > 2^27 for all x in S.
242+ // (2) abs(multiple) != 0 since (2^27 * inverse_pi_half_s26) > 1
243+ // (3) multiple is rounded down in 2^26 steps, so the rounding error is
244+ // at most max(ulp, 2^26).
245+ // (4) so for x > 2^27, we subtract at most (1+pi/4)x and at least
246+ // (1-pi/4)x
247+ // (5) The subtraction results in x' so that abs(x') <= abs(x)*pi/4.
248+ // Note that this difference cannot be simply rounded off.
249+ // Set S cannot exist since (5) violates (1). Loop must terminate.
250+ multiple = MathFloor ( x * kInversePiHalfS26 ) * kS26 ;
251+ x = x - multiple * kPiHalf1 - multiple * kPiHalf2 ;
271252 }
272- var double_index = x * index_convert ;
273- if ( phase & 1 ) double_index = samples - double_index ;
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 ;
284- }
285-
286- var MathSinInterpolation = function ( x ) {
287- x = x * 1 ; // Convert to number and deal with -0.
288- if ( % _IsMinusZero ( x ) ) return x ;
289- return Interpolation ( x , 0 ) ;
290- }
291-
292- // Cosine is sine with a phase offset.
293- var MathCosInterpolation = function ( x ) {
294- x = MathAbs ( x ) ; // Convert to number and get rid of -0.
295- return Interpolation ( x , 1 ) ;
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 ;
253+ multiple = MathFloor ( x * kInversePiHalf ) ;
254+ x = x - multiple * kPiHalf1 - multiple * kPiHalf2 ;
255+ phase += multiple ;
308256 }
257+ var double_index = x * kIndexConvert ;
258+ if ( phase & 1 ) double_index = kSamples - double_index ;
259+ var index = double_index | 0 ;
260+ var t1 = double_index - index ;
261+ var t2 = 1 - t1 ;
262+ var y1 = kSinTable [ index ] ;
263+ var y2 = kSinTable [ index + 1 ] ;
264+ var dy = y2 - y1 ;
265+ return ( t2 * y1 + t1 * y2 +
266+ t1 * t2 * ( ( kCosXIntervalTable [ index ] - dy ) * t2 +
267+ ( dy - kCosXIntervalTable [ index + 1 ] ) * t1 ) )
268+ * ( 1 - ( phase & 2 ) ) + 0 ;
309269}
310270
311- SetupTrigonometricFunctions ( ) ;
312-
313-
314271// -------------------------------------------------------------------
315272
316273function SetUpMath ( ) {
@@ -387,6 +344,7 @@ function SetUpMath() {
387344 % SetInlineBuiltinFlag ( MathSin ) ;
388345 % SetInlineBuiltinFlag ( MathCos ) ;
389346 % SetInlineBuiltinFlag ( MathTan ) ;
347+ % SetInlineBuiltinFlag ( TrigonometricInterpolation ) ;
390348}
391349
392350SetUpMath ( ) ;
0 commit comments