peak_partials = 0  # peak number of partials

def msum(iterable):
    "Full precision summation using multiple floats for intermediate values"

    # Rounded x+y stored in hi with the round-off stored in lo.  Together
    # hi+lo are exactly equal to x+y.  The inner loop applies hi/lo summation
    # to each partial so that the list of partial sums remains exact.

    # Depends on IEEE-754 arithmetic guarantees.  See proof of correctness at:
    # www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps

    global peak_partials
    partials = []  # sorted, non-overlapping partial sums
    for x in iterable:
        i = 0
        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]
        if peak_partials < i:
           peak_partials = i
    return sum(partials, 0.0)

from math import sum as math_sum
try:
    from cmath import sum as cmath_sum
except ImportError:
	cmath_sum = None

from random import random, gauss, shuffle

def test():

    for j in xrange(1000):
        vals = [7, 1e100, -7, -1e100, -9e-20, 8e-20] * 10
        s = 0
        for i in range(200):
            v = gauss(0, random()) ** 7 - s
            s += v
            vals.append(v)
        shuffle(vals)

         # compare the C result with the result
         # from Python function msum() above
        s = msum(vals)
        assert s == math_sum(vals)

        if cmath_sum:  # check complex too
            c = cmath_sum([complex(v, -v) for v in vals])
            assert s == c.real == -c.imag

         # check that the full precision sum is
         # different/better than the built-in sum
        assert s != sum(vals)

        ##print '.',  ## fails on MacOS X

    if cmath_sum:
    	print 'tests passed, incl. complex', '(%d partials)' % peak_partials
    else:
    	print 'tests passed', '(%d partials)' % peak_partials

if __name__ == '__main__':
    test()