Implement Math.log1p using port from fdlibm.

Port contributed by Raymond Toy <rtoy@google.com>.

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

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

git-svn-id: https://v8.googlecode.com/svn/branches/bleeding_edge@23082 ce2b1a6d-e550-0410-aec6-3dcde31c8c00

src/bootstrapper.cc

@@ -2655,19 +2655,17 @@ Genesis::Genesis(Isolate* isolate,
                                   NONE).Assert();
 
     // Initialize trigonometric lookup tables and constants.
-    const int constants_size =
-        ARRAY_SIZE(fdlibm::TrigonometricConstants::constants);
+    const int constants_size = ARRAY_SIZE(fdlibm::MathConstants::constants);
     const int table_num_bytes = constants_size * kDoubleSize;
     v8::Local<v8::ArrayBuffer> trig_buffer = v8::ArrayBuffer::New(
         reinterpret_cast<v8::Isolate*>(isolate),
-        const_cast<double*>(fdlibm::TrigonometricConstants::constants),
-        table_num_bytes);
+        const_cast<double*>(fdlibm::MathConstants::constants), table_num_bytes);
     v8::Local<v8::Float64Array> trig_table =
         v8::Float64Array::New(trig_buffer, 0, constants_size);
 
     Runtime::DefineObjectProperty(
         builtins,
-        factory()->InternalizeOneByteString(STATIC_ASCII_VECTOR("kTrig")),
+        factory()->InternalizeOneByteString(STATIC_ASCII_VECTOR("kMath")),
         Utils::OpenHandle(*trig_table), NONE).Assert();
   }
 

src/math.js

@@ -347,26 +347,6 @@ function MathExpm1(x) {
   }
 }
 
-// ES6 draft 09-27-13, section 20.2.2.20.
-// Use Taylor series to approximate. With y = x + 1;
-// log(y) at 1 == log(1) + log'(1)(y-1)/1! + log''(1)(y-1)^2/2! + ...
-//             == 0 + x - x^2/2 + x^3/3 ...
-// The closer x is to 0, the fewer terms are required.
-function MathLog1p(x) {
-  if (!IS_NUMBER(x)) x = NonNumberToNumber(x);
-  var xabs = MathAbs(x);
-  if (xabs < 1E-7) {
-    return x * (1 - x * (1/2));
-  } else if (xabs < 3E-5) {
-    return x * (1 - x * (1/2 - x * (1/3)));
-  } else if (xabs < 7E-3) {
-    return x * (1 - x * (1/2 - x * (1/3 - x * (1/4 -
-           x * (1/5 - x * (1/6 - x * (1/7)))))));
-  } else {  // Use regular log if not close enough to 0.
-    return MathLog(1 + x);
-  }
-}
-
 // -------------------------------------------------------------------
 
 function SetUpMath() {
@@ -428,7 +408,7 @@ function SetUpMath() {
     "fround", MathFroundJS,
     "clz32", MathClz32,
     "cbrt", MathCbrt,
-    "log1p", MathLog1p,
+    "log1p", MathLog1p,    // implemented by third_party/fdlibm
     "expm1", MathExpm1
   ));
 

test/mjsunit/es6/math-log1p.js

@@ -6,16 +6,16 @@ assertTrue(isNaN(Math.log1p(NaN)));
 assertTrue(isNaN(Math.log1p(function() {})));
 assertTrue(isNaN(Math.log1p({ toString: function() { return NaN; } })));
 assertTrue(isNaN(Math.log1p({ valueOf: function() { return "abc"; } })));
-assertEquals("Infinity", String(1/Math.log1p(0)));
-assertEquals("-Infinity", String(1/Math.log1p(-0)));
-assertEquals("Infinity", String(Math.log1p(Infinity)));
-assertEquals("-Infinity", String(Math.log1p(-1)));
+assertEquals(Infinity, 1/Math.log1p(0));
+assertEquals(-Infinity, 1/Math.log1p(-0));
+assertEquals(Infinity, Math.log1p(Infinity));
+assertEquals(-Infinity, Math.log1p(-1));
 assertTrue(isNaN(Math.log1p(-2)));
 assertTrue(isNaN(Math.log1p(-Infinity)));
 
