the queue can be expected to grow to about half the length of the original iterable by the time the original iterable is exhausted.
x = z mant, exp = frexp(x) exp2sum[exp] = x return sum(sorted(exp2sum.itervalues(), key=abs), 0.0)
The implementation can be tweaked to consume the error term right
away
so the queue won't grow by more than few pending items.
Theorem 10 in Shewchuk's "Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates" gives the obvious <wink> correct way to do that, although as a practical matter it greatly benifits from a bit more logic to eliminate zero entries as they're produced (Shewchuk doesn't because it's not his goal to add a bunch of same-precision floats). BTW, extracting binary exponents isn't necessary in his approach (largely specializations to "perfect" 754 arithmetic of Priest's earlier less-demanding methods).
The approach I'm using relies on being able to exactly multiply the 0 or 1 bit error term mantissas by an integer (a count of how many times the term occurs). With a Python dictionary keeping the counts, many steps are saved and the tool becomes much more memory friendly than with the previous queue approach. from math import frexp def summer3(iterable): exp2sum = {} # map to a value with a given exponent errdict = {} # map error terms to an occurrence count def addone(x, exppop=exp2sum.pop, errget=errdict.get, frexp=frexp): mant, exp = frexp(x) # extract exponent from float representation while exp in exp2sum: y = exppop(exp) # pair with a value having the same exp z = x + y # sum may be inexact by less than 1 ulp of z d = x - z + y # difference is the error term assert z + d == x + y # sum plus error term is exact! errdict[d] = errget(d, 0) + 1 # count occurrences of each term x = z # continue with new sum mant, exp = frexp(x) exp2sum[exp] = x for x in iterable: addone(x) while errdict: d, q = errdict.popitem() assert frexp(d)[0] in [-0.5, 0.0, 0.5] # error terms are 1 bit # product of 1 bit error term with an integer count is exact addone(d * q) dummy = errdict.pop(0.0, None) # At this point, the errdict is empty, all partial sums are exact, # and each partial sum has a different exponent. They can now be # added smallest absolute value to largest and only lose bits that # cannot fit in the final result. IOW, the final result is accurate # to full precision (assuming no representation error in the inputs). return sum(sorted(exp2sum.itervalues(), key=abs), 0.0) Aside from being a nice recipe, this was a good exercise in seeing what pure Python can do with IEEE-754 guarantees. The modest memory consumption and the typical O(n) runtime are a nice plus (no pun intended). Raymond