[builtins] Introduce proper Float64Log1p operator.

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.
 // ====================================================
@@ skipping 379 matching lines @@
 function MathTan(x) {
   x = x * 1;  // Convert to number.
   if ((%_DoubleHi(x) & 0x7fffffff) <= 0x3fe921fb) {
     // |x| < pi/4, approximately.  No reduction needed.
     return KernelTan(x, 0, 1);
   }
   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:
-//      Constants are found in fdlibm.cc. We assume the C++ compiler to convert
-//      from decimal to binary accurately enough to produce the intended values.
-//
-// 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.
-//
 define LN2_HI    = 6.93147180369123816490e-01;
 define LN2_LO    = 1.90821492927058770002e-10;
-define TWO_THIRD = 6.666666666666666666e-01;
-define LP1 = 6.666666666666735130e-01;
-define LP2 = 3.999999999940941908e-01;
-define LP3 = 2.857142874366239149e-01;
-define LP4 = 2.222219843214978396e-01;
-define LP5 = 1.818357216161805012e-01;
-define LP6 = 1.531383769920937332e-01;
-define LP7 = 1.479819860511658591e-01;
 
 // 2^54
 define TWO54 = 18014398509481984;
 
-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 * (LP1 + z * (LP2 + z * (LP3 + z * (LP4 +
-          z * (LP5 + z * (LP6 + z * LP7))))));
-  if (k === 0) {
-    return f - (hfsq - s * (hfsq + R));
-  } else {
-    return k * LN2_HI - ((hfsq - (s * (hfsq + R) + (k * LN2_LO + c))) - f);
-  }
-}
-
 // ES6 draft 09-27-13, section 20.2.2.14.
 // Math.expm1
 // Returns exp(x)-1, the exponential of x minus 1.
 //
 // Method
 //   1. Argument reduction:
 //      Given x, find r and integer k such that
 //
 //               x = k*ln2 + r,  |r| <= 0.5*ln2 ~ 0.34658  
 //
@@ skipping 526 matching lines @@
 
 utils.InstallFunctions(GlobalMath, DONT_ENUM, [
   "cos", MathCos,
   "sin", MathSin,
   "tan", MathTan,
   "sinh", MathSinh,
   "cosh", MathCosh,
   "tanh", MathTanh,
   "log10", MathLog10,
   "log2", MathLog2,
-  "log1p", MathLog1p,
   "expm1", MathExpm1
 ]);
 
 %SetForceInlineFlag(MathSin);
 %SetForceInlineFlag(MathCos);
 
 })