[Python-Dev] sum() (original) (raw)

Raymond Hettinger python at rcn.com
Sat Mar 12 23:32:39 CET 2005

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

More information about the Python-Dev mailing list