fdlibm: reduce accuracy slightly in favor of performance.

src/math.js

 // Copyright 2012 the V8 project authors. All rights reserved.
 // Use of this source code is governed by a BSD-style license that can be
 // found in the LICENSE file.
 
 "use strict";
 
 // This file relies on the fact that the following declarations have been made
 // in runtime.js:
 // var $Object = global.Object;
 
@@ skipping 170 matching lines @@
 // software is freely granted, provided that this notice
 // is preserved.
 // ====================================================
 
 var invpio2 = 6.36619772367581382433e-01;
 var pio2_1 = 1.57079632673412561417e+00;
 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) {
     // |x| ~<= 2^19*(pi/2), medium size
     var t = MathAbs(x);
     n = MathRound(t * invpio2);
     var fn = n;
     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
 
 // Sine for [-pi/4, pi/4], pi/4 ~ 0.7854
 
 var S1 = -1.66666666666666324348e-01;
 var S2 = 8.33333333332248946124e-03;
 var S3 = -1.98412698298579493134e-04;
 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
 
 var C1 = 4.16666666666666019037e-02;
 var C2 = -1.38888888888741095749e-03;
 var C3 = 2.48015872894767294178e-05;
 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) {
       qx = 0.28125;
     } 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));
   }
 }
 
 // Tangent for [-pi/4, pi/4], pi/4 ~ 0.7854
 
 var T0 = 3.33333333333334091986e-01;
 var T1 = 1.33333333333201242699e-01;
 var T2 = 5.39682539762260521377e-02;
 var T3 = 2.18694882948595424599e-02;
 var T4 = 8.86323982359930005737e-03;
 var T5 = 3.59207910759131235356e-03;
 var T6 = 1.45620945432529025516e-03;
 var T7 = 5.88041240820264096874e-04;
 var T8 = 2.46463134818469906812e-04;
 var T9 = 7.81794442939557092300e-05;
 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);
   var ix = hx & 0x7fffffff;
 
   if (ix < 0x3e300000) {
     // x < 2^-28
     // We don't try to generate inexact.
     if (((ix | %_DoubleLo(x)) | (returnTan + 1)) == 0) {
       return 1 / MathAbs(x);
     } else {
       if (returnTan == 1) {
         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;
         return t + a * (s + t * v);
       }
     }
   }
   if (ix >= 0x3fe59429) {
     // |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;
   //
   // 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]))
   //
   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) {
     return (1 - ((hx >> 30) & 2)) *
       (returnTan - 2.0 * (x - (w * w / (w + returnTan) - r)));
   }
   if (returnTan == 1) {
     return w;
   } else {
     // Compute -1/(x+r) accurately
     z = %_ConstructDouble( %_DoubleHi(w), 0);
     v = r - (z - x); // z+v = r+x
     var a = -1 / w;
     var t = %_ConstructDouble( %_DoubleHi(a), 0);
     s = 1 + t * z;
     return t + a * (s + t * v);
   }
 }
 
 
 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);
     }
   }
 }
 
 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);
     }
   }
 }
 
 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
 function MathSin(x) {
   x = x * 1;  // Convert to number and deal with -0.
   if (%_IsMinusZero(x)) return x;
   if ((%_DoubleHi(x) & 0x7fffffff) <= 0x3fe921fb) {
     // |x| < pi/4, approximately.  No reduction needed.
     var z = x * x;
     var v = z * x;
     var r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6)));
     return x + v * (S1 + z * r);
   }
   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);
+    return KernelCos(x);
   }
   return MathCosSlow(x);
 }
 
 //ECMA 262 - 15.8.2.18
 function MathTan(x) {
   x = x * 1;  // Convert to number and deal with -0.
   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);
 }
 
 // -------------------------------------------------------------------
 
 function SetUpMath() {
   %CheckIsBootstrapping();
 
   %SetPrototype($Math, $Object.prototype);
@@ skipping 44 matching lines @@
   %SetInlineBuiltinFlag(MathRandom);
   %SetInlineBuiltinFlag(MathSin);
   %SetInlineBuiltinFlag(MathCos);
   %SetInlineBuiltinFlag(MathTan);
   %SetInlineBuiltinFlag(KernelSin);
   %SetInlineBuiltinFlag(KernelCos);
   %SetInlineBuiltinFlag(KernelTan);
 }
 
 SetUpMath();