Merged revisions 65258,65292,65299,65308-65309,65315,65326 via svnmerge from
svn+ssh://pythondev@svn.python.org/python/trunk ........ r65258 | mark.dickinson | 2008-07-27 08:15:29 +0100 (Sun, 27 Jul 2008) | 4 lines Remove math.sum tests related to overflow, special values, and behaviour near the extremes of the floating-point range. (The behaviour of math.sum should be regarded as undefined in these cases.) ........ r65292 | mark.dickinson | 2008-07-29 19:45:38 +0100 (Tue, 29 Jul 2008) | 4 lines More modifications to tests for math.sum: replace the Python version of msum by a version using a different algorithm, and use the new float.fromhex method to specify test results exactly. ........ r65299 | mark.dickinson | 2008-07-30 13:01:41 +0100 (Wed, 30 Jul 2008) | 5 lines Fix special-value handling for math.sum. Also minor cleanups to the code: fix tabbing, remove trailing whitespace, and reformat to fit into 80 columns. ........ r65308 | mark.dickinson | 2008-07-30 17:20:10 +0100 (Wed, 30 Jul 2008) | 2 lines Rename math.sum to math.fsum ........ r65309 | mark.dickinson | 2008-07-30 17:25:16 +0100 (Wed, 30 Jul 2008) | 3 lines Replace math.sum with math.fsum in a couple of comments that were missed by r65308 ........ r65315 | mark.dickinson | 2008-07-30 21:23:15 +0100 (Wed, 30 Jul 2008) | 2 lines Add note about problems with math.fsum on x86 hardware. ........ r65326 | mark.dickinson | 2008-07-31 15:48:32 +0100 (Thu, 31 Jul 2008) | 2 lines Rename testSum to testFsum and move it to proper place in test_math.py ........
This commit is contained in:
Mark Dickinson
parent
c57df32319
commit
aa7633ab94
|
|
@ -76,6 +76,42 @@ Number-theoretic and representation functions:
|
||||||
apart" the internal representation of a float in a portable way.
|
apart" the internal representation of a float in a portable way.
|
||||||
|
|
||||||
|
|
||||||
|
.. function:: fsum(iterable)
|
||||||
|
|
||||||
|
Return an accurate floating point sum of values in the iterable. Avoids
|
||||||
|
loss of precision by tracking multiple intermediate partial sums. The
|
||||||
|
algorithm's accuracy depends on IEEE-754 arithmetic guarantees and the
|
||||||
|
typical case where the rounding mode is half-even.
|
||||||
|
|
||||||
|
.. note::
|
||||||
|
|
||||||
|
On platforms where arithmetic results are not correctly rounded,
|
||||||
|
:func:`fsum` may occasionally produce incorrect results; these
|
||||||
|
results should be no less accurate than those from the builtin
|
||||||
|
:func:`sum` function, but nevertheless may have arbitrarily
|
||||||
|
large relative error.
|
||||||
|
|
||||||
|
In particular, this affects some older Intel hardware (for
|
||||||
|
example Pentium and earlier x86 processors) that makes use of
|
||||||
|
'extended precision' floating-point registers with 64 bits of
|
||||||
|
precision instead of the 53 bits of precision provided by a C
|
||||||
|
double. Arithmetic operations using these registers may be
|
||||||
|
doubly rounded (rounded first to 64 bits, and then rerounded to
|
||||||
|
53 bits), leading to incorrectly rounded results. To test
|
||||||
|
whether your machine is one of those affected, try the following
|
||||||
|
at a Python prompt::
|
||||||
|
|
||||||
|
>>> 1e16 + 2.9999
|
||||||
|
10000000000000002.0
|
||||||
|
|
||||||
|
Machines subject to the double-rounding problem described above
|
||||||
|
are likely to print ``10000000000000004.0`` instead of
|
||||||
|
``10000000000000002.0``.
|
||||||
|
|
||||||
|
|
||||||
|
.. versionadded:: 2.6
|
||||||
|
|
||||||
|
|
||||||
.. function:: isinf(x)
|
.. function:: isinf(x)
|
||||||
|
|
||||||
Checks if the float *x* is positive or negative infinite.
|
Checks if the float *x* is positive or negative infinite.
|
||||||
|
|
@ -100,12 +136,6 @@ Number-theoretic and representation functions:
|
||||||
Return the fractional and integer parts of *x*. Both results carry the sign of
|
Return the fractional and integer parts of *x*. Both results carry the sign of
|
||||||
*x*, and both are floats.
|
*x*, and both are floats.
|
||||||
|
|
||||||
.. function:: sum(iterable)
|
|
||||||
|
|
||||||
Return an accurate floating point sum of values in the iterable. Avoids
|
|
||||||
loss of precision by tracking multiple intermediate partial sums. The
|
|
||||||
algorithm's accuracy depends on IEEE-754 arithmetic guarantees and the
|
|
||||||
typical case where the rounding mode is half-even.
|
|
||||||
|
|
||||||
.. function:: trunc(x)
|
.. function:: trunc(x)
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -1537,7 +1537,7 @@ Here are all of the changes that Python 2.6 makes to the core Python language.
|
||||||
* :func:`~math.factorial` computes the factorial of a number.
|
* :func:`~math.factorial` computes the factorial of a number.
|
||||||
(Contributed by Raymond Hettinger; :issue:`2138`.)
|
(Contributed by Raymond Hettinger; :issue:`2138`.)
|
||||||
|
|
||||||
* :func:`~math.sum` adds up the stream of numbers from an iterable,
|
* :func:`~math.fsum` adds up the stream of numbers from an iterable,
|
||||||
and is careful to avoid loss of precision by calculating partial sums.
|
and is careful to avoid loss of precision by calculating partial sums.
|
||||||
(Contributed by Jean Brouwers, Raymond Hettinger, and Mark Dickinson;
|
(Contributed by Jean Brouwers, Raymond Hettinger, and Mark Dickinson;
|
||||||
:issue:`2819`.)
|
:issue:`2819`.)
|
||||||
|
|
|
||||||
|
|
@ -359,6 +359,102 @@ class MathTests(unittest.TestCase):
|
||||||
self.assertEquals(math.frexp(NINF)[0], NINF)
|
self.assertEquals(math.frexp(NINF)[0], NINF)
|
||||||
self.assert_(math.isnan(math.frexp(NAN)[0]))
|
self.assert_(math.isnan(math.frexp(NAN)[0]))
|
||||||
|
|
||||||
|
def testFsum(self):
|
||||||
|
# math.fsum relies on exact rounding for correct operation.
|
||||||
|
# There's a known problem with IA32 floating-point that causes
|
||||||
|
# inexact rounding in some situations, and will cause the
|
||||||
|
# math.fsum tests below to fail; see issue #2937. On non IEEE
|
||||||
|
# 754 platforms, and on IEEE 754 platforms that exhibit the
|
||||||
|
# problem described in issue #2937, we simply skip the whole
|
||||||
|
# test.
|
||||||
|
|
||||||
|
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.0 != 1e16+2.9999:
|
||||||
|
return
|
||||||
|
|
||||||
|
# Python version of math.fsum, for comparison. Uses a
|
||||||
|
# different algorithm based on frexp, ldexp and integer
|
||||||
|
# arithmetic.
|
||||||
|
from sys import float_info
|
||||||
|
mant_dig = float_info.mant_dig
|
||||||
|
etiny = float_info.min_exp - mant_dig
|
||||||
|
|
||||||
|
def msum(iterable):
|
||||||
|
"""Full precision summation. Compute sum(iterable) without any
|
||||||
|
intermediate accumulation of error. Based on the 'lsum' function
|
||||||
|
at http://code.activestate.com/recipes/393090/
|
||||||
|
|
||||||
|
"""
|
||||||
|
tmant, texp = 0, 0
|
||||||
|
for x in iterable:
|
||||||
|
mant, exp = math.frexp(x)
|
||||||
|
mant, exp = int(math.ldexp(mant, mant_dig)), exp - mant_dig
|
||||||
|
if texp > exp:
|
||||||
|
tmant <<= texp-exp
|
||||||
|
texp = exp
|
||||||
|
else:
|
||||||
|
mant <<= exp-texp
|
||||||
|
tmant += mant
|
||||||
|
# Round tmant * 2**texp to a float. The original recipe
|
||||||
|
# used float(str(tmant)) * 2.0**texp for this, but that's
|
||||||
|
# a little unsafe because str -> float conversion can't be
|
||||||
|
# relied upon to do correct rounding on all platforms.
|
||||||
|
tail = max(len(bin(abs(tmant)))-2 - mant_dig, etiny - texp)
|
||||||
|
if tail > 0:
|
||||||
|
h = 1 << (tail-1)
|
||||||
|
tmant = tmant // (2*h) + bool(tmant & h and tmant & 3*h-1)
|
||||||
|
texp += tail
|
||||||
|
return math.ldexp(tmant, texp)
|
||||||
|
|
||||||
|
test_values = [
|
||||||
|
([], 0.0),
|
||||||
|
([0.0], 0.0),
|
||||||
|
([1e100, 1.0, -1e100, 1e-100, 1e50, -1.0, -1e50], 1e-100),
|
||||||
|
([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),
|
||||||
|
([1./n for n in range(1, 1001)],
|
||||||
|
float.fromhex('0x1.df11f45f4e61ap+2')),
|
||||||
|
([(-1.)**n/n for n in range(1, 1001)],
|
||||||
|
float.fromhex('-0x1.62a2af1bd3624p-1')),
|
||||||
|
([1.7**(i+1)-1.7**i for i in range(1000)] + [-1.7**1000], -1.0),
|
||||||
|
([1e16, 1., 1e-16], 10000000000000002.0),
|
||||||
|
([1e16-2., 1.-2.**-53, -(1e16-2.), -(1.-2.**-53)], 0.0),
|
||||||
|
# exercise code for resizing partials array
|
||||||
|
([2.**n - 2.**(n+50) + 2.**(n+52) for n in range(-1074, 972, 2)] +
|
||||||
|
[-2.**1022],
|
||||||
|
float.fromhex('0x1.5555555555555p+970')),
|
||||||
|
]
|
||||||
|
|
||||||
|
for i, (vals, expected) in enumerate(test_values):
|
||||||
|
try:
|
||||||
|
actual = math.fsum(vals)
|
||||||
|
except OverflowError:
|
||||||
|
self.fail("test %d failed: got OverflowError, expected %r "
|
||||||
|
"for math.fsum(%.100r)" % (i, expected, vals))
|
||||||
|
except ValueError:
|
||||||
|
self.fail("test %d failed: got ValueError, expected %r "
|
||||||
|
"for math.fsum(%.100r)" % (i, expected, vals))
|
||||||
|
self.assertEqual(actual, expected)
|
||||||
|
|
||||||
|
from random import random, gauss, shuffle
|
||||||
|
for j in range(1000):
|
||||||
|
vals = [7, 1e100, -7, -1e100, -9e-20, 8e-20] * 10
|
||||||
|
s = 0
|
||||||
|
for i in range(200):
|
||||||
|
v = gauss(0, random()) ** 7 - s
|
||||||
|
s += v
|
||||||
|
vals.append(v)
|
||||||
|
shuffle(vals)
|
||||||
|
|
||||||
|
s = msum(vals)
|
||||||
|
self.assertEqual(msum(vals), math.fsum(vals))
|
||||||
|
|
||||||
def testHypot(self):
|
def testHypot(self):
|
||||||
self.assertRaises(TypeError, math.hypot)
|
self.assertRaises(TypeError, math.hypot)
|
||||||
self.ftest('hypot(0,0)', math.hypot(0,0), 0)
|
self.ftest('hypot(0,0)', math.hypot(0,0), 0)
|
||||||
|
|
@ -641,158 +737,6 @@ class MathTests(unittest.TestCase):
|
||||||
self.assertRaises(ValueError, math.sqrt, NINF)
|
self.assertRaises(ValueError, math.sqrt, NINF)
|
||||||
self.assert_(math.isnan(math.sqrt(NAN)))
|
self.assert_(math.isnan(math.sqrt(NAN)))
|
||||||
|
|
||||||
def testSum(self):
|
|
||||||
# math.sum relies on exact rounding for correct operation.
|
|
||||||
# There's a known problem with IA32 floating-point that causes
|
|
||||||
# inexact rounding in some situations, and will cause the
|
|
||||||
# math.sum tests below to fail; see issue #2937. On non IEEE
|
|
||||||
# 754 platforms, and on IEEE 754 platforms that exhibit the
|
|
||||||
# problem described in issue #2937, we simply skip the whole
|
|
||||||
# test.
|
|
||||||
|
|
||||||
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)
|
|
||||||
|
|
||||||
test_values = [
|
|
||||||
([], 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.**-53], OverflowError),
|
|
||||||
([1./n for n in range(1, 1001)], 7.4854708605503451),
|
|
||||||
([(-1.)**n/n for n in range(1, 1001)], -0.69264743055982025),
|
|
||||||
([1.7**(i+1)-1.7**i for i in range(1000)] + [-1.7**1000], -1.0),
|
|
||||||
([INF, -INF, NAN], ValueError),
|
|
||||||
([NAN, INF, -INF], ValueError),
|
|
||||||
([INF, NAN, INF], ValueError),
|
|
||||||
|
|
||||||
([INF, INF], OverflowError),
|
|
||||||
([INF, -INF], ValueError),
|
|
||||||
([-INF, 1e308, 1e308, -INF], OverflowError),
|
|
||||||
([2.0**1023-2.0**970, 0.0, 2.0**1023], OverflowError),
|
|
||||||
([2.0**1023-2.0**970, 1.0, 2.0**1023], OverflowError),
|
|
||||||
([2.0**1023, 2.0**1023], OverflowError),
|
|
||||||
([2.0**1023, 2.0**1023, -1.0], OverflowError),
|
|
||||||
([twopow, twopow, twopow, twopow, -twopow, -twopow],
|
|
||||||
OverflowError),
|
|
||||||
([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),
|
|
||||||
([-2.**1023, 2.**971, -2.**1023], -maxfloat),
|
|
||||||
([-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),
|
|
||||||
]
|
|
||||||
|
|
||||||
for i, (vals, s) in enumerate(test_values):
|
|
||||||
if isinstance(s, type) and issubclass(s, Exception):
|
|
||||||
try:
|
|
||||||
m = math.sum(vals)
|
|
||||||
except s:
|
|
||||||
pass
|
|
||||||
else:
|
|
||||||
self.fail("test %d failed: got %r, expected %r "
|
|
||||||
"for math.sum(%.100r)" %
|
|
||||||
(i, m, s.__name__, vals))
|
|
||||||
else:
|
|
||||||
try:
|
|
||||||
self.assertEqual(math.sum(vals), s)
|
|
||||||
except OverflowError:
|
|
||||||
self.fail("test %d failed: got OverflowError, expected %r "
|
|
||||||
"for math.sum(%.100r)" % (i, s, vals))
|
|
||||||
except ValueError:
|
|
||||||
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 range(1000):
|
|
||||||
vals = [7, 1e100, -7, -1e100, -9e-20, 8e-20] * 10
|
|
||||||
s = 0
|
|
||||||
for i in range(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):
|
def testTan(self):
|
||||||
self.assertRaises(TypeError, math.tan)
|
self.assertRaises(TypeError, math.tan)
|
||||||
self.ftest('tan(0)', math.tan(0), 0)
|
self.ftest('tan(0)', math.tan(0), 0)
|
||||||
|
|
|
||||||
|
|
@ -5,7 +5,7 @@ import random
|
||||||
import time
|
import time
|
||||||
import pickle
|
import pickle
|
||||||
import warnings
|
import warnings
|
||||||
from math import log, exp, sqrt, pi, sum as msum
|
from math import log, exp, sqrt, pi, fsum as msum
|
||||||
from test import support
|
from test import support
|
||||||
|
|
||||||
class TestBasicOps(unittest.TestCase):
|
class TestBasicOps(unittest.TestCase):
|
||||||
|
|
|
||||||
|
|
@ -396,7 +396,7 @@ FUNC1(tanh, tanh, 0,
|
||||||
Note 4: A similar implementation is in Modules/cmathmodule.c.
|
Note 4: A similar implementation is in Modules/cmathmodule.c.
|
||||||
Be sure to update both when making changes.
|
Be sure to update both when making changes.
|
||||||
|
|
||||||
Note 5: The signature of math.sum() differs from __builtin__.sum()
|
Note 5: The signature of math.fsum() differs from __builtin__.sum()
|
||||||
because the start argument doesn't make sense in the context of
|
because the start argument doesn't make sense in the context of
|
||||||
accurate summation. Since the partials table is collapsed before
|
accurate summation. Since the partials table is collapsed before
|
||||||
returning a result, sum(seq2, start=sum(seq1)) may not equal the
|
returning a result, sum(seq2, start=sum(seq1)) may not equal the
|
||||||
|
|
@ -407,7 +407,7 @@ FUNC1(tanh, tanh, 0,
|
||||||
|
|
||||||
/* Extend the partials array p[] by doubling its size. */
|
/* Extend the partials array p[] by doubling its size. */
|
||||||
static int /* non-zero on error */
|
static int /* non-zero on error */
|
||||||
_sum_realloc(double **p_ptr, Py_ssize_t n,
|
_fsum_realloc(double **p_ptr, Py_ssize_t n,
|
||||||
double *ps, Py_ssize_t *m_ptr)
|
double *ps, Py_ssize_t *m_ptr)
|
||||||
{
|
{
|
||||||
void *v = NULL;
|
void *v = NULL;
|
||||||
|
|
@ -425,7 +425,7 @@ _sum_realloc(double **p_ptr, Py_ssize_t n,
|
||||||
v = PyMem_Realloc(p, sizeof(double) * m);
|
v = PyMem_Realloc(p, sizeof(double) * m);
|
||||||
}
|
}
|
||||||
if (v == NULL) { /* size overflow or no memory */
|
if (v == NULL) { /* size overflow or no memory */
|
||||||
PyErr_SetString(PyExc_MemoryError, "math sum partials");
|
PyErr_SetString(PyExc_MemoryError, "math.fsum partials");
|
||||||
return 1;
|
return 1;
|
||||||
}
|
}
|
||||||
*p_ptr = (double*) v;
|
*p_ptr = (double*) v;
|
||||||
|
|
@ -464,18 +464,19 @@ _sum_realloc(double **p_ptr, Py_ssize_t n,
|
||||||
*/
|
*/
|
||||||
|
|
||||||
static PyObject*
|
static PyObject*
|
||||||
math_sum(PyObject *self, PyObject *seq)
|
math_fsum(PyObject *self, PyObject *seq)
|
||||||
{
|
{
|
||||||
PyObject *item, *iter, *sum = NULL;
|
PyObject *item, *iter, *sum = NULL;
|
||||||
Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
|
Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
|
||||||
double x, y, t, ps[NUM_PARTIALS], *p = ps;
|
double x, y, t, ps[NUM_PARTIALS], *p = ps;
|
||||||
|
double xsave, special_sum = 0.0, inf_sum = 0.0;
|
||||||
volatile double hi, yr, lo;
|
volatile double hi, yr, lo;
|
||||||
|
|
||||||
iter = PyObject_GetIter(seq);
|
iter = PyObject_GetIter(seq);
|
||||||
if (iter == NULL)
|
if (iter == NULL)
|
||||||
return NULL;
|
return NULL;
|
||||||
|
|
||||||
PyFPE_START_PROTECT("sum", Py_DECREF(iter); return NULL)
|
PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL)
|
||||||
|
|
||||||
for(;;) { /* for x in iterable */
|
for(;;) { /* for x in iterable */
|
||||||
assert(0 <= n && n <= m);
|
assert(0 <= n && n <= m);
|
||||||
|
|
@ -485,18 +486,19 @@ math_sum(PyObject *self, PyObject *seq)
|
||||||
item = PyIter_Next(iter);
|
item = PyIter_Next(iter);
|
||||||
if (item == NULL) {
|
if (item == NULL) {
|
||||||
if (PyErr_Occurred())
|
if (PyErr_Occurred())
|
||||||
goto _sum_error;
|
goto _fsum_error;
|
||||||
break;
|
break;
|
||||||
}
|
}
|
||||||
x = PyFloat_AsDouble(item);
|
x = PyFloat_AsDouble(item);
|
||||||
Py_DECREF(item);
|
Py_DECREF(item);
|
||||||
if (PyErr_Occurred())
|
if (PyErr_Occurred())
|
||||||
goto _sum_error;
|
goto _fsum_error;
|
||||||
|
|
||||||
|
xsave = x;
|
||||||
for (i = j = 0; j < n; j++) { /* for y in partials */
|
for (i = j = 0; j < n; j++) { /* for y in partials */
|
||||||
y = p[j];
|
y = p[j];
|
||||||
if (fabs(x) < fabs(y)) {
|
if (fabs(x) < fabs(y)) {
|
||||||
t = x; x = y; y = t;
|
t = x; x = y; y = t;
|
||||||
}
|
}
|
||||||
hi = x + y;
|
hi = x + y;
|
||||||
yr = hi - x;
|
yr = hi - x;
|
||||||
|
|
@ -505,59 +507,73 @@ math_sum(PyObject *self, PyObject *seq)
|
||||||
p[i++] = lo;
|
p[i++] = lo;
|
||||||
x = hi;
|
x = hi;
|
||||||
}
|
}
|
||||||
|
|
||||||
n = i; /* ps[i:] = [x] */
|
n = i; /* ps[i:] = [x] */
|
||||||
if (x != 0.0) {
|
if (x != 0.0) {
|
||||||
/* If non-finite, reset partials, effectively
|
if (! Py_IS_FINITE(x)) {
|
||||||
adding subsequent items without roundoff
|
/* a nonfinite x could arise either as
|
||||||
and yielding correct non-finite results,
|
a result of intermediate overflow, or
|
||||||
provided IEEE 754 rules are observed */
|
as a result of a nan or inf in the
|
||||||
if (! Py_IS_FINITE(x))
|
summands */
|
||||||
|
if (Py_IS_FINITE(xsave)) {
|
||||||
|
PyErr_SetString(PyExc_OverflowError,
|
||||||
|
"intermediate overflow in fsum");
|
||||||
|
goto _fsum_error;
|
||||||
|
}
|
||||||
|
if (Py_IS_INFINITY(xsave))
|
||||||
|
inf_sum += xsave;
|
||||||
|
special_sum += xsave;
|
||||||
|
/* reset partials */
|
||||||
n = 0;
|
n = 0;
|
||||||
else if (n >= m && _sum_realloc(&p, n, ps, &m))
|
}
|
||||||
goto _sum_error;
|
else if (n >= m && _fsum_realloc(&p, n, ps, &m))
|
||||||
p[n++] = x;
|
goto _fsum_error;
|
||||||
|
else
|
||||||
|
p[n++] = x;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
if (special_sum != 0.0) {
|
||||||
|
if (Py_IS_NAN(inf_sum))
|
||||||
|
PyErr_SetString(PyExc_ValueError,
|
||||||
|
"-inf + inf in fsum");
|
||||||
|
else
|
||||||
|
sum = PyFloat_FromDouble(special_sum);
|
||||||
|
goto _fsum_error;
|
||||||
|
}
|
||||||
|
|
||||||
hi = 0.0;
|
hi = 0.0;
|
||||||
if (n > 0) {
|
if (n > 0) {
|
||||||
hi = p[--n];
|
hi = p[--n];
|
||||||
if (Py_IS_FINITE(hi)) {
|
/* sum_exact(ps, hi) from the top, stop when the sum becomes
|
||||||
/* sum_exact(ps, hi) from the top, stop when the sum becomes inexact. */
|
inexact. */
|
||||||
while (n > 0) {
|
while (n > 0) {
|
||||||
x = hi;
|
x = hi;
|
||||||
y = p[--n];
|
y = p[--n];
|
||||||
assert(fabs(y) < fabs(x));
|
assert(fabs(y) < fabs(x));
|
||||||
hi = x + y;
|
hi = x + y;
|
||||||
yr = hi - x;
|
yr = hi - x;
|
||||||
lo = y - yr;
|
lo = y - yr;
|
||||||
if (lo != 0.0)
|
if (lo != 0.0)
|
||||||
break;
|
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;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
else { /* raise exception corresponding to a special value */
|
/* Make half-even rounding work across multiple partials.
|
||||||
errno = Py_IS_NAN(hi) ? EDOM : ERANGE;
|
Needed so that sum([1e-16, 1, 1e16]) will round-up the last
|
||||||
if (is_error(hi))
|
digit to two instead of down to zero (the 1e-16 makes the 1
|
||||||
goto _sum_error;
|
slightly closer to two). With a potential 1 ULP rounding
|
||||||
|
error fixed-up, math.fsum() 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;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
sum = PyFloat_FromDouble(hi);
|
sum = PyFloat_FromDouble(hi);
|
||||||
|
|
||||||
_sum_error:
|
_fsum_error:
|
||||||
PyFPE_END_PROTECT(hi)
|
PyFPE_END_PROTECT(hi)
|
||||||
Py_DECREF(iter);
|
Py_DECREF(iter);
|
||||||
if (p != ps)
|
if (p != ps)
|
||||||
|
|
@ -567,7 +583,7 @@ _sum_error:
|
||||||
|
|
||||||
#undef NUM_PARTIALS
|
#undef NUM_PARTIALS
|
||||||
|
|
||||||
PyDoc_STRVAR(math_sum_doc,
|
PyDoc_STRVAR(math_fsum_doc,
|
||||||
"sum(iterable)\n\n\
|
"sum(iterable)\n\n\
|
||||||
Return an accurate floating point sum of values in the iterable.\n\
|
Return an accurate floating point sum of values in the iterable.\n\
|
||||||
Assumes IEEE-754 floating point arithmetic.");
|
Assumes IEEE-754 floating point arithmetic.");
|
||||||
|
|
@ -1078,6 +1094,7 @@ static PyMethodDef math_methods[] = {
|
||||||
{"floor", math_floor, METH_O, math_floor_doc},
|
{"floor", math_floor, METH_O, math_floor_doc},
|
||||||
{"fmod", math_fmod, METH_VARARGS, math_fmod_doc},
|
{"fmod", math_fmod, METH_VARARGS, math_fmod_doc},
|
||||||
{"frexp", math_frexp, METH_O, math_frexp_doc},
|
{"frexp", math_frexp, METH_O, math_frexp_doc},
|
||||||
|
{"fsum", math_fsum, METH_O, math_fsum_doc},
|
||||||
{"hypot", math_hypot, METH_VARARGS, math_hypot_doc},
|
{"hypot", math_hypot, METH_VARARGS, math_hypot_doc},
|
||||||
{"isinf", math_isinf, METH_O, math_isinf_doc},
|
{"isinf", math_isinf, METH_O, math_isinf_doc},
|
||||||
{"isnan", math_isnan, METH_O, math_isnan_doc},
|
{"isnan", math_isnan, METH_O, math_isnan_doc},
|
||||||
|
|
@ -1091,10 +1108,9 @@ static PyMethodDef math_methods[] = {
|
||||||
{"sin", math_sin, METH_O, math_sin_doc},
|
{"sin", math_sin, METH_O, math_sin_doc},
|
||||||
{"sinh", math_sinh, METH_O, math_sinh_doc},
|
{"sinh", math_sinh, METH_O, math_sinh_doc},
|
||||||
{"sqrt", math_sqrt, METH_O, math_sqrt_doc},
|
{"sqrt", math_sqrt, METH_O, math_sqrt_doc},
|
||||||
{"sum", math_sum, METH_O, math_sum_doc},
|
|
||||||
{"tan", math_tan, METH_O, math_tan_doc},
|
{"tan", math_tan, METH_O, math_tan_doc},
|
||||||
{"tanh", math_tanh, METH_O, math_tanh_doc},
|
{"tanh", math_tanh, METH_O, math_tanh_doc},
|
||||||
{"trunc", math_trunc, METH_O, math_trunc_doc},
|
{"trunc", math_trunc, METH_O, math_trunc_doc},
|
||||||
{NULL, NULL} /* sentinel */
|
{NULL, NULL} /* sentinel */
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue