Implement Math.log2 via ported extract from fdlibm.

src/third_party/fdlibm/fdlibm.js

 // The following is adapted from fdlibm (http://www.netlib.org/fdlibm),
 //
 // ====================================================
 // Copyright (C) 1993-2004 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.
 // ====================================================
 //
 // The original source code covered by the above license above has been
 // modified significantly by Google Inc.
 // Copyright 2014 the V8 project authors. All rights reserved.
 //
 // The following is a straightforward translation of fdlibm routines
 // by Raymond Toy (rtoy@google.com).
 
 // Double constants that do not have empty lower 32 bits are found in fdlibm.cc
 // 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];
 const PIO2_1  = kMath[1];
 const PIO2_1T = kMath[2];
 const PIO2_2  = kMath[3];
 const PIO2_2T = kMath[4];
 const PIO2_3  = kMath[5];
 const PIO2_3T = kMath[6];
 const PIO4    = kMath[32];
@@ skipping 835 matching lines @@
 
   k += (hx >> 20) - 1023;
   i = (k & 0x80000000) >> 31;
   hx = (hx & 0x000fffff) | ((0x3ff - i) << 20);
   y = k + i;
   x = %_ConstructDouble(hx, lx);
 
   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[65];
+const DP_L = kMath[66];
+
+// Polynomial coefficients for (3/2)*(log2(x) - 2*s - 2/3*s^3)
+macro KLOG2(x)
+(kMath[56+x])
+endmacro
+
+// cp = 2/(3*ln(2)). Note that cp_h + cp_l is cp, but with more accuracy.
+const CP = kMath[62];
+const CP_H = kMath[63];
+const CP_L = kMath[64];
+
+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 *= TWO54;
+    n -= 54;

[Raymond Toy, 2014/12/09 20:37:38]

Why 54? The original code multiplied by 2^53.  I suspect this is probably ok,
though, but would prefer 53 since that's what the fdlibm pow code did.

[Yang, 2014/12/10 10:01:41]

Done.

+    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;
+  var th_offset = 0x00080000;
+  if (j > 0x3988e) {  // |x| > sqrt(3/2)
+    if (j < 0xbb67a) {  // |x| < sqrt(3)
+      bp = 1.5;
+      dp_h = DP_H;
+      dp_l = DP_L;
+      th_offset = 0x000c0000;
+    } 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(((ix >> 1) | 0x20000000) + th_offset, 0);

[Raymond Toy, 2014/12/09 20:37:38]

If the comment on line 965 is correct, and I believe it is, then line 966 can be
replaced by

var t_h = %_ConstructDouble(%_DoubleHi(ax + bp), 0)

I ran the tests and this passes.  This is certainly a lot clearer than the
current bit-banging code, and th_offset can be removed.

[Yang, 2014/12/10 10:01:41]

Done.

+  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;
+}