mdickinson
commented
1 year ago Very nice! It's surprising how little code it takes to make this work.
A few questions:
- If this went in, would it be considered a CPython implementation detail, or would other implementations be expected to follow suit?
- I'd expect at least some breakage in downstream unit tests (and perhaps doctests) from the changed numbers. People shouldn't be relying on down-to-the-last-bit reproducibility from
sumin unit tests, but I'd like to bet that at least some are, and in practice our currentsumdoes effectively give reproducible results across platforms, especially now that x87 FPUs and the associated double-rounding issues are fading into obscurity. What level of breakage would we consider acceptable here? - Did you try timing the variant that uses 2Sum in place of Fast2Sum for computing the individual errors? It involves more floating-point operations, but eliminates the branch in the absolute value test.
- Would we need to keep a way to get the naive, uncompensated sum? It's a need that I've occasionally had, and I've written code that relies on
sumdoing the obvious thing. I'll admit that it's likely not a common use-case, and people who need it can still dofunctools.reduce(operator.add, values, start).
On performance, I'm not seeing the "zero cost" on macOS / Intel. Here are some timings on my machine (comparing commit af7613ebcf1d118291320f07e9fab8100266c13c on this branch - "cpython" below - with commit e0b4d966a8d1867c4b535b043e08288ca49b3548 on main - "cpython-compare" below); both clones built with --enable-optimizations. tl;dr - this is around 25% slower on my machine.
mdickinson@lovelace python % ./cpython/python.exe -m timeit -s 'from random import expovariate as r' -s 'd=[r(1.0) for i in range(10_000)]' 'sum(d)'
5000 loops, best of 5: 44.9 usec per loop
mdickinson@lovelace python % ./cpython/python.exe -m timeit -s 'from random import expovariate as r' -s 'd=[r(1.0) for i in range(10_000)]' 'sum(d)'
5000 loops, best of 5: 47.6 usec per loop
mdickinson@lovelace python % ./cpython/python.exe -m timeit -s 'from random import expovariate as r' -s 'd=[r(1.0) for i in range(10_000)]' 'sum(d)'
5000 loops, best of 5: 45 usec per loop
mdickinson@lovelace python % ./cpython-compare/python.exe -m timeit -s 'from random import expovariate as r' -s 'd=[r(1.0) for i in range(10_000)]' 'sum(d)'
10000 loops, best of 5: 35.3 usec per loop
mdickinson@lovelace python % ./cpython-compare/python.exe -m timeit -s 'from random import expovariate as r' -s 'd=[r(1.0) for i in range(10_000)]' 'sum(d)'
10000 loops, best of 5: 34.5 usec per loop
mdickinson@lovelace python % ./cpython-compare/python.exe -m timeit -s 'from random import expovariate as r' -s 'd=[r(1.0) for i in range(10_000)]' 'sum(d)'
10000 loops, best of 5: 34.4 usec per loopHere's another test case that involves a decent amount of cancellation - we're computing (the imaginary part of) the Gauss sum for the prime 609359; mathematically, it should give exactly √609359. Here's the code:
import fractions
import math
import sys
import timeit
# Test case: Gauss sum for the prime 609359; mathematically, result
# should be exactly sqrt(p).
p = 609359
terms = [math.sin(2*math.pi*r*r / p) for r in range(p)]
# Exact sum of the actual terms we're using, and a way to compute the
# ulps error.
true_sum = sum(map(fractions.Fraction, terms))
def ulps_error(x: float) -> float:
return (fractions.Fraction(x) - true_sum) / math.ulp(true_sum)
# Error incurred by sum.
error = ulps_error(sum(terms))
# Timing
number = 1000
time = timeit.timeit(stmt="sum(terms)", globals=globals(), number=number)
print(sys.version)
print(f"error: {ulps_error(sum(terms))} ulps")
print(f"time: {1000.0 * time/number:.3f} ms per run")Results on my machine show around a 20% slowdown:
mdickinson@lovelace python % ./cpython-compare/python.exe ~/Desktop/compensated.py
3.12.0a3+ (tags/v3.12.0a3-114-ge0b4d966a8:e0b4d966a8, Dec 22 2022, 13:23:36) [Clang 14.0.0 (clang-1400.0.29.202)]
error: 58.52302999049425 ulps
time: 3.066 ms per run
mdickinson@lovelace python % ./cpython/python.exe ~/Desktop/compensated.py
3.12.0a3+ (heads/compensated_sum:af7613ebcf, Dec 22 2022, 13:34:38) [Clang 14.0.0 (clang-1400.0.29.202)]
error: -0.47697000950574875 ulps
time: 3.799 ms per run
mdickinson@lovelace python % ./cpython-compare/python.exe ~/Desktop/compensated.py
3.12.0a3+ (tags/v3.12.0a3-114-ge0b4d966a8:e0b4d966a8, Dec 22 2022, 13:23:36) [Clang 14.0.0 (clang-1400.0.29.202)]
error: 58.52302999049425 ulps
time: 3.063 ms per run
mdickinson@lovelace python % ./cpython/python.exe ~/Desktop/compensated.py
3.12.0a3+ (heads/compensated_sum:af7613ebcf, Dec 22 2022, 13:34:38) [Clang 14.0.0 (clang-1400.0.29.202)]
error: -0.47697000950574875 ulps
time: 3.789 ms per runOn pairwise summation: agreed that this (i.e., compensated summation) seems superior for our use-case, though it might be interesting to code up pairwise summation just to see what the performance tradeoffs look like. It doesn't need much in the way of auxiliary memory: at most a few hundred bytes for maintaining a stack of partial values of size O(log_2(length)). E.g., in Python:
def _ruler_seq():
"""Sequence 0, 1, 0, 2, 0, 1, 0, 3, ... (A007814)"""
i = 0
while True:
yield (~(i + 1) & i).bit_length()
i += 1
def sum_pairwise(iterable, /, start=0.0):
"""Sum of an iterable of floats, via pairwise summation."""
partials = [start]
for x, pop_count in zip(iterable, _ruler_seq()):
for _ in range(pop_count):
x += partials.pop()
partials.append(x)
return sum(reversed(partials))(The _ruler_seq logic would of course be much more streamlined in C.)
rhettinger
pochmann
CAM-Gerlach
hauntsaninja
gpshead
tim-one
franklinvp
skirpichev
pochmann3
Currently
sum()makes no efforts to improve accuracy over a simple running total. We do havemath.fsum()that makes extreme efforts to be almost perfect; however, that function isn't well known, it runs 10x slower than regularsum(), and it always coerces to a float.I suggest switching the builtin
sum()handling of float inputs to Arnold Neumaier's improved variation of compensated summation. Per his paper, this algorithm has excellent error bounds (though not as perfect asfsum():The compensation tracking runs in parallel to the main accumulation. And except for pathological cases, the branch is predictable, making the test essentially free. Thanks to the magic of instruction level parallelism and branch prediction, this improvement has zero cost on my Apple M1 build. Timed with:
./python.exe -m timeit -s 'from random import expovariate as r' -s 'd=[r(1.0) for i in range(10_000)]' 'sum(d)'N.B. Numpy switched from a simple running total to partial pairwise summation. That isn't as accurate as what is being proposed here, but it made more sense for them because the extra work of Neumaier summation isn't masked by the overhead of fetching values from an iterator as we do here. Also with an iterator, we can't do pairwise summation without using auxiliary memory.
Linked PRs