[Python-Dev] sum()

[Raymond Hettinger]

> 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