peak_partials = 0 # peak number of partials
def msum(iterable):
"Full precision summation using multiple floats for intermediate values"
# Rounded x+y stored in hi with the round-off stored in lo. Together
# hi+lo are exactly equal to x+y. The inner loop applies hi/lo summation
# to each partial so that the list of partial sums remains exact.
# Depends on IEEE-754 arithmetic guarantees. See proof of correctness at:
# www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps
global peak_partials
partials = [] # sorted, non-overlapping partial sums
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 peak_partials < i:
peak_partials = i
return sum(partials, 0.0)
from math import sum as math_sum
try:
from cmath import sum as cmath_sum
except ImportError:
cmath_sum = None
from random import random, gauss, shuffle
def test():
for j in xrange(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)
# compare the C result with the result
# from Python function msum() above
s = msum(vals)
assert s == math_sum(vals)
if cmath_sum: # check complex too
c = cmath_sum([complex(v, -v) for v in vals])
assert s == c.real == -c.imag
# check that the full precision sum is
# different/better than the built-in sum
assert s != sum(vals)
##print '.', ## fails on MacOS X
if cmath_sum:
print 'tests passed, incl. complex', '(%d partials)' % peak_partials
else:
print 'tests passed', '(%d partials)' % peak_partials
if __name__ == '__main__':
test()