[Python-Dev] sum()
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 <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
More information about the Python-Dev
mailing list