Implement trigonometric functions using a fdlibm port.

src/math.js

Index: src/math.js
diff --git a/src/math.js b/src/math.js
index 9dc4b37d0ce2115ed9e9b5078f204fd060f926b9..370bff6870fb04beaf25c321a18859b1ca1c5371 100644
--- a/src/math.js
+++ b/src/math.js
@@ -56,12 +56,6 @@ function MathCeil(x) {
   return -MathFloor(-x);
 }
 
-// ECMA 262 - 15.8.2.7
-function MathCos(x) {
-  x = MathAbs(x);  // Convert to number and get rid of -0.
-  return TrigonometricInterpolation(x, 1);
-}
-
 // ECMA 262 - 15.8.2.8
 function MathExp(x) {
   return %MathExpRT(TO_NUMBER_INLINE(x));
@@ -164,94 +158,252 @@ function MathRound(x) {
   return %RoundNumber(TO_NUMBER_INLINE(x));
 }
 
-// ECMA 262 - 15.8.2.16
-function MathSin(x) {
-  x = x * 1;  // Convert to number and deal with -0.
-  if (%_IsMinusZero(x)) return x;
-  return TrigonometricInterpolation(x, 0);
-}
-
 // ECMA 262 - 15.8.2.17
 function MathSqrt(x) {
   return %_MathSqrtRT(TO_NUMBER_INLINE(x));
 }
 
-// ECMA 262 - 15.8.2.18
-function MathTan(x) {
-  return MathSin(x) / MathCos(x);
-}
-
 // Non-standard extension.
 function MathImul(x, y) {
   return %NumberImul(TO_NUMBER_INLINE(x), TO_NUMBER_INLINE(y));
 }
 
+// -------------------------------------------------------------------
+
+// A straightforward translation of fdlibm routines for sin, cos, and
+// tan, by Raymond Toy (rtoy@google.com).
+
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunSoft, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+
+// Initialized to a Float64Array during genesis and is not writable.
+var kTrig;
+
+macro REMPIO2(X)

[Raymond Toy, 2014/07/24 18:08:44]

Does this have to be a macro? Dropping this big chunk of code all over would
seem to hurt the instruction cache on mobile devices.  Is there a huge loss in
performance if this were a regular function?

[Yang, 2014/07/28 11:28:59]

Yes. Calling a function with double arguments in V8 requires boxing it into a
HeapNumber. We also have three return values instead of one. So inlining it as a
macro is indeed (I tested) faster.

[Raymond Toy, 2014/07/30 19:55:04]

Ok, that's not too surprising if your tests stress trig functions.

I was curious on your thoughts about blowing out the instruction cache on low
end devices, but I don't have any good tests for this.

+  var n, y0, y1;
+  var hx = %_DoubleHi(X);
+  var ix = hx & 0x7fffffff;
+
+  if (ix < 0x4002d97c) {
+    // |X| ~< 3*pi/4, special case with n = +/- 1
+    if (hx > 0) {
+      var z = X - kTrig[1];

[Raymond Toy, 2014/07/24 16:51:11]

Readability would be improved if the kTrig values used here were actual
constants. Now it's pretty hard to know that kTrig[2] was pio2_1t. (That's not
really clear either, but at least it hints that the value is probably related to
pi/2, which it is.)

[Yang, 2014/07/28 11:28:59]

That is true. However, given that we currently don't have constant pools for
ia32 and x64, this speeds up loading double constants a lot. I would like to
keep it this way. I'll add a comment to explain where to look it up.

[Raymond Toy, 2014/07/30 19:55:04]

Could you give symbolic names for the indices? So instead of kTrig[2], use
kTrig[pio2_1t]? Doesn't solve all readability issues, but does make it a bit
more obvious.

[Yang, 2014/08/01 07:29:55]

Done with macros

+      if (ix != 0x3ff921fb) {
+        // 33+53 bit pi is good enough
+        y0 = z - kTrig[2];
+        y1 = (z - y0) - kTrig[2];
+      } else {
+        // near pi/2, use 33+33+53 bit pi
+        z -= kTrig[3];
+        y0 = z - kTrig[4];
+        y1 = (z - y0) - kTrig[4];
+      }
+      n = 1;
+    } else {
+      // Negative X
+      var z = X + kTrig[1];
+      if (ix != 0x3ff921fb) {
+        // 33+53 bit pi is good enough
+        y0 = z + kTrig[2];
+        y1 = (z - y0) + kTrig[2];
+      } else {
+        // near pi/2, use 33+33+53 bit pi
+        z += kTrig[3];
+        y0 = z + kTrig[4];
+        y1 = (z - y0) + kTrig[4];
+      }
+      n = -1;
+    }
+  } else if (ix <= 0x413921fb) {
+    // |X| ~<= 2^19*(pi/2), medium size
+    var t = MathAbs(X);
+    n = (t * kTrig[0] + 0.5) | 0;
+    var r = t - n * kTrig[1];
+    var w = n * kTrig[2];
+    // First round good to 85 bit
+    y0 = r - w;
+    if (ix - (%_DoubleHi(y0) & 0x7ff00000) > 0x1000000) {

[Raymond Toy, 2014/07/30 19:55:04]

It would certainly be good to document what this is doing. I know the original
didn't, but that doesn't mean we shouldn't. I can help with this.

[Yang, 2014/08/01 07:29:55]

I don't think we should be spending time documenting third_party code, let alone
one with one-letter variable names.

Feel free to contribute, but I'd like to not to have this CL blocked by this.

[Raymond Toy, 2014/08/01 14:06:00]

I would agree if this were in fact the original code, but it's now a translation
of the C code to Javascript.
Fair enough.

FWIW, I think the %_DoubleHi & 0x7ff00000 is extracting the (biased) exponent
from the y0, and subtracting that from ix gives a rough idea of how close ix is
to a multiple of pi.  But I'm not sure.

+      // 2nd iteration needed, good to 118
+      t = r;
+      w = n * kTrig[3];
+      r = t - w;
+      w = n * kTrig[4] - ((t - r) - w);
+      y0 = r - w;
+      if (ix - (%_DoubleHi(y0) & 0x7ff00000) > 0x3100000) {

[Raymond Toy, 2014/07/30 19:55:03]

Same comment as for line 231.

[Yang, 2014/08/01 07:29:55]

Acknowledged.

+        // 3rd iteration needed. 151 bits accuracy
+        t = r;
+        w = n * kTrig[5];
+        r = t - w;
+        w = n * kTrig[6] - ((t - r) - w);
+        y0 = r - w;
+      }
+    }
+    y1 = (r - y0) - w;
+    if (hx < 0) {
+      n = -n;
+      y0 = -y0;
+      y1 = -y1;
+    }
+  } else {
+    // Need to do full Payne-Hanek reduction here!
+    var r = %RemPiO2(X);
+    n = r[0];
+    y0 = r[1];
+    y1 = r[2];
+  }
+endmacro
+
+// Sine for [-pi/4, pi/4], pi/4 ~ 0.7854

