fdlibm: reduce accuracy slightly in favor of performance.

src/math.js

Index: src/math.js
diff --git a/src/math.js b/src/math.js
index dcff5cb287248fa4e477b5f66a934ef96d73df04..19e97bf5ba6fc9bf3d485106fba5341b2992799d 100644
--- a/src/math.js
+++ b/src/math.js
@@ -188,53 +188,18 @@ var pio2_1t = 6.07710050650619224932e-11;
 var pio2_2 = 6.07710050630396597660e-11;
 var pio2_2t = 2.02226624879595063154e-21;
 
-// Table of values of multiples of pi/2 from pi/2 to 50*pi/2. This is
-// used as a quick check to see if an argument is close to a multiple
-// of pi/2 and needs extra bits for reduction.  This array contains
-// the high word the multiple of pi/2.
-
-var npio2_hw = [
-  0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
-  0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
-  0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
-  0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
-  0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
-  0x404858EB, 0x404921FB
-];
-
 macro REMPIO2(x)
-  var n, y0, y1;
+  var n, y;
   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 - pio2_1;
-      if (ix != 0x3ff921fb) {
-        // 33+53 bit pi is good enough
-        y0 = z - pio2_1t;
-        y1 = (z - y0) - pio2_1t;
-      } else {
-        // near pi/2, use 33+33+53 bit pi
-        z -= pio2_2;
-        y0 = z - pio2_2t;
-        y1 = (z - y0) - pio2_2t;
-      }
+      y = x - pio2_1 - pio2_2 - pio2_2t;

[Raymond Toy, 2014/06/02 16:55:53]

I don't think this is computing the same values as before.

What happened to the if statement? The value of y0 is different on the two
branches.  Can you prove that y is the same value as y0 on both branches.

And I think you need to be careful about the order of operations in computing y.
 It's not clear to me that the compiler won't notice that pio2_1, pio2_2, and
pio2_2t are constants that the compiler won't just coalesce into one constant,
thereby destroying all of the extra bits you actually want to use.

       n = 1;
     } else {
-      // Negative x
-      var z = x + pio2_1;
-      if (ix != 0x3ff921fb) {
-        // 33+53 bit pi is good enough
-        y0 = z + pio2_1t;
-        y1 = (z - y0) + pio2_1t;
-      } else {
-        // near pi/2, use 33+33+53 bit pi
-        z += pio2_2;
-        y0 = z + pio2_2t;
-        y1 = (z - y0) + pio2_2t;
-      }
+      y = x + pio2_1 + pio2_2 + pio2_2t;
       n = -1;
     }
   } else if (ix <= 0x413921fb) {
@@ -245,38 +210,28 @@ macro REMPIO2(x)
     var r = t - fn * pio2_1;
     var w = fn * pio2_1t;
     // First round good to 85 bit
-    if (n < 32 && ix != npio2_hw[n - 1]) {
-      // Quick check for cancellation
-      y0 = r - w;
-    } else {
-      var j = ix >> 20;
-      y0 = r - w;
-      var i = j - (( %_DoubleHi(y0) >> 20) & 0x7ff);
-      if (i > 16) {
-        // 2nd iteration needed, good to 118
-        t = r;
-        w = fn * pio2_2;
-        r = t - w;
-        w = fn * pio2_2t - ((t - r) - w);
-        y0 = r - w;
-      }
+    y = r - w;
+    var i = (ix - (%_DoubleHi(y) & 0x7ff00000)) >> 20;
+    if (i > 16) {
+      // 2nd iteration needed, good to 118
+      t = r;
+      w = fn * pio2_2;
+      r = t - w;
+      w = fn * pio2_2t - ((t - r) - w);
+      y = r - w;
     }
-    y1 = (r - y0) - w;
     if (hx < 0) {
       n = -n;
-      y0 = -y0;
-      y1 = -y1;
+      y = -y;
     }
   } else if (ix >= 0x7ff00000) {
     n = 0;
-    y0 = NAN;
-    y1 = NAN;
+    y = NAN;
   } else {
     // Need to do full Payne-Hanek reduction here!
     var r = %RemPiO2(x);
     n = r[0];
-    y0 = r[1];
-    y1 = r[2];
+    y = r[1];
   }
 endmacro
 
@@ -289,11 +244,11 @@ var S4 = 2.75573137070700676789e-06;
 var S5 = -2.50507602534068634195e-08;
 var S6 = 1.58969099521155010221e-10;
 
-function KernelSin(x, y) {
+function KernelSin(x) {
   var z = x * x;
   var v = z * x;
   var r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6)));
-  return x - ((z * (0.5 * y - v * r) - y) - v * S1);
+  return v * (z * r + S1) + x;
 }
 
 // Cosine for [-pi/4, pi/4], pi/4 ~ 0.7854
