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 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] return sum(partials, 0.0) from math import frexp def lsum(iterable): "Full precision summation using long integers for intermediate values" # Transform (exactly) a float to m * 2 ** e where m and e are integers. # Adjust (tmant,texp) and (mant,exp) to make texp the common exponent. # Given a common exponent, the mantissas can be summed directly. tmant, texp = 0L, 0 for x in iterable: mant, exp = frexp(x) mant, exp = long(mant * 2.0 ** 53), exp-53 if texp > exp: tmant <<= texp - exp texp = exp else: mant <<= exp - texp tmant += mant return float(str(tmant)) * 2.0 ** texp from itertools import imap def lsum_26(iterable): "Full precision summation using long integers for intermediate values" # Py2.6 and later version of lsum() relying on float.as_integer_ratio(). # Saves work by accumulating separate sums for each power-of-two # denominator and then doing the bits shifts on the combined totals. fractions = {} fractions_get = fractions.get for n, d in imap(float.as_integer_ratio, iterable): fractions[d] = fractions_get(d, 0) + n tn, td = 0, 1 for d, n in sorted(fractions.items()): while td < d: td <<= 1 tn <<= 1 tn += n return tn / td from decimal import getcontext, Decimal, Inexact getcontext().traps[Inexact] = True def dsum(iterable): "Full precision summation using Decimal objects for intermediate values" # Transform (exactly) a float to m * 2 ** e where m and e are integers. # Convert (mant, exp) to a Decimal and add to the cumulative sum. # If the precision is too small for exact conversion and addition, # then retry with a larger precision. total = Decimal(0) for x in iterable: mant, exp = frexp(x) mant, exp = int(mant * 2.0 ** 53), exp-53 while True: try: total += mant * Decimal(2) ** exp break except Inexact: getcontext().prec += 1 return float(total) from fractions import Fraction def frsum(iterable): "Full precision summation using fractions for itermediate values" return float(sum(imap(Fraction.from_float, iterable))) 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) assert msum(vals) == lsum(vals) == dsum(vals) == frsum(vals) == lsum_26(vals) print '.', print 'Tests Passed' if __name__ == '__main__': test()