[builtins] Make sure the Math functions and constants agree.

src/base/ieee754.cc

 // The following is adapted from fdlibm (http://www.netlib.org/fdlibm).
 //
 // ====================================================
 // Copyright (C) 1993 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 1232 matching lines @@
       u_threshold = -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */
       ln2HI[2] = {6.93147180369123816490e-01,    /* 0x3fe62e42, 0xfee00000 */
                   -6.93147180369123816490e-01},  /* 0xbfe62e42, 0xfee00000 */
       ln2LO[2] = {1.90821492927058770002e-10,    /* 0x3dea39ef, 0x35793c76 */
                   -1.90821492927058770002e-10},  /* 0xbdea39ef, 0x35793c76 */
       invln2 = 1.44269504088896338700e+00,       /* 0x3ff71547, 0x652b82fe */
       P1 = 1.66666666666666019037e-01,           /* 0x3FC55555, 0x5555553E */
       P2 = -2.77777777770155933842e-03,          /* 0xBF66C16C, 0x16BEBD93 */
       P3 = 6.61375632143793436117e-05,           /* 0x3F11566A, 0xAF25DE2C */
       P4 = -1.65339022054652515390e-06,          /* 0xBEBBBD41, 0xC5D26BF1 */
-      P5 = 4.13813679705723846039e-08;           /* 0x3E663769, 0x72BEA4D0 */
+      P5 = 4.13813679705723846039e-08,           /* 0x3E663769, 0x72BEA4D0 */
+      E = 2.718281828459045;                     /* 0x4005bf0a, 0x8b145769 */
 
   static volatile double
       huge = 1.0e+300,
       twom1000 = 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0*/
       two1023 = 8.988465674311579539e307;     /* 0x1p1023 */
 
   double y, hi = 0.0, lo = 0.0, c, t, twopk;
   int32_t k = 0, xsb;
   uint32_t hx;
 
@@ skipping 11 matching lines @@
       else
         return (xsb == 0) ? x : 0.0; /* exp(+-inf)={inf,0} */
     }
     if (x > o_threshold) return huge * huge;         /* overflow */
     if (x < u_threshold) return twom1000 * twom1000; /* underflow */
   }
 
   /* argument reduction */
   if (hx > 0x3fd62e42) {   /* if  |x| > 0.5 ln2 */
     if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
+      /* TODO(rtoy): We special case exp(1) here to return the correct
+       * value of E, as the computation below would get the last bit
+       * wrong. We should probably fix the algorithm instead.
+       */
+      if (x == 1.0) return E;

[Raymond Toy, 2016/06/21 16:20:12]

It's probably not possible to fix the algorithm easily to make this round
correctly.  fdlibm does not produce correctly rounded results; it's only
accurate to about 1 ulp.

I think you have implemented some of my tests to ensure that exp(x) is still
monotonic near x = 1.  This is an important property that should be maintained.

[Benedikt Meurer, 2016/06/22 03:55:59]

Since exp is not performance critical or anything, would it make sense to import
an algorithm that produces correctly rounded results? Is there one w/ managable
complexity?

[Raymond Toy, 2016/06/22 22:11:48]

Yes, this is possible, and something that I would really like to see.  But I
think it's a fairly large undertaking and should probably done all at "once". 
It turns out that some of the correctly rounded algorithms for some functions
are actually faster than fdlibm.

This has been a long-term goal for me and I really need to get back on it!

[Benedikt Meurer, 2016/06/23 03:30:58]

This would be super awesome and highly appreciated if you could take over the
improvements to the Math functions (again). :-)
We'll hopefully finish the move of the remaining Math functions to this file,
and then it should be easier to experiment and change, since everything has only
a single C++ implementation.

       hi = x - ln2HI[xsb];
       lo = ln2LO[xsb];
       k = 1 - xsb - xsb;
     } else {
       k = static_cast<int>(invln2 * x + halF[xsb]);
       t = k;
       hi = x - t * ln2HI[0]; /* t*ln2HI is exact here */
       lo = t * ln2LO[0];
     }
     STRICT_ASSIGN(double, x, hi - lo);
@@ skipping 1003 matching lines @@
     /* argument reduction needed */
     n = __ieee754_rem_pio2(x, y);
     /* 1 -> n even, -1 -> n odd */
     return __kernel_tan(y[0], y[1], 1 - ((n & 1) << 1));
   }
 }
 
 }  // namespace ieee754
 }  // namespace base
 }  // namespace v8