Fix oflow bound in exp and pow.

Issue #600 highlighted a limit of current
implementation where early overflow is noticed,
which diverges from the C99 standard.

This patch adds test to catch this type of cases.

It only checks a few points therefore we don't
know for sure if overflow occurs at the right time
but we know if something is wrong in the overflow
region.

Fix both scalar and simd implementation (double precision)

Fixes #523.

src/common/misc.h

@@ -13,10 +13,15 @@
 #include <string.h>
 #endif
 
+
 #ifndef M_PI
 #define M_PI 3.141592653589793238462643383279502884
 #endif
 
+#ifndef M_PIf
+# define M_PIf ((float)M_PI)
+#endif
+
 #ifndef M_PIl
 #define M_PIl 3.141592653589793238462643383279502884L
 #endif
@@ -137,9 +142,10 @@
 #define L2Lf 1.428606765330187045e-06f
 
 #define R_LN2f 1.442695040888963407359924681001892137426645954152985934135449406931f
-#ifndef M_PIf
-# define M_PIf ((float)M_PI)
-#endif
+
+// Overflow bound for exp and pow
+
+#define LOG_DBL_MAX 0x1.62e42fefa39efp+9 /* 709.782712893384 */
 
 //
 

src/libm-tester/tester.c

@@ -3903,6 +3903,8 @@ void do_test() {
     fprintf(stderr, "exp : ");
     for(d = -10;d < 10 && success;d += 0.002) checkAccuracy_d(mpfr_exp, child_exp, d, 1.0);
     for(d = -1000;d < 1000 && success;d += 1.1) checkAccuracy_d(mpfr_exp, child_exp, d, 1.0);
+    // Test for early or late overflow, e.g before or after x = LOG_DBL_MAX
+    for(d = LOG_DBL_MAX - 0.0001;(d < LOG_DBL_MAX + 0.0001) && success;d += 0.00001) checkAccuracy_d(mpfr_exp, child_exp, d, 1.0);
     showResult(success);
 
     //
@@ -3914,6 +3916,8 @@ void do_test() {
       }
     }
     for(y = -1000;y < 1000 && success;y += 0.1) checkAccuracy_d_d(mpfr_pow, child_pow, 2.1, y, 1.0);
+    // Test for early or late overflow (test limited to x = e)
+    for(d = LOG_DBL_MAX - 0.0001;(d < LOG_DBL_MAX + 0.0001) && success;d += 0.00001) checkAccuracy_d_d(mpfr_pow, child_pow, exp(1.0), d, 1.0);
     showResult(success);
 
     //

src/libm/sleefdp.c

@@ -1490,9 +1490,9 @@ EXPORT CONST double xexp(double d) {
   u = s * s * u + s + 1;
   u = ldexp2k(u, q);
 
-  if (d > 709.78271114955742909217217426) u = SLEEF_INFINITY;
+  if (d > LOG_DBL_MAX) u = SLEEF_INFINITY;
   if (d < -1000) u = 0;
-  
+
   return u;
 }
 
@@ -1641,7 +1641,7 @@ EXPORT CONST double xpow(double x, double y) {
   Sleef_double2 d = ddmul_d2_d2_d(logk(fabsk(x)), y);
   double result = expk(d);
 
-  result = (d.x > 709.78271114955742909217217426 || xisnan(result)) ? SLEEF_INFINITY : result;
+  result = (d.x > LOG_DBL_MAX || xisnan(result)) ? SLEEF_INFINITY : result;
   result *= (x > 0 ? 1 : (yisint ? (yisodd ? -1 : 1) : SLEEF_NAN));
 
   double efx = mulsign(fabsk(x) - 1, y);

src/libm/sleefsimddp.c

@@ -2167,10 +2167,11 @@ EXPORT CONST VECTOR_CC vdouble xexp(vdouble d) {
 
   u = vadd_vd_vd_vd(vcast_vd_d(1), vmla_vd_vd_vd_vd(vmul_vd_vd_vd(s, s), u, s));
 #endif // #ifdef ENABLE_FMA_DP
-  
+
   u = vldexp2_vd_vd_vi(u, q);
 
-  u = vsel_vd_vo_vd_vd(vgt_vo_vd_vd(d, vcast_vd_d(709.78271114955742909217217426)), vcast_vd_d(SLEEF_INFINITY), u);
+  vopmask o = vgt_vo_vd_vd(d, vcast_vd_d(LOG_DBL_MAX));
+  u = vsel_vd_vo_vd_vd(o, vcast_vd_d(SLEEF_INFINITY), u);
   u = vreinterpret_vd_vm(vandnot_vm_vo64_vm(vlt_vo_vd_vd(d, vcast_vd_d(-1000)), vreinterpret_vm_vd(u)));
 
   return u;
@@ -2340,13 +2341,13 @@ static INLINE CONST VECTOR_CC vdouble expk(vdouble2 d) {
 
 #if !defined(DETERMINISTIC)
 EXPORT CONST VECTOR_CC vdouble xpow(vdouble x, vdouble y) {
-#if 1
   vopmask yisint = visint_vo_vd(y);
   vopmask yisodd = vand_vo_vo_vo(visodd_vo_vd(y), yisint);
 
   vdouble2 d = ddmul_vd2_vd2_vd(logk(vabs_vd_vd(x)), y);
   vdouble result = expk(d);
-  result = vsel_vd_vo_vd_vd(vgt_vo_vd_vd(vd2getx_vd_vd2(d), vcast_vd_d(709.78271114955742909217217426)), vcast_vd_d(SLEEF_INFINITY), result);
+  vopmask o = vgt_vo_vd_vd(vd2getx_vd_vd2(d), vcast_vd_d(LOG_DBL_MAX));
+  result = vsel_vd_vo_vd_vd(o, vcast_vd_d(SLEEF_INFINITY), result);
 
   result = vmul_vd_vd_vd(result,
 			 vsel_vd_vo_vd_vd(vgt_vo_vd_vd(x, vcast_vd_d(0)),
@@ -2372,9 +2373,6 @@ EXPORT CONST VECTOR_CC vdouble xpow(vdouble x, vdouble y) {
   result = vsel_vd_vo_vd_vd(vor_vo_vo_vo(veq_vo_vd_vd(y, vcast_vd_d(0)), veq_vo_vd_vd(x, vcast_vd_d(1))), vcast_vd_d(1), result);
 
   return result;
-#else
-  return expk(ddmul_vd2_vd2_vd(logk(x), y));
-#endif
 }
 #endif // #if !defined(DETERMINISTIC)