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()