-for (var x = 1E300; x > 1E-1; x *= 0.8) {
+for (var x = 1E300; x > 1E16; x *= 0.8) {
   var expected = Math.log(x + 1);
-  assertEqualsDelta(expected, Math.log1p(x), expected * 1E-14);
+  assertEqualsDelta(expected, Math.log1p(x), expected * 1E-16);
 }
 
 // Values close to 0:
@@ -35,5 +35,36 @@ function log1p(x) {
 
 for (var x = 1E-1; x > 1E-300; x *= 0.8) {
   var expected = log1p(x);
-  assertEqualsDelta(expected, Math.log1p(x), expected * 1E-14);
+  assertEqualsDelta(expected, Math.log1p(x), expected * 1E-16);
 }
+
+// Issue 3481.
+assertEquals(6.9756137364252422e-03,
+             Math.log1p(8070450532247929/Math.pow(2,60)));
+
+// Tests related to the fdlibm implementation.
+// Test largest double value.
+assertEquals(709.782712893384, Math.log1p(1.7976931348623157e308));
+// Test small values.
+assertEquals(Math.pow(2, -55), Math.log1p(Math.pow(2, -55)));
+assertEquals(9.313225741817976e-10, Math.log1p(Math.pow(2, -30)));
+// Cover various code paths.
+// -.2929 < x < .41422, k = 0
+assertEquals(-0.2876820724517809, Math.log1p(-0.25));
+assertEquals(0.22314355131420976, Math.log1p(0.25));
+// 0.41422 < x < 9.007e15
+assertEquals(2.3978952727983707, Math.log1p(10));
+// x > 9.007e15
+assertEquals(36.841361487904734, Math.log1p(10e15));
+// Normalize u.
+assertEquals(37.08337388996168, Math.log1p(12738099905822720));
+// Normalize u/2.
+assertEquals(37.08336444902049, Math.log1p(12737979646738432));
+// |f| = 0, k != 0
+assertEquals(1.3862943611198906, Math.log1p(3));
+// |f| != 0, k != 0
+assertEquals(1.3862945995384413, Math.log1p(3 + Math.pow(2,-20)));
+// final if-clause: k = 0
+assertEquals(0.5596157879354227, Math.log1p(0.75));
+// final if-clause: k != 0
+assertEquals(0.8109302162163288, Math.log1p(1.25));

third_party/fdlibm/fdlibm.cc

@@ -26,7 +26,7 @@ namespace fdlibm {
 inline double scalbn(double x, int y) { return _scalb(x, y); }
 #endif  // _MSC_VER
 
-const double TrigonometricConstants::constants[] = {
+const double MathConstants::constants[] = {
     6.36619772367581382433e-01,   // invpio2   0
     1.57079632673412561417e+00,   // pio2_1    1
     6.07710050650619224932e-11,   // pio2_1t   2
@@ -61,6 +61,17 @@ const double TrigonometricConstants::constants[] = {
     2.59073051863633712884e-05,   // T12      31
     7.85398163397448278999e-01,   // pio4     32
     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
+    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
 };
 
 

third_party/fdlibm/fdlibm.h

@@ -22,8 +22,8 @@ namespace fdlibm {
 int rempio2(double x, double* y);
 
 // Constants to be exposed to builtins via Float64Array.
-struct TrigonometricConstants {
-  static const double constants[34];
+struct MathConstants {
+  static const double constants[45];
 };
 }
 }  // namespace v8::internal

third_party/fdlibm/fdlibm.js

@@ -13,21 +13,21 @@
 // modified significantly by Google Inc.
 // Copyright 2014 the V8 project authors. All rights reserved.
 //
-// The following is a straightforward translation of fdlibm routines for
-// sin, cos, and tan, by Raymond Toy (rtoy@google.com).
+// The following is a straightforward translation of fdlibm routines
+// by Raymond Toy (rtoy@google.com).
 
 
-var kTrig;  // Initialized to a Float64Array during genesis and is not writable.
+var kMath;  // Initialized to a Float64Array during genesis and is not writable.
 
-const INVPIO2 = kTrig[0];
-const PIO2_1  = kTrig[1];
-const PIO2_1T = kTrig[2];
-const PIO2_2  = kTrig[3];
-const PIO2_2T = kTrig[4];
-const PIO2_3  = kTrig[5];
-const PIO2_3T = kTrig[6];
-const PIO4    = kTrig[32];
-const PIO4LO  = kTrig[33];
+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];
+const PIO4LO  = kMath[33];
 
 // Compute k and r such that x - k*pi/2 = r where |r| < pi/4. For
 // precision, r is returned as two values y0 and y1 such that r = y0 + y1
@@ -133,7 +133,7 @@ endmacro
 //          sin(x) = X + (S1*X + (X *(r-Y/2)+Y))
 //
 macro KSIN(x)
-kTrig[7+x]
+kMath[7+x]
 endmacro
 
 macro RETURN_KERNELSIN(X, Y, SIGN)
@@ -177,7 +177,7 @@ endmacro
 //     thus, reducing the rounding error in the subtraction.
 //
 macro KCOS(x)
-kTrig[13+x]
+kMath[13+x]
 endmacro
 
 macro RETURN_KERNELCOS(X, Y, SIGN)
@@ -199,6 +199,7 @@ macro RETURN_KERNELCOS(X, Y, SIGN)
   }
 endmacro
 
+
 // kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
 // Input x is assumed to be bounded by ~pi/4 in magnitude.
 // Input y is the tail of x.
@@ -235,7 +236,7 @@ endmacro
 // and will cause incorrect results.
 //
 macro KTAN(x)
-kTrig[19+x]
+kMath[19+x]
 endmacro
 
 function KernelTan(x, y, returnTan) {
@@ -354,3 +355,164 @@ function MathTan(x) {
   REMPIO2(x);
   return KernelTan(y0, y1, (n & 1) ? -1 : 1);
 }
+
+// ES6 draft 09-27-13, section 20.2.2.20.
+// Math.log1p
+//
+// Method :                  
+//   1. Argument Reduction: find k and f such that 
+//                      1+x = 2^k * (1+f), 
+//         where  sqrt(2)/2 < 1+f < sqrt(2) .
+//
+//      Note. If k=0, then f=x is exact. However, if k!=0, then f
+//      may not be representable exactly. In that case, a correction
+//      term is need. Let u=1+x rounded. Let c = (1+x)-u, then
+//      log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u),
+//      and add back the correction term c/u.
+//      (Note: when x > 2**53, one can simply return log(x))
+//
+//   2. Approximation of log1p(f).
+//      Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+//            = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+//            = 2s + s*R
+//      We use a special Reme algorithm on [0,0.1716] to generate 
+//      a polynomial of degree 14 to approximate R The maximum error 
+//      of this polynomial approximation is bounded by 2**-58.45. In
+//      other words,
+//                      2      4      6      8      10      12      14
+//          R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s  +Lp6*s  +Lp7*s
+//      (the values of Lp1 to Lp7 are listed in the program)
+//      and
+//          |      2          14          |     -58.45
+//          | Lp1*s +...+Lp7*s    -  R(z) | <= 2 
+//          |                             |
+//      Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+//      In order to guarantee error in log below 1ulp, we compute log
+//      by
+//              log1p(f) = f - (hfsq - s*(hfsq+R)).
+//
+//      3. Finally, log1p(x) = k*ln2 + log1p(f).  
+//                           = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
+//         Here ln2 is split into two floating point number: 
+//                      ln2_hi + ln2_lo,
+//         where n*ln2_hi is always exact for |n| < 2000.
+//
+// Special cases:
+//      log1p(x) is NaN with signal if x < -1 (including -INF) ; 
+//      log1p(+INF) is +INF; log1p(-1) is -INF with signal;
+//      log1p(NaN) is that NaN with no signal.
+//
+// Accuracy:
+//      according to an error analysis, the error is always less than
+//      1 ulp (unit in the last place).
+//
+// Constants:
+// The hexadecimal values are the intended ones for the following 
+// constants. The decimal values may be used, provided that the 
+// compiler will convert from decimal to binary accurately enough 
+// to produce the hexadecimal values shown.
+//
+// Note: Assuming log() return accurate answer, the following
+//       algorithm can be used to compute log1p(x) to within a few ULP:
+//
+//              u = 1+x;
+//              if (u==1.0) return x ; else
+//                          return log(u)*(x/(u-1.0));
+//
+//       See HP-15C Advanced Functions Handbook, p.193.
+//
+const LN2_HI    = kMath[34];
+const LN2_LO    = kMath[35];
+const TWO54     = kMath[36];
+const TWO_THIRD = kMath[37];
+macro KLOGP1(x)
+(kMath[38+x])
+endmacro
+
+function MathLog1p(x) {
+  x = x * 1;  // Convert to number.
+  var hx = %_DoubleHi(x);
+  var ax = hx & 0x7fffffff;
+  var k = 1;
+  var f = x;
+  var hu = 1;
+  var c = 0;
+  var u = x;
+
+  if (hx < 0x3fda827a) {
+    // x < 0.41422
+    if (ax >= 0x3ff00000) {  // |x| >= 1
+      if (x === -1) {
+        return -INFINITY;  // log1p(-1) = -inf
+      } else {
+        return NAN;  // log1p(x<-1) = NaN
+      }
+    } else if (ax < 0x3c900000)  {
+      // For |x| < 2^-54 we can return x.
+      return x;
+    } else if (ax < 0x3e200000) {
+      // For |x| < 2^-29 we can use a simple two-term Taylor series.
+      return x - x * x * 0.5;
+    }
+
+    if ((hx > 0) || (hx <= -0x402D413D)) {  // (int) 0xbfd2bec3 = -0x402d413d
+      // -.2929 < x < 0.41422
+      k = 0;
+    }
+  }
+
+  // Handle Infinity and NAN
+  if (hx >= 0x7ff00000) return x;
+
+  if (k !== 0) {
+    if (hx < 0x43400000) {
+      // x < 2^53
+      u = 1 + x;
+      hu = %_DoubleHi(u);
+      k = (hu >> 20) - 1023;
+      c = (k > 0) ? 1 - (u - x) : x - (u - 1);
+      c = c / u;
+    } else {
+      hu = %_DoubleHi(u);
+      k = (hu >> 20) - 1023;
+    }
+    hu = hu & 0xfffff;
+    if (hu < 0x6a09e) {
+      u = %_ConstructDouble(hu | 0x3ff00000, %_DoubleLo(u));  // Normalize u.
+    } else {
+      ++k;
+      u = %_ConstructDouble(hu | 0x3fe00000, %_DoubleLo(u));  // Normalize u/2.
+      hu = (0x00100000 - hu) >> 2;
+    }
+    f = u - 1;
+  }
+
+  var hfsq = 0.5 * f * f;
+  if (hu === 0) {
+    // |f| < 2^-20;
+    if (f === 0) {
+      if (k === 0) {
+        return 0.0;
+      } else {
+        return k * LN2_HI + (c + k * LN2_LO);
+      }
+    }
+    var R = hfsq * (1 - TWO_THIRD * f);
+    if (k === 0) {
+      return f - R;
+    } else {
+      return k * LN2_HI - ((R - (k * LN2_LO + c)) - f);
+    }
+  }
+
+  var s = f / (2 + f); 
+  var z = s * s;
+  var R = z * (KLOGP1(0) + z * (KLOGP1(1) + z *
+              (KLOGP1(2) + z * (KLOGP1(3) + z *
+              (KLOGP1(4) + z * (KLOGP1(5) + z * KLOGP1(6)))))));
+  if (k === 0) {
+    return f - (hfsq - s * (hfsq + R));
+  } else {
+    return k * LN2_HI - ((hfsq - (s * (hfsq + R) + (k * LN2_LO + c))) - f);
+  }
+}

tools/generate-runtime-tests.py

@@ -51,7 +51,7 @@ EXPECTED_FUNCTION_COUNT = 428
 EXPECTED_FUZZABLE_COUNT = 331
 EXPECTED_CCTEST_COUNT = 7
 EXPECTED_UNKNOWN_COUNT = 16
-EXPECTED_BUILTINS_COUNT = 809
+EXPECTED_BUILTINS_COUNT = 808
 
 
 # Don't call these at all.