[Raymond Toy, 2014/07/24 16:51:11]

I think it would be beneficial to add the original fdlibm comment on how this
routine works.

[Yang, 2014/07/28 11:28:59]

Done.

+macro RETURN_KERNELSIN(X0, X1, SIGN)
+  var z = X0 * X0;
+  var v = z * X0;
+  var r = kTrig[8] + z * (kTrig[9] + z * (kTrig[10] +
+                     z * (kTrig[11] + z * kTrig[12])));
+  return (X0 - ((z * (0.5 * X1 - v * r) - X1) - v * kTrig[7])) SIGN;

[Raymond Toy, 2014/07/24 16:51:11]

The original kernel_sin function had a slight optimization where if x1 were
known to be zero, some of the computations could be removed. It doesn't affect
the result to do the extra computations, but there are extra computations.

Is it not needed?

[Yang, 2014/07/28 11:28:59]

I think I checked this. And I rechecked. In neither micro-benchmarks nor popular
JS Benchmarks this buys anything.

[Raymond Toy, 2014/07/30 19:55:04]

Fair enough.

+endmacro

[Raymond Toy, 2014/07/24 16:51:12]

I think readability would be enhanced if you had a separate array for the
coefficients here instead of using kTrig[7-12].  Same goes for the cosine and
tangent coefficients below.

Is it really necessary to use just one array?

[Yang, 2014/07/28 11:28:59]

