diff --git a/src/common/misc.h b/src/common/misc.h index 472cae68..258d39c6 100644 --- a/src/common/misc.h +++ b/src/common/misc.h @@ -13,10 +13,15 @@ #include #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 */ // diff --git a/src/libm-tester/tester.c b/src/libm-tester/tester.c index 5e2cf2ac..2831056f 100644 --- a/src/libm-tester/tester.c +++ b/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); // diff --git a/src/libm/sleefdp.c b/src/libm/sleefdp.c index e5ffd857..09a40de1 100644 --- a/src/libm/sleefdp.c +++ b/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); diff --git a/src/libm/sleefsimddp.c b/src/libm/sleefsimddp.c index 1600969a..4dfc434e 100644 --- a/src/libm/sleefsimddp.c +++ b/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)