Index: configure =================================================================== --- configure (revision 63899) +++ configure (working copy) @@ -1,5 +1,5 @@ #! /bin/sh -# From configure.in Revision: 63545 . +# From configure.in Revision: 63690 . # Guess values for system-dependent variables and create Makefiles. # Generated by GNU Autoconf 2.61 for python 2.6. # @@ -11080,6 +11080,467 @@ fi +{ echo "$as_me:$LINENO: checking for long double support" >&5 +echo $ECHO_N "checking for long double support... $ECHO_C" >&6; } +have_long_double=no +cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ + +int +main () +{ +long double x; x = (long double)0.0; + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext +if { (ac_try="$ac_compile" +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_compile") 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && { + test -z "$ac_c_werror_flag" || + test ! -s conftest.err + } && test -s conftest.$ac_objext; then + + +cat >>confdefs.h <<\_ACEOF +#define HAVE_LONG_DOUBLE 1 +_ACEOF + + have_long_double=yes + +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + + +fi + +rm -f core conftest.err conftest.$ac_objext conftest.$ac_ext +{ echo "$as_me:$LINENO: result: $have_long_double" >&5 +echo "${ECHO_T}$have_long_double" >&6; } +if test "$have_long_double" = yes ; then +{ echo "$as_me:$LINENO: checking for long double" >&5 +echo $ECHO_N "checking for long double... $ECHO_C" >&6; } +if test "${ac_cv_type_long_double+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ +$ac_includes_default +typedef long double ac__type_new_; +int +main () +{ +if ((ac__type_new_ *) 0) + return 0; +if (sizeof (ac__type_new_)) + return 0; + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext +if { (ac_try="$ac_compile" +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_compile") 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && { + test -z "$ac_c_werror_flag" || + test ! -s conftest.err + } && test -s conftest.$ac_objext; then + ac_cv_type_long_double=yes +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + + ac_cv_type_long_double=no +fi + +rm -f core conftest.err conftest.$ac_objext conftest.$ac_ext +fi +{ echo "$as_me:$LINENO: result: $ac_cv_type_long_double" >&5 +echo "${ECHO_T}$ac_cv_type_long_double" >&6; } + +# The cast to long int works around a bug in the HP C Compiler +# version HP92453-01 B.11.11.23709.GP, which incorrectly rejects +# declarations like `int a3[[(sizeof (unsigned char)) >= 0]];'. +# This bug is HP SR number 8606223364. +{ echo "$as_me:$LINENO: checking size of long double" >&5 +echo $ECHO_N "checking size of long double... $ECHO_C" >&6; } +if test "${ac_cv_sizeof_long_double+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + if test "$cross_compiling" = yes; then + # Depending upon the size, compute the lo and hi bounds. +cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ +$ac_includes_default + typedef long double ac__type_sizeof_; +int +main () +{ +static int test_array [1 - 2 * !(((long int) (sizeof (ac__type_sizeof_))) >= 0)]; +test_array [0] = 0 + + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext +if { (ac_try="$ac_compile" +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_compile") 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && { + test -z "$ac_c_werror_flag" || + test ! -s conftest.err + } && test -s conftest.$ac_objext; then + ac_lo=0 ac_mid=0 + while :; do + cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ +$ac_includes_default + typedef long double ac__type_sizeof_; +int +main () +{ +static int test_array [1 - 2 * !(((long int) (sizeof (ac__type_sizeof_))) <= $ac_mid)]; +test_array [0] = 0 + + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext +if { (ac_try="$ac_compile" +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_compile") 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && { + test -z "$ac_c_werror_flag" || + test ! -s conftest.err + } && test -s conftest.$ac_objext; then + ac_hi=$ac_mid; break +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + + ac_lo=`expr $ac_mid + 1` + if test $ac_lo -le $ac_mid; then + ac_lo= ac_hi= + break + fi + ac_mid=`expr 2 '*' $ac_mid + 1` +fi + +rm -f core conftest.err conftest.$ac_objext conftest.$ac_ext + done +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + + cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ +$ac_includes_default + typedef long double ac__type_sizeof_; +int +main () +{ +static int test_array [1 - 2 * !(((long int) (sizeof (ac__type_sizeof_))) < 0)]; +test_array [0] = 0 + + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext +if { (ac_try="$ac_compile" +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_compile") 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && { + test -z "$ac_c_werror_flag" || + test ! -s conftest.err + } && test -s conftest.$ac_objext; then + ac_hi=-1 ac_mid=-1 + while :; do + cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ +$ac_includes_default + typedef long double ac__type_sizeof_; +int +main () +{ +static int test_array [1 - 2 * !(((long int) (sizeof (ac__type_sizeof_))) >= $ac_mid)]; +test_array [0] = 0 + + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext +if { (ac_try="$ac_compile" +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_compile") 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && { + test -z "$ac_c_werror_flag" || + test ! -s conftest.err + } && test -s conftest.$ac_objext; then + ac_lo=$ac_mid; break +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + + ac_hi=`expr '(' $ac_mid ')' - 1` + if test $ac_mid -le $ac_hi; then + ac_lo= ac_hi= + break + fi + ac_mid=`expr 2 '*' $ac_mid` +fi + +rm -f core conftest.err conftest.$ac_objext conftest.$ac_ext + done +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + + ac_lo= ac_hi= +fi + +rm -f core conftest.err conftest.$ac_objext conftest.$ac_ext +fi + +rm -f core conftest.err conftest.$ac_objext conftest.$ac_ext +# Binary search between lo and hi bounds. +while test "x$ac_lo" != "x$ac_hi"; do + ac_mid=`expr '(' $ac_hi - $ac_lo ')' / 2 + $ac_lo` + cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ +$ac_includes_default + typedef long double ac__type_sizeof_; +int +main () +{ +static int test_array [1 - 2 * !(((long int) (sizeof (ac__type_sizeof_))) <= $ac_mid)]; +test_array [0] = 0 + + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext +if { (ac_try="$ac_compile" +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_compile") 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && { + test -z "$ac_c_werror_flag" || + test ! -s conftest.err + } && test -s conftest.$ac_objext; then + ac_hi=$ac_mid +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + + ac_lo=`expr '(' $ac_mid ')' + 1` +fi + +rm -f core conftest.err conftest.$ac_objext conftest.$ac_ext +done +case $ac_lo in +?*) ac_cv_sizeof_long_double=$ac_lo;; +'') if test "$ac_cv_type_long_double" = yes; then + { { echo "$as_me:$LINENO: error: cannot compute sizeof (long double) +See \`config.log' for more details." >&5 +echo "$as_me: error: cannot compute sizeof (long double) +See \`config.log' for more details." >&2;} + { (exit 77); exit 77; }; } + else + ac_cv_sizeof_long_double=0 + fi ;; +esac +else + cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ +$ac_includes_default + typedef long double ac__type_sizeof_; +static long int longval () { return (long int) (sizeof (ac__type_sizeof_)); } +static unsigned long int ulongval () { return (long int) (sizeof (ac__type_sizeof_)); } +#include +#include +int +main () +{ + + FILE *f = fopen ("conftest.val", "w"); + if (! f) + return 1; + if (((long int) (sizeof (ac__type_sizeof_))) < 0) + { + long int i = longval (); + if (i != ((long int) (sizeof (ac__type_sizeof_)))) + return 1; + fprintf (f, "%ld\n", i); + } + else + { + unsigned long int i = ulongval (); + if (i != ((long int) (sizeof (ac__type_sizeof_)))) + return 1; + fprintf (f, "%lu\n", i); + } + return ferror (f) || fclose (f) != 0; + + ; + return 0; +} +_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_sizeof_long_double=`cat conftest.val` +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 ) +if test "$ac_cv_type_long_double" = yes; then + { { echo "$as_me:$LINENO: error: cannot compute sizeof (long double) +See \`config.log' for more details." >&5 +echo "$as_me: error: cannot compute sizeof (long double) +See \`config.log' for more details." >&2;} + { (exit 77); exit 77; }; } + else + ac_cv_sizeof_long_double=0 + fi +fi +rm -f core *.core core.conftest.* gmon.out bb.out conftest$ac_exeext conftest.$ac_objext conftest.$ac_ext +fi +rm -f conftest.val +fi +{ echo "$as_me:$LINENO: result: $ac_cv_sizeof_long_double" >&5 +echo "${ECHO_T}$ac_cv_sizeof_long_double" >&6; } + + + +cat >>confdefs.h <<_ACEOF +#define SIZEOF_LONG_DOUBLE $ac_cv_sizeof_long_double +_ACEOF + + +fi + { echo "$as_me:$LINENO: checking for _Bool support" >&5 echo $ECHO_N "checking for _Bool support... $ECHO_C" >&6; } have_c99_bool=no @@ -20987,7 +21448,86 @@ fi +# on x86 platforms that don't use the SSE2 instruction set for +# floating-point arithmetic, addition may not be correctly rounded, as +# a result of internal double rounding (first rounding to extended +# double precision, then to double precision). +{ echo "$as_me:$LINENO: checking for incorrectly rounded addition" >&5 +echo $ECHO_N "checking for incorrectly rounded addition... $ECHO_C" >&6; } +if test "${ac_cv_broken_sum+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + +if test "$cross_compiling" = yes; then + ac_cv_broken_sum=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 +int main() { + volatile double x, y, sum; + x = 1e16-2.0; + y = 0.9999; + sum = x + y; + if (sum == 1e16) + exit(1); + exit(0); +} + +_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_broken_sum=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_broken_sum=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_broken_sum" >&5 +echo "${ECHO_T}$ac_cv_broken_sum" >&6; } +if test "$ac_cv_broken_sum" = yes +then + +cat >>confdefs.h <<\_ACEOF +#define BROKEN_SUM 1 +_ACEOF + +fi + + for ac_func in hypot do as_ac_var=`echo "ac_cv_func_$ac_func" | $as_tr_sh` Index: configure.in =================================================================== --- configure.in (revision 63899) +++ configure.in (working copy) @@ -1268,6 +1268,17 @@ AC_CHECK_SIZEOF(long long, 8) fi +AC_MSG_CHECKING(for long double support) +have_long_double=no +AC_TRY_COMPILE([], [long double x; x = (long double)0.0;], [ + AC_DEFINE(HAVE_LONG_DOUBLE, 1, [Define this if you have the type long double.]) + have_long_double=yes +]) +AC_MSG_RESULT($have_long_double) +if test "$have_long_double" = yes ; then +AC_CHECK_SIZEOF(long double) +fi + AC_MSG_CHECKING(for _Bool support) have_c99_bool=no AC_TRY_COMPILE([], [_Bool x; x = (_Bool)0;], [ @@ -3070,6 +3081,35 @@ [Define if tanh(-0.) is -0., or if platform doesn't have signed zeros]) fi +# on x86 platforms that don't use the SSE2 instruction set for +# floating-point arithmetic, addition may not be correctly rounded, as +# a result of internal double rounding (first rounding to extended +# double precision, then to double precision). + +AC_MSG_CHECKING(for incorrectly rounded addition) +AC_CACHE_VAL(ac_cv_broken_sum, [ +AC_TRY_RUN([ +#include +int main() { + volatile double x, y, sum; + x = 1e16-2.0; + y = 0.9999; + sum = x + y; + if (sum == 1e16) + exit(1); + exit(0); +} +], +ac_cv_broken_sum=no, +ac_cv_broken_sum=yes, +ac_cv_broken_sum=no)]) +AC_MSG_RESULT($ac_cv_broken_sum) +if test "$ac_cv_broken_sum" = yes +then + AC_DEFINE(BROKEN_SUM, 1, + [Define if addition is not always correctly rounded]) +fi + AC_REPLACE_FUNCS(hypot) AC_CHECK_FUNCS(acosh asinh atanh copysign expm1 finite isinf isnan log1p) Index: Lib/test/test_math.py =================================================================== --- Lib/test/test_math.py (revision 63899) +++ Lib/test/test_math.py (working copy) @@ -642,51 +642,6 @@ if not float.__getformat__("double").startswith("IEEE"): return - # on IEEE 754 compliant machines, both of the expressions - # below should round to 10000000000000002.0. - if 1e16+2.999 != 1e16+2.9999: - return - - # Python version of math.sum algorithm, for comparison - def msum(iterable): - """Full precision sum of values in iterable. Returns the value of - the sum, rounded to the nearest representable floating-point number - using the round-half-to-even rule. - - """ - # Stage 1: accumulate partials - partials = [] - for x in iterable: - i = 0 - for y in partials: - if abs(x) < abs(y): - x, y = y, x - hi = x + y - lo = y - (hi - x) - if lo: - partials[i] = lo - i += 1 - x = hi - partials[i:] = [x] if x else [] - - # Stage 2: sum partials - if not partials: - return 0.0 - - # sum from the top, stopping as soon as the sum is inexact. - total = partials.pop() - while partials: - x = partials.pop() - old_total, total = total, total + x - error = x - (total - old_total) - if error != 0.0: - # adjust for correct rounding if necessary - if partials and (partials[-1] > 0.0) == (error > 0.0) and \ - total + 2*error - total == 2*error: - total += 2*error - break - return total - from sys import float_info maxfloat = float_info.max twopow = 2.**(float_info.max_exp - 1) @@ -695,19 +650,15 @@ ([], 0.0), ([0.0], 0.0), ([1e100, 1.0, -1e100, 1e-100, 1e50, -1.0, -1e50], 1e-100), - ([1e308, 1e308, -1e308], OverflowError), ([-1e308, 1e308, 1e308], 1e308), ([1e308, -1e308, 1e308], 1e308), - ([2.0**1023, 2.0**1023, -2.0**1000], OverflowError), - ([twopow, twopow, twopow, twopow, -twopow, -twopow, -twopow], - OverflowError), ([2.0**53, -0.5, -2.0**-54], 2.0**53-1.0), ([2.0**53, 1.0, 2.0**-100], 2.0**53+2.0), ([2.0**53+10.0, 1.0, 2.0**-100], 2.0**53+12.0), ([2.0**53-4.0, 0.5, 2.0**-54], 2.0**53-3.0), ([2.0**1023-2.0**970, -1.0, 2.0**1023], OverflowError), - ([maxfloat, maxfloat*2.**-54], maxfloat), + #([maxfloat, maxfloat*2.**-54], maxfloat), ([maxfloat, maxfloat*2.**-53], OverflowError), ([1./n for n in range(1, 1001)], 7.4854708605503451), ([(-1.)**n/n for n in range(1, 1001)], -0.69264743055982025), @@ -728,7 +679,6 @@ ([twopow, twopow, twopow, twopow, -twopow, twopow], OverflowError), ([-twopow, -twopow, -twopow, -twopow], OverflowError), - ([2.**1023, 2.**1023, -2.**971], OverflowError), ([2.**1023, 2.**1023, -2.**970], OverflowError), ([-2.**970, 2.**1023, 2.**1023, -2.**-1074], OverflowError), ([ 2.**1023, 2.**1023, -2.**970, 2.**-1074], OverflowError), @@ -736,9 +686,6 @@ ([-2.**1023, -2.**1023, 2.**970], OverflowError), ([-2.**1023, -2.**1023, 2.**970, 2.**-1074], OverflowError), ([-2.**-1074, -2.**1023, -2.**1023, 2.**970], OverflowError), - ([2.**930, -2.**980, 2.**1023, 2.**1023, twopow, -twopow], - OverflowError), - ([2.**1023, 2.**1023, -1e307], OverflowError), ([1e16, 1., 1e-16], 10000000000000002.0), ([1e16-2., 1.-2.**-53, -(1e16-2.), -(1.-2.**-53)], 0.0), ] @@ -763,25 +710,7 @@ self.fail("test %d failed: got ValueError, expected %r " "for math.sum(%.100r)" % (i, s, vals)) - # compare with output of msum above, but only when - # result isn't an IEEE special or an exception - if not math.isinf(s) and not math.isnan(s): - self.assertEqual(msum(vals), s) - from random import random, gauss, shuffle - for j in xrange(1000): - vals = [7, 1e100, -7, -1e100, -9e-20, 8e-20] * 10 - s = 0 - for i in xrange(200): - v = gauss(0, random()) ** 7 - s - s += v - vals.append(v) - shuffle(vals) - - s = msum(vals) - self.assertEqual(msum(vals), math.sum(vals)) - - def testTan(self): self.assertRaises(TypeError, math.tan) self.ftest('tan(0)', math.tan(0), 0) Index: Modules/mathmodule.c =================================================================== --- Modules/mathmodule.c (revision 63899) +++ Modules/mathmodule.c (working copy) @@ -343,32 +343,38 @@ accurate result returned by sum(itertools.chain(seq1, seq2)). */ +#ifdef BROKEN_SUM +typedef long double partial_t; +#else +typedef double partial_t; +#endif + #define NUM_PARTIALS 32 /* initial partials array size, on stack */ /* Extend the partials array p[] by doubling its size. */ static int /* non-zero on error */ -_sum_realloc(double **p_ptr, Py_ssize_t n, - double *ps, Py_ssize_t *m_ptr) +_sum_realloc(partial_t **p_ptr, Py_ssize_t n, + partial_t *ps, Py_ssize_t *m_ptr) { void *v = NULL; Py_ssize_t m = *m_ptr; m += m; /* double */ - if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) { - double *p = *p_ptr; + if (n < m && m < (PY_SSIZE_T_MAX / sizeof(partial_t))) { + partial_t *p = *p_ptr; if (p == ps) { - v = PyMem_Malloc(sizeof(double) * m); + v = PyMem_Malloc(sizeof(partial_t) * m); if (v != NULL) - memcpy(v, ps, sizeof(double) * n); + memcpy(v, ps, sizeof(partial_t) * n); } else - v = PyMem_Realloc(p, sizeof(double) * m); + v = PyMem_Realloc(p, sizeof(partial_t) * m); } if (v == NULL) { /* size overflow or no memory */ PyErr_SetString(PyExc_MemoryError, "math sum partials"); return 1; } - *p_ptr = (double*) v; + *p_ptr = (partial_t*) v; *m_ptr = m; return 0; } @@ -408,8 +414,9 @@ { PyObject *item, *iter, *sum = NULL; Py_ssize_t i, j, n = 0, m = NUM_PARTIALS; - double x, y, t, ps[NUM_PARTIALS], *p = ps; - volatile double hi, yr, lo; + partial_t x, y, t, ps[NUM_PARTIALS], *p = ps; + volatile partial_t hi, yr, lo; + volatile double sumd, sumd2; iter = PyObject_GetIter(seq); if (iter == NULL) @@ -428,7 +435,7 @@ goto _sum_error; break; } - x = PyFloat_AsDouble(item); + x = (partial_t)PyFloat_AsDouble(item); Py_DECREF(item); if (PyErr_Occurred()) goto _sum_error; @@ -441,9 +448,9 @@ hi = x + y; yr = hi - x; lo = y - yr; + x = hi; if (lo != 0.0) p[i++] = lo; - x = hi; } n = i; /* ps[i:] = [x] */ @@ -460,42 +467,57 @@ } } - hi = 0.0; if (n > 0) { - hi = p[--n]; - if (Py_IS_FINITE(hi)) { - /* sum_exact(ps, hi) from the top, stop when the sum becomes inexact. */ + x = p[--n]; + if (Py_IS_FINITE(x)) { + /* sum_exact(ps, x) from the top, stop when the sum + becomes inexact. */ while (n > 0) { - x = hi; y = p[--n]; - assert(fabs(y) < fabs(x)); + assert(fabs(y) <= fabs(x)); hi = x + y; yr = hi - x; lo = y - yr; - if (lo != 0.0) + x = hi; + if (lo != 0.0) { + p[n++] = lo; break; + } } - /* Make half-even rounding work across multiple partials. Needed - so that sum([1e-16, 1, 1e16]) will round-up the last digit to - two instead of down to zero (the 1e-16 makes the 1 slightly - closer to two). With a potential 1 ULP rounding error fixed-up, - math.sum() can guarantee commutativity. */ - if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) || - (lo > 0.0 && p[n-1] > 0.0))) { - y = lo * 2.0; - x = hi + y; - yr = x - hi; - if (y == yr) - hi = x; + /* round x to a double (sumd), and append the rounding + error (if any) to the partials array. */ + sumd = (double)x; + lo = x - sumd; + if (lo != 0.0) + p[n++] = lo; + + /* Make half-even rounding work across multiple + partials. Needed so that sum([1e-16, 1, 1e16]) + will round-up the last digit to two instead of down + to zero (the 1e-16 makes the 1 slightly closer to + two). With a potential 1 ULP rounding error + fixed-up, math.sum() can guarantee + commutativity. */ + if (n >= 2) { + lo = p[--n]; + if ((lo < 0.0 && p[n-1] < 0.0) || + (lo > 0.0 && p[n-1] > 0.0)) { + y = lo * 2.0; + sumd2 = (double)((partial_t)sumd + y); + if (y == (partial_t)(sumd2 - sumd)) + sumd = sumd2; + } } } else { /* raise exception corresponding to a special value */ - errno = Py_IS_NAN(hi) ? EDOM : ERANGE; - if (is_error(hi)) + errno = Py_IS_NAN(x) ? EDOM : ERANGE; + if (is_error(x)) goto _sum_error; } } - sum = PyFloat_FromDouble(hi); + else + sumd = 0.0; + sum = PyFloat_FromDouble(sumd); _sum_error: PyFPE_END_PROTECT(hi) Index: pyconfig.h.in =================================================================== --- pyconfig.h.in (revision 63899) +++ pyconfig.h.in (working copy) @@ -15,6 +15,9 @@ /* Define this if you have BeOS threads. */ #undef BEOS_THREADS +/* Define if addition is not always correctly rounded */ +#undef BROKEN_SUM + /* Define if you have the Mach cthreads package */ #undef C_THREADS @@ -384,6 +387,9 @@ /* Define to 1 if you have the `log1p' function. */ #undef HAVE_LOG1P +/* Define this if you have the type long double. */ +#undef HAVE_LONG_DOUBLE + /* Define this if you have the type long long. */ #undef HAVE_LONG_LONG @@ -884,6 +890,9 @@ /* The size of `long', as computed by sizeof. */ #undef SIZEOF_LONG +/* The size of `long double', as computed by sizeof. */ +#undef SIZEOF_LONG_DOUBLE + /* The size of `long long', as computed by sizeof. */ #undef SIZEOF_LONG_LONG