Acknowledged.

 
-var kInversePiHalf      = 0.636619772367581343;      // 2 / pi
-var kInversePiHalfS26   = 9.48637384723993156e-9;    // 2 / pi / (2^26)
-var kS26                = 1 << 26;
-var kTwoStepThreshold   = 1 << 27;
-// pi / 2 rounded up
-var kPiHalf             = 1.570796326794896780;      // 0x192d4454fb21f93f
-// We use two parts for pi/2 to emulate a higher precision.
-// pi_half_1 only has 26 significant bits for mantissa.
-// Note that pi_half > pi_half_1 + pi_half_2
-var kPiHalf1            = 1.570796325802803040;      // 0x00000054fb21f93f
-var kPiHalf2            = 9.920935796805404252e-10;  // 0x3326a611460b113e
-
-var kSamples;            // Initialized to a number during genesis.
-var kIndexConvert;       // Initialized to kSamples / (pi/2) during genesis.
-var kSinTable;           // Initialized to a Float64Array during genesis.
-var kCosXIntervalTable;  // Initialized to a Float64Array during genesis.
-
-// This implements sine using the following algorithm.
-// 1) Multiplication takes care of to-number conversion.
-// 2) Reduce x to the first quadrant [0, pi/2].
-//    Conveniently enough, in case of +/-Infinity, we get NaN.
-//    Note that we try to use only 26 instead of 52 significant bits for
-//    mantissa to avoid rounding errors when multiplying.  For very large
-//    input we therefore have additional steps.
-// 3) Replace x by (pi/2-x) if x was in the 2nd or 4th quadrant.
-// 4) Do a table lookup for the closest samples to the left and right of x.
-// 5) Find the derivatives at those sampling points by table lookup:
-//    dsin(x)/dx = cos(x) = sin(pi/2-x) for x in [0, pi/2].
-// 6) Use cubic spline interpolation to approximate sin(x).
-// 7) Negate the result if x was in the 3rd or 4th quadrant.
-// 8) Get rid of -0 by adding 0.
-function TrigonometricInterpolation(x, phase) {
-  if (x < 0 || x > kPiHalf) {
-    var multiple;
-    while (x < -kTwoStepThreshold || x > kTwoStepThreshold) {
-      // Let's assume this loop does not terminate.
-      // All numbers x in each loop forms a set S.
-      // (1) abs(x) > 2^27 for all x in S.
-      // (2) abs(multiple) != 0 since (2^27 * inverse_pi_half_s26) > 1
-      // (3) multiple is rounded down in 2^26 steps, so the rounding error is
-      //     at most max(ulp, 2^26).
-      // (4) so for x > 2^27, we subtract at most (1+pi/4)x and at least
-      //     (1-pi/4)x
-      // (5) The subtraction results in x' so that abs(x') <= abs(x)*pi/4.
-      //     Note that this difference cannot be simply rounded off.
-      // Set S cannot exist since (5) violates (1).  Loop must terminate.
-      multiple = MathFloor(x * kInversePiHalfS26) * kS26;
-      x = x - multiple * kPiHalf1 - multiple * kPiHalf2;
+// Cosine for [-pi/4, pi/4], pi/4 ~ 0.7854
+macro RETURN_KERNELCOS(X0, X1, SIGN)
+  var ix = %_DoubleHi(X0) & 0x7fffffff;
+  var z = X0 * X0;
+  var r = z * (kTrig[13] + z * (kTrig[14] + z * (kTrig[15] +

[Raymond Toy, 2014/07/24 16:51:11]

Can we use a separate array for kTrig values here to indicate that these are the
coefficients of the cos polynomial approximation?

[Yang, 2014/07/28 11:28:59]

There is not much readability to be gained here imo. I'd rather clean this up
once we have constant pools for intel platforms as well instead of splitting it
up into smaller constant arrays.

[Raymond Toy, 2014/07/30 19:55:04]

I have no problems with the various random constants placed in kTrig. I think,
however, that readability is vastly improved if the coefficients of the
approximations are put in their own arrays. 

When doing the first pass over this, I had to go hunt around to figure out what
kTrig[13] was. If you had just had C[0] (or is it C[1])? then it would have been
pretty obvious what is happening here.

[Yang, 2014/08/01 07:29:55]

Done with macros.

+          z * (kTrig[16] + z * (kTrig[17] + z * kTrig[18])))));
+  if (ix < 0x3fd33333) {

[Raymond Toy, 2014/07/30 19:55:03]

Add comment that ix < 0x3fd33333 is testing for |x| ~< 0.3, so it's easier to
match the function documentation.

[Yang, 2014/08/01 07:29:55]

Done.

+    return (1 - (0.5 * z - (z * r - X0 * X1))) SIGN;
+  } else {
+    var qx;
+    if (ix > 0x3fe90000) {

[Raymond Toy, 2014/07/24 16:51:12]

I think it's really useful to include the fdlibm comment about how this function
works. Without it, it's really impossible to understand why there's a test here
for 0x3fe90000 and what the code is trying to do with that.

[Yang, 2014/07/28 11:28:59]

imo the comment in your port doesn't explain much either except for pointing to
the original fdlibm. So I think adding that line doesn't add much value.

[Raymond Toy, 2014/07/30 19:55:03]

No, but at least it makes it clear that ix > 0x3fe90000 is testing the case for
x > 0.78125. And that is certainly easier to understand than the former and
makes it easier to match up with the function documentation.

[Yang, 2014/08/01 07:29:55]

Done.

+      qx = 0.28125;
+    } else {
+      qx = %_ConstructDouble(%_DoubleHi(0.25 * X0), 0);
     }
-    multiple = MathFloor(x * kInversePiHalf);
-    x = x - multiple * kPiHalf1 - multiple * kPiHalf2;
-    phase += multiple;
+    var hz = 0.5 * z - qx;
+    return (1 - qx - (hz - (z * r - X0 * X1))) SIGN;
+  }
+endmacro
+
+// Tangent for [-pi/4, pi/4], pi/4 ~ 0.7854

[Raymond Toy, 2014/07/24 16:51:12]

Include original fdlibm comment. Without that it's really hard to figure out
that even though this is called KernelTan, it can return either tan(x) or cot(x)
depending on the value of returnTan.

[Yang, 2014/07/28 11:28:59]

Done.

+function KernelTan(x, y, returnTan) {
+  var z;
+  var w;
+  var hx = %_DoubleHi(x);
+  var ix = hx & 0x7fffffff;
+
+  if (ix < 0x3e300000) {
+    // x < 2^-28
+    if (((ix | %_DoubleLo(x)) | (returnTan + 1)) == 0) {

[Raymond Toy, 2014/07/30 19:55:04]

This might be more obvious if we said abs(x) == 0 and returnTan = -1.  I think
that's what this is trying to test so we can return Infinity for cot(0).

[Yang, 2014/08/01 07:29:55]

I don't think I changed anything here from your port. I added a comment to
clarify.

+      return 1 / MathAbs(x);
+    } else {
+      if (returnTan == 1) {
+        return x;
+      } else {
+        // Compute -1/(x + y) carefully
+        var w = x + y;
+        var z = %_ConstructDouble(%_DoubleHi(w), 0);
+        var v = y - (z - x);
+        var a = -1 / w;
+        var t = %_ConstructDouble(%_DoubleHi(a), 0);
+        var s = 1 + t * z;
+        return t + a * (s + t * v);
+      }
+    }
+  }
+  if (ix >= 0x3fe59429) {
+    // |x| > .6744
+    if (x < 0) {
+      x = -x;
+      y = -y;
+    }
+    z = kTrig[32] - x;
+    w = kTrig[33] - y;
+    x = z + w;
+    y = 0;
+  }
+  z = x * x;
+  w = z * z;
+
+  var r = kTrig[20] + w * (kTrig[22] + w * (kTrig[24] +

[Raymond Toy, 2014/07/24 16:51:11]

Can we use a separate array for the coefficients instead of kTrig?

[Yang, 2014/07/28 11:28:59]

Acknowledged.

[Raymond Toy, 2014/07/30 19:55:04]

Add back the comment. This would be unexpected from the function description. 
The comment was:

    //
    // Break x^5*(T[1]+x^2*T[2]+...) into
    // x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
    // x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))

I do not know why fdlibm does this.

+                      w * (kTrig[26] + w * (kTrig[28] + w * kTrig[30]))));
+  var v = z * (kTrig[21] + w * (kTrig[23] + w * (kTrig[25] +
+                           w * (kTrig[27] + w * (kTrig[29] + w * kTrig[31])))));
+  var s = z * x;
+  r = y + z * (s * (r + v) + y);
+  r = r + kTrig[19] * s;
+  w = x + r;
+  if (ix >= 0x3fe59428) {

[Raymond Toy, 2014/07/30 19:55:04]

Document what 0x3fe59428 is.

+    return (1 - ((hx >> 30) & 2)) *

[Raymond Toy, 2014/07/30 19:55:04]

Describe what ((hx >> 30) & 2) is doing.

+      (returnTan - 2.0 * (x - (w * w / (w + returnTan) - r)));
+  }
+  if (returnTan == 1) {
+    return w;
+  } else {
+    z = %_ConstructDouble(%_DoubleHi(w), 0);
+    v = r - (z - x);
+    var a = -1 / w;
+    var t = %_ConstructDouble(%_DoubleHi(a), 0);
+    s = 1 + t * z;
+    return t + a * (s + t * v);
+  }
+}
+
+function MathSinSlow(x) {

[Raymond Toy, 2014/07/24 16:51:12]

I think MathSinAccurate is a better name. Or maybe MathSinWithPIReduction.

[Yang, 2014/07/28 11:28:59]

This function is not more accurate. The difference is that MathSin handles the
fast case, which is inlined, while this handles the slow case. We have this
naming scheme elsewhere, so I'd like to keep it this way.

[Raymond Toy, 2014/07/30 19:55:03]

Sure, that makes sense.

+  REMPIO2(x);
+  var sign = 1 - (n & 2);
+  if (n & 1) {
+    RETURN_KERNELCOS(y0, y1, * sign);
+  } else {
+    RETURN_KERNELSIN(y0, y1, * sign);
+  }
+}
+
+function MathCosSlow(x) {

[Raymond Toy, 2014/07/24 16:51:11]

Similar comment as MathSinSlow.

[Yang, 2014/07/28 11:28:59]

Acknowledged.

+  REMPIO2(x);
+  if (n & 1) {
+    var sign = (n & 2) - 1;
+    RETURN_KERNELSIN(y0, y1, * sign);
+  } else {
+    var sign = 1 - (n & 2);
+    RETURN_KERNELCOS(y0, y1, * sign);
+  }
+}
+
+// ECMA 262 - 15.8.2.16
+function MathSin(x) {
+  x = x * 1;  // Convert to number.
+  if ((%_DoubleHi(x) & 0x7fffffff) <= 0x3fe921fb) {
+    // |x| < pi/4, approximately.  No reduction needed.
+    if (%_IsMinusZero(x)) return x;

[Raymond Toy, 2014/07/24 18:08:44]

Is this test for -0 necessary? KERNELSIN should be able to handle this.

[Yang, 2014/07/28 11:28:59]

You are right. Removed.

+    RETURN_KERNELSIN(x, 0, /* empty */);
+  }
+  return MathSinSlow(x);
+}
+
+// ECMA 262 - 15.8.2.7
+function MathCos(x) {
+  x = x * 1;  // Convert to number.
+  if ((%_DoubleHi(x) & 0x7fffffff) <= 0x3fe921fb) {
+    // |x| < pi/4, approximately.  No reduction needed.
+    RETURN_KERNELCOS(x, 0, /* empty */);
+  }
+  return MathCosSlow(x);
+}
+
+// ECMA 262 - 15.8.2.18
+function MathTan(x) {
+  x = x * 1;  // Convert to number.
+  if ((%_DoubleHi(x) & 0x7fffffff) <= 0x3fe921fb) {
+    // |x| < pi/4, approximately.  No reduction needed.
+    if (%_IsMinusZero(x)) return x;

[Raymond Toy, 2014/07/24 18:08:43]

Is this test for -0 really necessary? If KernelTan doesn't work for -0, then
KernelTan is probably broken.

[Yang, 2014/07/28 11:28:59]

Done.

+    return KernelTan(x, 0, 1);
   }
-  var double_index = x * kIndexConvert;
-  if (phase & 1) double_index = kSamples - double_index;
-  var index = double_index | 0;
-  var t1 = double_index - index;
-  var t2 = 1 - t1;
-  var y1 = kSinTable[index];
-  var y2 = kSinTable[index + 1];
-  var dy = y2 - y1;
-  return (t2 * y1 + t1 * y2 +
-              t1 * t2 * ((kCosXIntervalTable[index] - dy) * t2 +
-                         (dy - kCosXIntervalTable[index + 1]) * t1))
-         * (1 - (phase & 2)) + 0;
+  REMPIO2(x);
+  return KernelTan(y0, y1, (n & 1) ? -1 : 1);
 }
 
 
@@ -537,8 +689,6 @@ function SetUpMath() {
   %SetInlineBuiltinFlag(MathRandom);
   %SetInlineBuiltinFlag(MathSin);
   %SetInlineBuiltinFlag(MathCos);
-  %SetInlineBuiltinFlag(MathTan);
-  %SetInlineBuiltinFlag(TrigonometricInterpolation);
 }
 
 SetUpMath();