@@ -305,12 +260,12 @@ var C4 = -2.75573143513906633035e-07;
 var C5 = 2.08757232129817482790e-09;
 var C6 = -1.13596475577881948265e-11;
 
-function KernelCos(x, y) {
+function KernelCos(x) {
   var ix = %_DoubleHi(x) & 0x7fffffff;
   var z = x * x;
   var r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6)))));
   if (ix < 0x3fd33333) {
-    return 1 - (0.5 * z - (z * r - x * y));
+    return (r - 0.5) * z + 1;
   } else {
     var qx;
     if (ix > 0x3fe90000) {
@@ -318,9 +273,7 @@ function KernelCos(x, y) {
     } else {
       qx = %_ConstructDouble(%_DoubleHi(0.25 * x), 0);
     }
-    var hz = 0.5 * z - qx;
-    var a = 1 - qx;
-    return a - (hz - (z * r - x * y));
+    return (1 - qx) - ((0.5 * z - qx) - (z * r));
   }
 }
 
@@ -340,7 +293,7 @@ var T10 = 7.14072491382608190305e-05;
 var T11 = -1.85586374855275456654e-05;
 var T12 = 2.59073051863633712884e-05;
 
-function KernelTan(x, y, returnTan) {
+function KernelTan(x, returnTan) {
   var z;
   var w;
   var hx = %_DoubleHi(x);
@@ -356,9 +309,9 @@ function KernelTan(x, y, returnTan) {
         return x;
       } else {
         // Compute -1/(x + y) carefully
-        var w = x + y;
+        var w = x;
         var z = %_ConstructDouble( %_DoubleHi(w), 0);
-        var v = y - (z - x);
+        var v = x - z;
         var a = -1 / w;
         var t = %_ConstructDouble( %_DoubleHi(a), 0);
         var s = 1 + t * z;
@@ -370,14 +323,12 @@ function KernelTan(x, y, returnTan) {
     // |x| > .6744
     if (x < 0) {
       x = -x;
-      y = -y;
     }
     var pio4 = 7.85398163397448278999e-01;
     var pio4lo = 3.06161699786838301793e-17;
     z = pio4 - x;
-    w = pio4lo - y;
+    w = pio4lo;
     x = z + w;
-    y = 0;
   }
   z = x * x;
   w = z * z;
@@ -389,7 +340,7 @@ function KernelTan(x, y, returnTan) {
   var r = T1 + w * (T3 + w * (T5 + w * (T7 + w * (T9 + w * T11))));
   var v = z * (T2 + w * (T4 + w * (T6 + w * (T8 + w * (T10 + w * T12)))));
   var s = z * x;
-  r = y + z * (s * (r + v) + y);
+  r = z * (s * (r + v));
   r = r + T0 * s;
   w = x + r;
   if (ix >= 0x3fe59428) {
@@ -414,15 +365,15 @@ function MathSinSlow(x) {
   REMPIO2(x);
   if (n & 2) {
     if (n & 1) {
-      return -KernelCos(y0, y1);
+      return -KernelCos(y);
     } else {
-      return -KernelSin(y0, y1);
+      return -KernelSin(y);
     }
   } else {
     if (n & 1) {
-      return KernelCos(y0, y1);
+      return KernelCos(y);
     } else {
-      return KernelSin(y0, y1);
+      return KernelSin(y);
     }
   }
 }
@@ -431,15 +382,15 @@ function MathCosSlow(x) {
   REMPIO2(x);
   if (n & 2) {
     if (n & 1) {
-      return KernelSin(y0, y1);
+      return KernelSin(y);
     } else {
-      return -KernelCos(y0, y1);
+      return -KernelCos(y);
     }
   } else {
     if (n & 1) {
-      return -KernelSin(y0, y1);
+      return -KernelSin(y);
     } else {
-      return KernelCos(y0, y1);
+      return KernelCos(y);
     }
   }
 }
@@ -448,7 +399,7 @@ function MathTanSlow(x) {
   REMPIO2(x);
   // flag is 1 if n is even and -1 if n is odd
   var flag = (n & 1) ? -1 : 1;
-  return KernelTan(y0, y1, flag)
+  return KernelTan(y, flag)
 }
 
 //ECMA 262 - 15.8.2.16
@@ -470,7 +421,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);
+    return KernelCos(x);
   }
   return MathCosSlow(x);
 }
@@ -481,7 +432,7 @@ function MathTan(x) {
   if (%_IsMinusZero(x)) return x;
   if ((%_DoubleHi(x) & 0x7fffffff) <= 0x3fe921fb) {
     // |x| < pi/4, approximately.  No reduction needed.
-    return KernelTan(x, 0, 1);
+    return KernelTan(x, 1);
   }
   return MathTanSlow(x);
 }