Implement Math.log2 via ported extract from fdlibm.

Adapted from Raymond Toy's (rtoy@chromium.org) port, extracted from fdlibm's pow implementation.

R=rtoy@chromium.org
BUG=v8:3579
LOG=N

Review URL: https://codereview.chromium.org/786823003

Cr-Commit-Position: refs/heads/master@{#25768}

src/math.js

@@ -227,11 +227,6 @@ function MathAtanh(x) {
   return 0.5 * MathLog((1 + x) / (1 - x));
 }
 
-// ES6 draft 09-27-13, section 20.2.2.22.
-function MathLog2(x) {
-  return MathLog(x) * 1.442695040888963407;  // log2(x) = log(x)/log(2).
-}
-
 // ES6 draft 09-27-13, section 20.2.2.17.
 function MathHypot(x, y) {  // Function length is 2.
   // We may want to introduce fast paths for two arguments and when
@@ -364,7 +359,7 @@ function SetUpMath() {
     "acosh", MathAcosh,
     "atanh", MathAtanh,
     "log10", MathLog10,   // implemented by third_party/fdlibm
-    "log2", MathLog2,
+    "log2", MathLog2,     // implemented by third_party/fdlibm
     "hypot", MathHypot,
     "fround", MathFroundJS,
     "clz32", MathClz32,

src/third_party/fdlibm/fdlibm.cc

@@ -63,26 +63,36 @@ const double MathConstants::constants[] = {
     3.06161699786838301793e-17,   // pio4lo    33
     6.93147180369123816490e-01,   // ln2_hi    34
     1.90821492927058770002e-10,   // ln2_lo    35
-    1.80143985094819840000e+16,   // 2^54      36
-    6.666666666666666666e-01,     // 2/3       37
-    6.666666666666735130e-01,     // LP1       38  coefficients for log1p
-    3.999999999940941908e-01,     //           39
-    2.857142874366239149e-01,     //           40
-    2.222219843214978396e-01,     //           41
-    1.818357216161805012e-01,     //           42
-    1.531383769920937332e-01,     //           43
-    1.479819860511658591e-01,     // LP7       44
-    7.09782712893383973096e+02,   //           45  overflow threshold for expm1
-    1.44269504088896338700e+00,   // 1/ln2     46
-    -3.33333333333331316428e-02,  // Q1        47  coefficients for expm1
-    1.58730158725481460165e-03,   //           48
-    -7.93650757867487942473e-05,  //           49
-    4.00821782732936239552e-06,   //           50
-    -2.01099218183624371326e-07,  // Q5        51
-    710.4758600739439,            //           52  overflow threshold sinh, cosh
-    4.34294481903251816668e-01,   // ivln10    53  coefficients for log10
-    3.01029995663611771306e-01,   // log10_2hi 54
-    3.69423907715893078616e-13    // log10_2lo 55
+    6.666666666666666666e-01,     // 2/3       36
+    6.666666666666735130e-01,     // LP1       37  coefficients for log1p
+    3.999999999940941908e-01,     //           38
+    2.857142874366239149e-01,     //           39
+    2.222219843214978396e-01,     //           40
+    1.818357216161805012e-01,     //           41
+    1.531383769920937332e-01,     //           42
+    1.479819860511658591e-01,     // LP7       43
+    7.09782712893383973096e+02,   //           44  overflow threshold for expm1
+    1.44269504088896338700e+00,   // 1/ln2     45
+    -3.33333333333331316428e-02,  // Q1        46  coefficients for expm1
+    1.58730158725481460165e-03,   //           47
+    -7.93650757867487942473e-05,  //           48
+    4.00821782732936239552e-06,   //           49
+    -2.01099218183624371326e-07,  // Q5        50
+    710.4758600739439,            //           51  overflow threshold sinh, cosh
+    4.34294481903251816668e-01,   // ivln10    52  coefficients for log10
+    3.01029995663611771306e-01,   // log10_2hi 53
+    3.69423907715893078616e-13,   // log10_2lo 54
+    5.99999999999994648725e-01,   // L1        55  coefficients for log2
+    4.28571428578550184252e-01,   //           56
+    3.33333329818377432918e-01,   //           57
+    2.72728123808534006489e-01,   //           58
+    2.30660745775561754067e-01,   //           59
+    2.06975017800338417784e-01,   // L6        60
+    9.61796693925975554329e-01,   // cp        61  2/(3*ln(2))
+    9.61796700954437255859e-01,   // cp_h      62
+    -7.02846165095275826516e-09,  // cp_l      63
+    5.84962487220764160156e-01,   // dp_h      64
+    1.35003920212974897128e-08    // dp_l      65
 };
 
 

src/third_party/fdlibm/fdlibm.h

@@ -23,7 +23,7 @@ int rempio2(double x, double* y);
 
 // Constants to be exposed to builtins via Float64Array.
 struct MathConstants {
-  static const double constants[56];
+  static const double constants[66];
 };
 }
 }  // namespace v8::internal

src/third_party/fdlibm/fdlibm.js

@@ -20,6 +20,9 @@
 // and exposed through kMath as typed array. We assume the compiler to convert
 // from decimal to binary accurately enough to produce the intended values.
 // kMath is initialized to a Float64Array during genesis and not writable.
+
+"use strict";
+
 var kMath;
 
 const INVPIO2 = kMath[0];
@@ -424,11 +427,12 @@ function MathTan(x) {
 //
 const LN2_HI    = kMath[34];
 const LN2_LO    = kMath[35];
-const TWO54     = kMath[36];
-const TWO_THIRD = kMath[37];
+const TWO_THIRD = kMath[36];
 macro KLOG1P(x)
-(kMath[38+x])
+(kMath[37+x])
 endmacro
+// 2^54
+const TWO54 = 18014398509481984;
 
 function MathLog1p(x) {
   x = x * 1;  // Convert to number.
@@ -607,10 +611,10 @@ function MathLog1p(x) {
 //      For IEEE double 
 //          if x > 7.09782712893383973096e+02 then expm1(x) overflow
 //
-const KEXPM1_OVERFLOW = kMath[45];
-const INVLN2          = kMath[46];
+const KEXPM1_OVERFLOW = kMath[44];
+const INVLN2          = kMath[45];
 macro KEXPM1(x)
-(kMath[47+x])
+(kMath[46+x])
 endmacro
 
 function MathExpm1(x) {
@@ -730,7 +734,7 @@ function MathExpm1(x) {
 //      sinh(x) is |x| if x is +Infinity, -Infinity, or NaN.
 //      only sinh(0)=0 is exact for finite x.
 //
-const KSINH_OVERFLOW = kMath[52];
+const KSINH_OVERFLOW = kMath[51];
 const TWO_M28 = 3.725290298461914e-9;  // 2^-28, empty lower half
 const LOG_MAXD = 709.7822265625;  // 0x40862e42 00000000, empty lower half
 
@@ -782,7 +786,7 @@ function MathSinh(x) {
 //      cosh(x) is |x| if x is +INF, -INF, or NaN.
 //      only cosh(0)=1 is exact for finite x.
 //
-const KCOSH_OVERFLOW = kMath[52];
+const KCOSH_OVERFLOW = kMath[51];
 
 function MathCosh(x) {
   x = x * 1;  // Convert to number.
@@ -840,9 +844,9 @@ function MathCosh(x) {
 //      log10(10**N) = N  for N=0,1,...,22.
 //
 
-const IVLN10 = kMath[53];
-const LOG10_2HI = kMath[54];
-const LOG10_2LO = kMath[55];
+const IVLN10 = kMath[52];
+const LOG10_2HI = kMath[53];
+const LOG10_2LO = kMath[54];
 
 function MathLog10(x) {
   x = x * 1;  // Convert to number.
@@ -875,3 +879,118 @@ function MathLog10(x) {
   z = y * LOG10_2LO + IVLN10 * MathLog(x);
   return z + y * LOG10_2HI;
 }
+
+
+// ES6 draft 09-27-13, section 20.2.2.22.
+// Return the base 2 logarithm of x
+//
+// fdlibm does not have an explicit log2 function, but fdlibm's pow
+// function does implement an accurate log2 function as part of the
+// pow implementation.  This extracts the core parts of that as a
+// separate log2 function.
+
+// Method:
+// Compute log2(x) in two pieces:
+// log2(x) = w1 + w2
+// where w1 has 53-24 = 29 bits of trailing zeroes.
+
+const DP_H = kMath[64];
+const DP_L = kMath[65];
+
+// Polynomial coefficients for (3/2)*(log2(x) - 2*s - 2/3*s^3)
+macro KLOG2(x)
+(kMath[55+x])
+endmacro
+
+// cp = 2/(3*ln(2)). Note that cp_h + cp_l is cp, but with more accuracy.
+const CP = kMath[61];
+const CP_H = kMath[62];
+const CP_L = kMath[63];
+// 2^53
+const TWO53 = 9007199254740992; 
+
+function MathLog2(x) {
+  x = x * 1;  // Convert to number.
+  var ax = MathAbs(x);
+  var hx = %_DoubleHi(x);
+  var lx = %_DoubleLo(x);
+  var ix = hx & 0x7fffffff;
+
+  // Handle special cases.
+  // log2(+/- 0) = -Infinity
+  if ((ix | lx) == 0) return -INFINITY;
+
+  // log(x) = NaN, if x < 0
+  if (hx < 0) return NAN;
+
+  // log2(Infinity) = Infinity, log2(NaN) = NaN
+  if (ix >= 0x7ff00000) return x;
+
+  var n = 0;
+
+  // Take care of subnormal number.
+  if (ix < 0x00100000) {
+    ax *= TWO53;
+    n -= 53;
+    ix = %_DoubleHi(ax);
+  }
+
+  n += (ix >> 20) - 0x3ff;
+  var j = ix & 0x000fffff;
+
+  // Determine interval.
+  ix = j | 0x3ff00000;  // normalize ix.
+
+  var bp = 1;
+  var dp_h = 0;
+  var dp_l = 0;
+  if (j > 0x3988e) {  // |x| > sqrt(3/2)
+    if (j < 0xbb67a) {  // |x| < sqrt(3)
+      bp = 1.5;
+      dp_h = DP_H;
+      dp_l = DP_L;
+    } else {
+      n += 1;
+      ix -= 0x00100000;
+    }
+  }
+ 
+  ax = %_ConstructDouble(ix, %_DoubleLo(ax));
+
+  // Compute ss = s_h + s_l = (x - 1)/(x+1) or (x - 1.5)/(x + 1.5)
+  var u = ax - bp;
+  var v = 1 / (ax + bp);
+  var ss = u * v;
+  var s_h = %_ConstructDouble(%_DoubleHi(ss), 0);
+
+  // t_h = ax + bp[k] High
+  var t_h = %_ConstructDouble(%_DoubleHi(ax + bp), 0)
+  var t_l = ax - (t_h - bp);
+  var s_l = v * ((u - s_h * t_h) - s_h * t_l);
+
+  // Compute log2(ax)
+  var s2 = ss * ss;
+  var r = s2 * s2 * (KLOG2(0) + s2 * (KLOG2(1) + s2 * (KLOG2(2) + s2 * (
+                     KLOG2(3) + s2 * (KLOG2(4) + s2 * KLOG2(5))))));
+  r += s_l * (s_h + ss);
+  s2  = s_h * s_h;
+  t_h = %_ConstructDouble(%_DoubleHi(3.0 + s2 + r), 0);
+  t_l = r - ((t_h - 3.0) - s2);
+  // u + v = ss * (1 + ...)
+  u = s_h * t_h;
+  v = s_l * t_h + t_l * ss;
+
+  // 2 / (3 * log(2)) * (ss + ...)
+  p_h = %_ConstructDouble(%_DoubleHi(u + v), 0);
+  p_l = v - (p_h - u);
+  z_h = CP_H * p_h;
+  z_l = CP_L * p_h + p_l * CP + dp_l;
+
+  // log2(ax) = (ss + ...) * 2 / (3 * log(2)) = n + dp_h + z_h + z_l
+  var t = n;
+  var t1 = %_ConstructDouble(%_DoubleHi(((z_h + z_l) + dp_h) + t), 0);
+  var t2 = z_l - (((t1 - t) - dp_h) - z_h);
+
+  // t1 + t2 = log2(ax), sum up because we do not care about extra precision.
+  return t1 + t2;
+}

test/mjsunit/es6/math-log2-log10.js

@@ -25,6 +25,8 @@
 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
+// Flags: --allow-natives-syntax
+
 [Math.log10, Math.log2].forEach( function(fun) {
   assertTrue(isNaN(fun(NaN)));
   assertTrue(isNaN(fun(fun)));
@@ -41,8 +43,57 @@
 
 for (var i = -310; i <= 308; i += 0.5) {
   assertEquals(i, Math.log10(Math.pow(10, i)));
-  assertEqualsDelta(i, Math.log2(Math.pow(2, i)), 1E-13);
+  // Square roots are tested below.
+  if (i != -0.5 && i != 0.5) assertEquals(i, Math.log2(Math.pow(2, i)));
 }
 
 // Test denormals.
 assertEquals(-307.77759430519706, Math.log10(1.5 * Math.pow(2, -1023)));
+
+// Test Math.log2(2^k) for -1074 <= k <= 1023.
+var n = -1074;
+// This loop covers n from -1074 to -1043
+for (var lowbits = 1; lowbits <= 0x80000000; lowbits *= 2) {
+  var x = %_ConstructDouble(0, lowbits);
+  assertEquals(n, Math.log2(x));
+  n++;
+}
+// This loop covers n from -1042 to -1023
+for (var hibits = 1; hibits <= 0x80000; hibits *= 2) {
+  var x = %_ConstructDouble(hibits, 0);
+  assertEquals(n, Math.log2(x));
+  n++;
+}
+// The rest of the normal values of 2^n
+var x = 1;
+for (var n = -1022; n <= 1023; ++n) {
+  var x = Math.pow(2, n);
+  assertEquals(n, Math.log2(x));
+}
+
+// Test special values.
+// Expectation isn't exactly 1/2 because Math.SQRT2 isn't exactly sqrt(2).
+assertEquals(0.5000000000000001, Math.log2(Math.SQRT2));
+
+// Expectation isn't exactly -1/2 because Math.SQRT1_2 isn't exactly sqrt(1/2).
+assertEquals(-0.4999999999999999, Math.log2(Math.SQRT1_2));
+
+assertEquals(3.321928094887362, Math.log2(10));
+assertEquals(6.643856189774724, Math.log2(100));
+
+// Test relationships
+x = 1;
+for (var k = 0; k < 1000; ++k) {
+  var y = Math.abs(Math.log2(x) + Math.log2(1/x));
+  assertEqualsDelta(0, y, 1.5e-14);
+  x *= 1.1;
+}
+
+x = Math.pow(2, -100);
+for (var k = 0; k < 1000; ++k) {
+  var y = Math.log2(x);
+  var expected = Math.log(x) / Math.LN2;
+  var err = Math.abs(y - expected) / expected;
+  assertEqualsDelta(0, err, 1e-15);
+  x *= 1.1;
+}