I would like to be able use fsum incrementally however it is not currently possible.

With sum() you can do:

subtotal = sum(nums)
subtotal = sum(othernums, subtotal)

This wouldn't work for fsum() because the returned float is not the same as the state maintained internally which is a list of floats. I propose instead that a fsum() could take an optional second argument which is a list of floats and in this case return the updated list of floats.

I have modified Raymond Hettinger's msum() recipe to do what I want:

def fsum(iterable, carry=None):
    "Full precision summation using multiple floats for intermediate values"
    partials = list(carry) if carry else []
    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 carry is None:
        return sum(partials, 0.0)
    else:
        return partials


Here's an interactive session showing how you might use it:

>>> from fsum import fsum
>>> fsum([1e20, 1, -1e20])
1.0
>>> fsum([1e20, 1, -1e20], [])
[1.0, 0.0]
>>> fsum([1e20, 1, -1e20], [])
[1.0, 0.0]
>>> fsum([1e20, 1, -1e20], [])
[1.0, 0.0]
>>> fsum([1e20, 1], [])
[1.0, 1e+20]
>>> carry = fsum([1e20, 1], [])
>>> fsum([-1e20], carry)
[1.0, 0.0]
>>> nums = [7, 1e100, -7, -1e100, -9e-20, 8e-20] * 10
>>> subtotal = []
>>> for n in nums:
...     subtotal = fsum([n], subtotal)
...
>>> subtotal
[-1.0000000000000007e-19]
>>> fsum(subtotal)
-1.0000000000000007e-19


My motivation for this is that I want to be able to write an efficient sum function that is accurate for a mix of ints and floats while being as fast as possible in the case that I have a list of only floats (or only ints). What I have so far looks a bit like:

def sum(numbers):
    exact_total = 0
    float_total = 0.0
    for T, nums in groupby(numbers):
        if T is int:
            exact_total = sum(nums, exact_total)
        elif T is float:
            # This doesn't really work:
            float_total += fsum(float_total)
        # ...


However fsum is only exact if it adds all the numbers in a single pass. The above would have order-dependent results given a mixed list of ints and floats e.g.:

[1e20, -1e20, 1, -1, 1.0, -1.0]
  vs
[1e20, 1.0, -1e20, 1, -1, -1.0]

Although fsum is internally very accurate it always discards information on output. Even if I build a list of all floats and fsum those at the end it can still be inaccurate if the exact_total cannot be exactly coerced to float and passed to fsum.

If I could get fsum to return the exact float expansion then I could use that in a number of different ways. Given the partials list I can combine all of the different subtotals with Fraction arithmetic and coerce to float only at the very end. If I can also get fsum to accept a float expansion on input then I can use it incrementally and there is no need to build a list of all floats.

I am prepared to write a patch for this if the idea is deemed acceptable.