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. |