[Python-Dev] sum()

[Raymond Hettinger]

[Tim]
> For contrast, here's one that doesn't use frexp(), and is probably
> faster because of that; internally, len(sums) probably won't exceed 5
> in real life (but could get as large as 2K for pathological inputs,
> spanning fp's full dynamic range):
> 
> def summer4(iterable):
>     sums = [0.0]
>     for x in iterable:
>         sums.append(x)
>         for i in xrange(len(sums)-2, -1, -1):
>             y = sums[i]
>             if abs(x) < abs(y):
>                 x, y = y, x
>             hi = x+y
>             lo = y - (hi - x)
>             if lo:
>                 sums[i+1] = lo
>             else:
>                 del sums[i+1]
>             x = hi
>         sums[0] = x
>     return sum(reversed(sums), 0.0)


Here's a rewrite with the partial sums ordered by increasing magnitude.
A cursor is used to write-out new, non-zero terms.  That makes the list
growth or shrinkage occur at the end of the list and it avoids reversed
indexing.  Those two changes made the code easier for me to follow.

Leaving off the 0.0 makes the routine generic.  It works equally well
with int, long, Decimal, float, and complex.


def summer5(iterable, start=0):
    "Binary or Decimal floating point summation accurate to full
precision."
    partials = []       # sorted, non-overlapping partial sums
    for x in iterable:
        i = 0           # cursor for writing-out new partials sums
        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, start)



The loop invariants for the list of partial sums are:

    assert partials == sorted(partials, key=abs)
    assert nonoverlapping(partials)

where a rough check for overlaps is:

def nonoverlapping(seqn,    offset=100):
    """Verify that sequence elements have no overlapping bits.

    Set offset to -Emin to handle the full range of floats.
    """
    cumv = 0L
    for elem in seqn:
        m, exp = frexp(abs(elem))
        v = int(m * 2 ** 53) * 2L ** (exp + offset)
        if cumv & v:
            return False
        cumv |= v
    return True


Raymond