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() {
"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);
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 + ...
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));
-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
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
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
};
// Constants to be exposed to builtins via Float64Array.
struct MathConstants {
- static const double constants[45];
+ static const double constants[52];
};
}
} // namespace v8::internal
// 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
// 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];
// 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:
const LN2_LO = kMath[35];
const TWO54 = kMath[36];
const TWO_THIRD = kMath[37];
-macro KLOGP1(x)
+macro KLOG1P(x)
(kMath[38+x])
endmacro
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;
+}
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.
projects / platform / upstream / v8.git / commitdiff