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.
This commit is contained in:
Mark Dickinson 2008-07-29 18:45:38 +00:00
parent cde8cdd898
commit ff5f16e4e5

View file

@ -659,48 +659,42 @@ class MathTests(unittest.TestCase):
# on IEEE 754 compliant machines, both of the expressions # on IEEE 754 compliant machines, both of the expressions
# below should round to 10000000000000002.0. # below should round to 10000000000000002.0.
if 1e16+2.999 != 1e16+2.9999: if 1e16+2.0 != 1e16+2.9999:
return return
# Python version of math.sum algorithm, for comparison # Python version of math.sum, 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): def msum(iterable):
"""Full precision sum of values in iterable. Returns the value of """Full precision summation. Compute sum(iterable) without any
the sum, rounded to the nearest representable floating-point number intermediate accumulation of error. Based on the 'lsum' function
using the round-half-to-even rule. at http://code.activestate.com/recipes/393090/
""" """
# Stage 1: accumulate partials tmant, texp = 0, 0
partials = []
for x in iterable: for x in iterable:
i = 0 mant, exp = math.frexp(x)
for y in partials: mant, exp = int(math.ldexp(mant, mant_dig)), exp - mant_dig
if abs(x) < abs(y): if texp > exp:
x, y = y, x tmant <<= texp-exp
hi = x + y texp = exp
lo = y - (hi - x) else:
if lo: mant <<= exp-texp
partials[i] = lo tmant += mant
i += 1 # Round tmant * 2**texp to a float. The original recipe
x = hi # used float(str(tmant)) * 2.0**texp for this, but that's
partials[i:] = [x] if x else [] # a little unsafe because str -> float conversion can't be
# relied upon to do correct rounding on all platforms.
# Stage 2: sum partials tail = max(len(bin(abs(tmant)))-2 - mant_dig, etiny - texp)
if not partials: if tail > 0:
return 0.0 h = 1 << (tail-1)
tmant = tmant // (2*h) + bool(tmant & h and tmant & 3*h-1)
# sum from the top, stopping as soon as the sum is inexact. texp += tail
total = partials.pop() return math.ldexp(tmant, texp)
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
test_values = [ test_values = [
([], 0.0), ([], 0.0),
@ -710,12 +704,18 @@ class MathTests(unittest.TestCase):
([2.0**53, 1.0, 2.0**-100], 2.0**53+2.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+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**53-4.0, 0.5, 2.0**-54], 2.0**53-3.0),
([1./n for n in range(1, 1001)], 7.4854708605503451), ([1./n for n in range(1, 1001)],
([(-1.)**n/n for n in range(1, 1001)], -0.69264743055982025), 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), ([1.7**(i+1)-1.7**i for i in range(1000)] + [-1.7**1000], -1.0),
([1e16, 1., 1e-16], 10000000000000002.0), ([1e16, 1., 1e-16], 10000000000000002.0),
([1e16-2., 1.-2.**-53, -(1e16-2.), -(1.-2.**-53)], 0.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): for i, (vals, expected) in enumerate(test_values):
try: try: