Index: Python/pymath.c =================================================================== --- Python/pymath.c (revision 67911) +++ Python/pymath.c (working copy) @@ -1,5 +1,16 @@ #include "Python.h" +/* On x86, call this function (preferably via the macro Py_FORCE_DOUBLE + defined in pymath.h) to force a floating-point number out of an x87 FPU + register and into memory, thus rounding from extended precision to double + precision. */ + +double _Py_force_double(double x) +{ + volatile double y = x; + return y; +} + #ifndef HAVE_HYPOT double hypot(double x, double y) { Index: configure =================================================================== --- configure (revision 67911) +++ configure (working copy) @@ -1,5 +1,5 @@ #! /bin/sh -# From configure.in Revision: 67227 . +# From configure.in Revision: 67463 . # Guess values for system-dependent variables and create Makefiles. # Generated by GNU Autoconf 2.61 for python 2.7. # @@ -13403,7 +13403,7 @@ fi # Dynamic linking for HP-UX -# only check for sem_ini if thread support is requested +# only check for sem_init if thread support is requested if test "$with_threads" = "yes" -o -z "$with_threads"; then { echo "$as_me:$LINENO: checking for library containing sem_init" >&5 echo $ECHO_N "checking for library containing sem_init... $ECHO_C" >&6; } @@ -21516,6 +21516,94 @@ LIBS_SAVE=$LIBS LIBS="$LIBS $LIBM" +# Detect whether system arithmetic is subject to x87-style double +# rounding issues. The result of this test has little meaning on non +# IEEE 754 platforms. (See http://bugs.python.org/issue2937.) +{ echo "$as_me:$LINENO: checking for x87-style double rounding" >&5 +echo $ECHO_N "checking for x87-style double rounding... $ECHO_C" >&6; } +if test "${ac_cv_x87_double_rounding+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + +if test "$cross_compiling" = yes; then + ac_cv_x87_double_rounding=no +else + cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ + +#include +#include +int main() { + /* Assume IEEE 754 double and round-to-nearest. Results in other + cases are undefined. */ + volatile double x, y, z; + /* 1./(1-2**-53) -> 1+2**-52 (correct), 1.0 (double rounding) */ + x = 0.99999999999999989; /* 1-2**-53 */ + y = 1./x; + if (y != 1.) + exit(0); + /* 1e16+2.99999 -> 1e16+2. (correct), 1e16+4. (double rounding) */ + x = 1e16; + y = 2.99999; + z = x + y; + if (z != 1e16+4.) + exit(0); + /* both tests show evidence of double rounding */ + exit(1); +} + +_ACEOF +rm -f conftest$ac_exeext +if { (ac_try="$ac_link" +case "(($ac_try" in + *\"* | *\`* | *\\*) ac_try_echo=\$ac_try;; + *) ac_try_echo=$ac_try;; +esac +eval "echo \"\$as_me:$LINENO: $ac_try_echo\"") >&5 + (eval "$ac_link") 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && { ac_try='./conftest$ac_exeext' + { (case "(($ac_try" in + *\"* | *\`* | *\\*) ac_try_echo=\$ac_try;; + *) ac_try_echo=$ac_try;; +esac +eval "echo \"\$as_me:$LINENO: $ac_try_echo\"") >&5 + (eval "$ac_try") 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; }; then + ac_cv_x87_double_rounding=no +else + echo "$as_me: program exited with status $ac_status" >&5 +echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + +( exit $ac_status ) +ac_cv_x87_double_rounding=yes +fi +rm -f core *.core core.conftest.* gmon.out bb.out conftest$ac_exeext conftest.$ac_objext conftest.$ac_ext +fi + + +fi + +{ echo "$as_me:$LINENO: result: $ac_cv_x87_double_rounding" >&5 +echo "${ECHO_T}$ac_cv_x87_double_rounding" >&6; } +if test "$ac_cv_x87_double_rounding" = yes +then + +cat >>confdefs.h <<\_ACEOF +#define X87_DOUBLE_ROUNDING 1 +_ACEOF + +fi + + # On FreeBSD 6.2, it appears that tanh(-0.) returns 0. instead of # -0. on some architectures. { echo "$as_me:$LINENO: checking whether tanh preserves the sign of zero" >&5 Index: Include/pymath.h =================================================================== --- Include/pymath.h (revision 67911) +++ Include/pymath.h (working copy) @@ -145,6 +145,23 @@ #define Py_NAN (Py_HUGE_VAL * 0.) #endif +/* On x86, Py_FORCE_DOUBLE forces a floating-point number out of an x87 FPU + register and into a 64-bit memory location, rounding from extended + precision to double precision in the process. It is particularly useful in + conjunction with Py_IS_INFINITY and Py_IS_FINITE, which otherwise can give + bogus results on x86 platforms that use the x87 FPU. */ + +/* take double rounding as evidence of x87 usage */ +#ifndef Py_FORCE_DOUBLE +#if defined(X87_DOUBLE_ROUNDING) +PyAPI_FUNC(double) _Py_force_double(double); +#define Py_FORCE_DOUBLE(x) (_Py_force_double(x)) +#else +#define Py_FORCE_DOUBLE(x) (x) +#endif +#endif + + /* Py_OVERFLOWED(X) * Return 1 iff a libm function overflowed. Set errno to 0 before calling * a libm function, and invoke this macro after, passing the function Index: configure.in =================================================================== --- configure.in (revision 67911) +++ configure.in (working copy) @@ -3136,6 +3136,44 @@ LIBS_SAVE=$LIBS LIBS="$LIBS $LIBM" +# Detect whether system arithmetic is subject to x87-style double +# rounding issues. The result of this test has little meaning on non +# IEEE 754 platforms. (See http://bugs.python.org/issue2937.) +AC_MSG_CHECKING(for x87-style double rounding) +AC_CACHE_VAL(ac_cv_x87_double_rounding, [ +AC_TRY_RUN([ +#include +#include +int main() { + /* Assume IEEE 754 double and round-to-nearest. Results in other + cases are undefined. */ + volatile double x, y, z; + /* 1./(1-2**-53) -> 1+2**-52 (correct), 1.0 (double rounding) */ + x = 0.99999999999999989; /* 1-2**-53 */ + y = 1./x; + if (y != 1.) + exit(0); + /* 1e16+2.99999 -> 1e16+2. (correct), 1e16+4. (double rounding) */ + x = 1e16; + y = 2.99999; + z = x + y; + if (z != 1e16+4.) + exit(0); + /* both tests show evidence of double rounding */ + exit(1); +} +], +ac_cv_x87_double_rounding=no, +ac_cv_x87_double_rounding=yes, +ac_cv_x87_double_rounding=no)]) +AC_MSG_RESULT($ac_cv_x87_double_rounding) +if test "$ac_cv_x87_double_rounding" = yes +then + AC_DEFINE(X87_DOUBLE_ROUNDING, 1, + [Define if arithmetic is subject to x87-style double rounding issue]) +fi + + # On FreeBSD 6.2, it appears that tanh(-0.) returns 0. instead of # -0. on some architectures. AC_MSG_CHECKING(whether tanh preserves the sign of zero) Index: Objects/complexobject.c =================================================================== --- Objects/complexobject.c (revision 67911) +++ Objects/complexobject.c (working copy) @@ -212,7 +212,7 @@ return Py_NAN; } result = hypot(z.real, z.imag); - if (!Py_IS_FINITE(result)) + if (!Py_IS_FINITE(Py_FORCE_DOUBLE(result))) errno = ERANGE; else errno = 0; Index: Modules/cmathmodule.c =================================================================== --- Modules/cmathmodule.c (revision 67911) +++ Modules/cmathmodule.c (working copy) @@ -418,7 +418,8 @@ r.imag = sin(z.imag) * sinh(z.real); } /* detect overflow, and set errno accordingly */ - if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag)) + if (Py_IS_INFINITY(Py_FORCE_DOUBLE(r.real)) || + Py_IS_INFINITY(Py_FORCE_DOUBLE(r.imag))) errno = ERANGE; else errno = 0; @@ -478,7 +479,8 @@ r.imag = l*sin(z.imag); } /* detect overflow, and set errno accordingly */ - if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag)) + if (Py_IS_INFINITY(Py_FORCE_DOUBLE(r.real)) || + Py_IS_INFINITY(Py_FORCE_DOUBLE(r.imag))) errno = ERANGE; else errno = 0; @@ -647,7 +649,8 @@ r.imag = sin(z.imag) * cosh(z.real); } /* detect overflow, and set errno accordingly */ - if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag)) + if (Py_IS_INFINITY(Py_FORCE_DOUBLE(r.real)) || + Py_IS_INFINITY(Py_FORCE_DOUBLE(r.imag))) errno = ERANGE; else errno = 0; Index: Modules/mathmodule.c =================================================================== --- Modules/mathmodule.c (revision 67911) +++ Modules/mathmodule.c (working copy) @@ -729,7 +729,7 @@ PyFPE_START_PROTECT("in math_ldexp", return 0); r = ldexp(x, (int)exp); PyFPE_END_PROTECT(r); - if (Py_IS_INFINITY(r)) + if (Py_IS_INFINITY(Py_FORCE_DOUBLE(r))) errno = ERANGE; } @@ -754,7 +754,7 @@ return Py_BuildValue("(dd)", copysign(0., x), x); else if (Py_IS_NAN(x)) return Py_BuildValue("(dd)", x, x); - } + } errno = 0; PyFPE_START_PROTECT("in math_modf", return 0); @@ -902,7 +902,7 @@ else errno = 0; } - else if (Py_IS_INFINITY(r)) { + else if (Py_IS_INFINITY(Py_FORCE_DOUBLE(r))) { if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) errno = ERANGE; else @@ -986,7 +986,7 @@ (A) (+/-0.)**negative (-> divide-by-zero) (B) overflow of x**y with x and y finite */ - else if (Py_IS_INFINITY(r)) { + else if (Py_IS_INFINITY(Py_FORCE_DOUBLE(r))) { if (x == 0.) errno = EDOM; else Index: pyconfig.h.in =================================================================== --- pyconfig.h.in (revision 67911) +++ pyconfig.h.in (working copy) @@ -983,6 +983,9 @@ first (like Motorola and SPARC, unlike Intel and VAX). */ #undef WORDS_BIGENDIAN +/* Define if arithmetic is subject to x87-style double rounding issue */ +#undef X87_DOUBLE_ROUNDING + /* Define to 1 if on AIX 3. System headers sometimes define this. We just want to avoid a redefinition error message. */