Skip to content

V
V8
Commit 930e5ccc authored by yangguo@chromium.org's avatar yangguo@chromium.org
Browse files
Options

Implement Math.expm1 using port from fdlibm. 

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

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

git-svn-id: https://v8.googlecode.com/svn/branches/bleeding_edge@23238 ce2b1a6d-e550-0410-aec6-3dcde31c8c00
parent 63134745
Hide whitespace changes
Inline Side-by-side
Showing with 269 additions and 48 deletions
src/math.js
View file @ 930e5ccc
......@@ -327,26 +327,6 @@ function CubeRoot(x) {
return NEWTON_ITERATION_CBRT(x, approx);
}
// ES6 draft 09-27-13, section 20.2.2.14.
// Use Taylor series to approximate.
// exp(x) - 1 at 0 == -1 + exp(0) + exp'(0)*x/1! + exp''(0)*x^2/2! + ...
// == x/1! + x^2/2! + x^3/3! + ...
// The closer x is to 0, the fewer terms are required.
function MathExpm1(x) {
if (!IS_NUMBER(x)) x = NonNumberToNumber(x);
var xabs = MathAbs(x);
if (xabs < 2E-7) {
return x * (1 + x * (1/2));
} else if (xabs < 6E-5) {
return x * (1 + x * (1/2 + x * (1/6)));
} else if (xabs < 2E-2) {
return x * (1 + x * (1/2 + x * (1/6 +
x * (1/24 + x * (1/120 + x * (1/720))))));
} else { // Use regular exp if not close enough to 0.
return MathExp(x) - 1;
}
}
// -------------------------------------------------------------------
function SetUpMath() {
......@@ -408,8 +388,8 @@ function SetUpMath() {
"fround", MathFroundJS,
"clz32", MathClz32,
"cbrt", MathCbrt,
"log1p", MathLog1p, // implemented by third_party/fdlibm
"expm1", MathExpm1
"log1p", MathLog1p, // implemented by third_party/fdlibm
"expm1", MathExpm1 // implemented by third_party/fdlibm
));
%SetInlineBuiltinFlag(MathCeil);
......
test/mjsunit/es6/math-expm1.js
View file @ 930e5ccc
......@@ -8,19 +8,22 @@ assertTrue(isNaN(Math.expm1(NaN)));
assertTrue(isNaN(Math.expm1(function() {})));
assertTrue(isNaN(Math.expm1({ toString: function() { return NaN; } })));
assertTrue(isNaN(Math.expm1({ valueOf: function() { return "abc"; } })));
assertEquals("Infinity", String(1/Math.expm1(0)));
assertEquals("-Infinity", String(1/Math.expm1(-0)));
assertEquals("Infinity", String(Math.expm1(Infinity)));
assertEquals(Infinity, 1/Math.expm1(0));
assertEquals(-Infinity, 1/Math.expm1(-0));
assertEquals(Infinity, Math.expm1(Infinity));
assertEquals(-1, Math.expm1(-Infinity));
for (var x = 0.1; x < 700; x += 0.1) {
// Sanity check:
// Math.expm1(x) stays reasonably close to Math.exp(x) - 1 for large values.
for (var x = 1; x < 700; x += 0.25) {
var expected = Math.exp(x) - 1;
assertEqualsDelta(expected, Math.expm1(x), expected * 1E-14);
assertEqualsDelta(expected, Math.expm1(x), expected * 1E-15);
expected = Math.exp(-x) - 1;
assertEqualsDelta(expected, Math.expm1(-x), -expected * 1E-14);
assertEqualsDelta(expected, Math.expm1(-x), -expected * 1E-15);
}
// Values close to 0:
// Approximation for values close to 0:
// Use six terms of Taylor expansion at 0 for exp(x) as test expectation:
// exp(x) - 1 == exp(0) + exp(0) * x + x * x / 2 + ... - 1
// == x + x * x / 2 + x * x * x / 6 + ...
......@@ -32,7 +35,44 @@ function expm1(x) {
1/362880 + x * (1/3628800))))))))));
}
// Sanity check:
// Math.expm1(x) stays reasonabliy close to the Taylor series for small values.
for (var x = 1E-1; x > 1E-300; x *= 0.8) {
var expected = expm1(x);
assertEqualsDelta(expected, Math.expm1(x), expected * 1E-14);
assertEqualsDelta(expected, Math.expm1(x), expected * 1E-15);
}
// Tests related to the fdlibm implementation.
// Test overflow.
assertEquals(Infinity, Math.expm1(709.8));
// Test largest double value.
assertEquals(Infinity, Math.exp(1.7976931348623157e308));
// Cover various code paths.
assertEquals(-1, Math.expm1(-56 * Math.LN2));
assertEquals(-1, Math.expm1(-50));
// Test most negative double value.
assertEquals(-1, Math.expm1(-1.7976931348623157e308));
// Test argument reduction.
// Cases for 0.5*log(2) < |x| < 1.5*log(2).
assertEquals(Math.E - 1, Math.expm1(1));
assertEquals(1/Math.E - 1, Math.expm1(-1));
// Cases for 1.5*log(2) < |x|.
assertEquals(6.38905609893065, Math.expm1(2));
assertEquals(-0.8646647167633873, Math.expm1(-2));
// Cases where Math.expm1(x) = x.
assertEquals(0, Math.expm1(0));
assertEquals(Math.pow(2,-55), Math.expm1(Math.pow(2,-55)));
// Tests for the case where argument reduction has x in the primary range.
// Test branch for k = 0.
assertEquals(0.18920711500272105, Math.expm1(0.25 * Math.LN2));
// Test branch for k = -1.
assertEquals(-0.5, Math.expm1(-Math.LN2));
// Test branch for k = 1.
assertEquals(1, Math.expm1(Math.LN2));
// Test branch for k <= -2 || k > 56. k = -3.
assertEquals(1.4411518807585582e17, Math.expm1(57 * Math.LN2));
// Test last branch for k < 20, k = 19.
assertEquals(524286.99999999994, Math.expm1(19 * Math.LN2));
// Test the else branch, k = 20.
assertEquals(1048575, Math.expm1(20 * Math.LN2));
third_party/fdlibm/LICENSE
View file @ 930e5ccc
Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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
......
third_party/fdlibm/fdlibm.cc
View file @ 930e5ccc
......@@ -34,19 +34,19 @@ const double MathConstants::constants[] = {
2.02226624879595063154e-21, // pio2_2t 4
2.02226624871116645580e-21, // pio2_3 5
8.47842766036889956997e-32, // pio2_3t 6
-1.66666666666666324348e-01, // S1 7
-1.66666666666666324348e-01, // S1 7 coefficients for sin
8.33333333332248946124e-03, // 8
-1.98412698298579493134e-04, // 9
2.75573137070700676789e-06, // 10
-2.50507602534068634195e-08, // 11
1.58969099521155010221e-10, // S6 12
4.16666666666666019037e-02, // C1 13
4.16666666666666019037e-02, // C1 13 coefficients for cos
-1.38888888888741095749e-03, // 14
2.48015872894767294178e-05, // 15
-2.75573143513906633035e-07, // 16
2.08757232129817482790e-09, // 17
-1.13596475577881948265e-11, // C6 18
3.33333333333334091986e-01, // T0 19
3.33333333333334091986e-01, // T0 19 coefficients for tan
1.33333333333201242699e-01, // 20
5.39682539762260521377e-02, // 21
2.18694882948595424599e-02, // 22
......@@ -65,13 +65,20 @@ const double MathConstants::constants[] = {
1.90821492927058770002e-10, // ln2_lo 35
1.80143985094819840000e+16, // 2^54 36
6.666666666666666666e-01, // 2/3 37
6.666666666666735130e-01, // LP1 38
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
};
......
third_party/fdlibm/fdlibm.h
View file @ 930e5ccc
......@@ -23,7 +23,7 @@ int rempio2(double x, double* y);
// Constants to be exposed to builtins via Float64Array.
struct MathConstants {
static const double constants[45];
static const double constants[52];
};
}
} // namespace v8::internal
......
third_party/fdlibm/fdlibm.js
View file @ 930e5ccc
// The following is adapted from fdlibm (http://www.netlib.org/fdlibm),
//
// ====================================================
// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
// 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
......@@ -16,8 +16,11 @@
// The following is a straightforward translation of fdlibm routines
// by Raymond Toy (rtoy@google.com).
var kMath; // Initialized to a Float64Array during genesis and is not writable.
// 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.
var kMath;
const INVPIO2 = kMath[0];
const PIO2_1 = kMath[1];
......@@ -407,10 +410,8 @@ function MathTan(x) {
// 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.
// 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:
......@@ -425,7 +426,7 @@ const LN2_HI = kMath[34];
const LN2_LO = kMath[35];
const TWO54 = kMath[36];
const TWO_THIRD = kMath[37];
macro KLOGP1(x)
macro KLOG1P(x)
(kMath[38+x])
endmacro
......@@ -507,12 +508,205 @@ function MathLog1p(x) {
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)))))));
var R = z * (KLOG1P(0) + z * (KLOG1P(1) + z *
(KLOG1P(2) + z * (KLOG1P(3) + z *
(KLOG1P(4) + z * (KLOG1P(5) + z * KLOG1P(6)))))));
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
//
// Here a correction term c will be computed to compensate
// the error in r when rounded to a floating-point number.
//
// 2. Approximating expm1(r) by a special rational function on
// the interval [0,0.34658]:
// Since
// r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ...
// we define R1(r*r) by
// r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r)
// That is,
// R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
// = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
// = 1 - r^2/60 + r^4/2520 - r^6/100800 + ...
// We use a special Remes algorithm on [0,0.347] to generate
// a polynomial of degree 5 in r*r to approximate R1. The
// maximum error of this polynomial approximation is bounded
// by 2**-61. In other words,
// R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5
// where Q1 = -1.6666666666666567384E-2,
// Q2 = 3.9682539681370365873E-4,
// Q3 = -9.9206344733435987357E-6,
// Q4 = 2.5051361420808517002E-7,
// Q5 = -6.2843505682382617102E-9;
// (where z=r*r, and the values of Q1 to Q5 are listed below)
// with error bounded by
// | 5 | -61
// | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2
// | |
//
// expm1(r) = exp(r)-1 is then computed by the following
// specific way which minimize the accumulation rounding error:
// 2 3
// r r [ 3 - (R1 + R1*r/2) ]
// expm1(r) = r + --- + --- * [--------------------]
// 2 2 [ 6 - r*(3 - R1*r/2) ]
//
// To compensate the error in the argument reduction, we use
// expm1(r+c) = expm1(r) + c + expm1(r)*c
// ~ expm1(r) + c + r*c
// Thus c+r*c will be added in as the correction terms for
// expm1(r+c). Now rearrange the term to avoid optimization
// screw up:
// ( 2 2 )
// ({ ( r [ R1 - (3 - R1*r/2) ] ) } r )
// expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- )
// ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 )
// ( )
//
// = r - E
// 3. Scale back to obtain expm1(x):
// From step 1, we have
// expm1(x) = either 2^k*[expm1(r)+1] - 1
// = or 2^k*[expm1(r) + (1-2^-k)]
// 4. Implementation notes:
// (A). To save one multiplication, we scale the coefficient Qi
// to Qi*2^i, and replace z by (x^2)/2.
// (B). To achieve maximum accuracy, we compute expm1(x) by
// (i) if x < -56*ln2, return -1.0, (raise inexact if x!=inf)
// (ii) if k=0, return r-E
// (iii) if k=-1, return 0.5*(r-E)-0.5
// (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E)
// else return 1.0+2.0*(r-E);
// (v) if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1)
// (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else
// (vii) return 2^k(1-((E+2^-k)-r))
//
// Special cases:
// expm1(INF) is INF, expm1(NaN) is NaN;
// expm1(-INF) is -1, and
// for finite argument, only expm1(0)=0 is exact.
//
// Accuracy:
// according to an error analysis, the error is always less than
// 1 ulp (unit in the last place).
//
// Misc. info.
// For IEEE double
// if x > 7.09782712893383973096e+02 then expm1(x) overflow
//
const KEXPM1_OVERFLOW = kMath[45];
const INVLN2 = kMath[46];
macro KEXPM1(x)
(kMath[47+x])
endmacro
function MathExpm1(x) {
x = x * 1; // Convert to number.
var y;
var hi;
var lo;
var k;
var t;
var c;
var hx = %_DoubleHi(x);
var xsb = hx & 0x80000000; // Sign bit of x
var y = (xsb === 0) ? x : -x; // y = |x|
hx &= 0x7fffffff; // High word of |x|
// Filter out huge and non-finite argument
if (hx >= 0x4043687a) { // if |x| ~=> 56 * ln2
if (hx >= 0x40862e42) { // if |x| >= 709.78
if (hx >= 0x7ff00000) {
// expm1(inf) = inf; expm1(-inf) = -1; expm1(nan) = nan;
return (x === -INFINITY) ? -1 : x;
}
if (x > KEXPM1_OVERFLOW) return INFINITY; // Overflow
}
if (xsb != 0) return -1; // x < -56 * ln2, return -1.
}
// Argument reduction
if (hx > 0x3fd62e42) { // if |x| > 0.5 * ln2
if (hx < 0x3ff0a2b2) { // and |x| < 1.5 * ln2
if (xsb === 0) {
hi = x - LN2_HI;
lo = LN2_LO;
k = 1;
} else {
hi = x + LN2_HI;
lo = -LN2_LO;
k = -1;
}
} else {
k = (INVLN2 * x + ((xsb === 0) ? 0.5 : -0.5)) | 0;
t = k;
// t * ln2_hi is exact here.
hi = x - t * LN2_HI;
lo = t * LN2_LO;
}
x = hi - lo;
c = (hi - x) - lo;
} else if (hx < 0x3c900000) {
// When |x| < 2^-54, we can return x.
return x;
} else {
// Fall through.
k = 0;
}
// x is now in primary range
var hfx = 0.5 * x;
var hxs = x * hfx;
var r1 = 1 + hxs * (KEXPM1(0) + hxs * (KEXPM1(1) + hxs *
(KEXPM1(2) + hxs * (KEXPM1(3) + hxs * KEXPM1(4)))));
t = 3 - r1 * hfx;
var e = hxs * ((r1 - t) / (6 - x * t));
if (k === 0) { // c is 0
return x - (x*e - hxs);
} else {
e = (x * (e - c) - c);
e -= hxs;
if (k === -1) return 0.5 * (x - e) - 0.5;
if (k === 1) {
if (x < -0.25) return -2 * (e - (x + 0.5));
return 1 + 2 * (x - e);
}
if (k <= -2 || k > 56) {
// suffice to return exp(x) + 1
y = 1 - (e - x);
// Add k to y's exponent
y = %_ConstructDouble(%_DoubleHi(y) + (k << 20), %_DoubleLo(y));
return y - 1;
}
if (k < 20) {
// t = 1 - 2^k
t = %_ConstructDouble(0x3ff00000 - (0x200000 >> k), 0);
y = t - (e - x);
// Add k to y's exponent
y = %_ConstructDouble(%_DoubleHi(y) + (k << 20), %_DoubleLo(y));
} else {
// t = 2^-k
t = %_ConstructDouble((0x3ff - k) << 20, 0);
y = x - (e + t);
y += 1;
// Add k to y's exponent
y = %_ConstructDouble(%_DoubleHi(y) + (k << 20), %_DoubleLo(y));
}
}
return y;
}
tools/generate-runtime-tests.py
View file @ 930e5ccc
......@@ -51,7 +51,7 @@ EXPECTED_FUNCTION_COUNT = 429
EXPECTED_FUZZABLE_COUNT = 330
EXPECTED_CCTEST_COUNT = 7
EXPECTED_UNKNOWN_COUNT = 17
EXPECTED_BUILTINS_COUNT = 809
EXPECTED_BUILTINS_COUNT = 808
# Don't call these at all.
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment