From 0f81d7698a394c8f231f439847bde8f831ec4a73 Mon Sep 17 00:00:00 2001 From: "yangguo@chromium.org" Date: Tue, 12 Aug 2014 13:36:33 +0000 Subject: [PATCH] Implement Math.log1p using port from fdlibm. Port contributed by Raymond Toy . 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 | 8 +- src/math.js | 22 +---- test/mjsunit/es6/math-log1p.js | 45 ++++++++-- third_party/fdlibm/fdlibm.cc | 13 ++- third_party/fdlibm/fdlibm.h | 4 +- third_party/fdlibm/fdlibm.js | 192 ++++++++++++++++++++++++++++++++++++---- tools/generate-runtime-tests.py | 2 +- 7 files changed, 234 insertions(+), 52 deletions(-) diff --git a/src/bootstrapper.cc b/src/bootstrapper.cc index 361960d..240be71 100644 --- a/src/bootstrapper.cc +++ b/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 trig_buffer = v8::ArrayBuffer::New( reinterpret_cast(isolate), - const_cast(fdlibm::TrigonometricConstants::constants), - table_num_bytes); + const_cast(fdlibm::MathConstants::constants), table_num_bytes); v8::Local 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(); } diff --git a/src/math.js b/src/math.js index f9139c6..13cdb31 100644 --- a/src/math.js +++ b/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 )); diff --git a/test/mjsunit/es6/math-log1p.js b/test/mjsunit/es6/math-log1p.js index 01f3533..5468444 100644 --- a/test/mjsunit/es6/math-log1p.js +++ b/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)); diff --git a/third_party/fdlibm/fdlibm.cc b/third_party/fdlibm/fdlibm.cc index fe7e231..2f6eab1 100644 --- a/third_party/fdlibm/fdlibm.cc +++ b/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 }; diff --git a/third_party/fdlibm/fdlibm.h b/third_party/fdlibm/fdlibm.h index 39c4b50..7985c3a 100644 --- a/third_party/fdlibm/fdlibm.h +++ b/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 diff --git a/third_party/fdlibm/fdlibm.js b/third_party/fdlibm/fdlibm.js index d5dbb72..a55b7c7 100644 --- a/third_party/fdlibm/fdlibm.js +++ b/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); + } +} diff --git a/tools/generate-runtime-tests.py b/tools/generate-runtime-tests.py index a41df67..0fcd987 100755 --- a/tools/generate-runtime-tests.py +++ b/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. -- 2.7.4