Pairwise summation in Python
Recursively sums halves, giving error growth logarithmic rather than linear in the number of terms.
- File
pairwise_summation.py, 82 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/pairwise-summation/. The authors are listed in prov/QGCJX6.jsonl.
"""Pairwise (cascade) summation of binary64 floats.
The terms are split into two halves, each half is summed the same way, and
the two partial sums are added. Every term then passes through at most
ceil(log2 n) additions, so the error bound grows like log2 n instead of n.
"""
from collections.abc import Sequence
def recursive_sum(xs: Sequence[float]) -> float:
"""Plain left-to-right summation, one rounding per term. Cost: n - 1 additions.
Python's built-in sum() uses compensated summation for floats since
Python 3.12, so it is not used here as the naive reference.
"""
s = 0.0
for x in xs:
s += x
return s
def pairwise_sum(xs: Sequence[float], block: int = 1) -> float:
"""Pairwise sum of xs; slices of at most `block` terms are summed left to right.
With block = 1 the tree is a full binary split down to single terms.
The recursion depth is ceil(log2(n / block)) + 1, so no deep recursion
occurs. Cost: n - 1 additions and O(log n) stack frames.
"""
if block < 1:
raise ValueError("block must be at least 1")
n = len(xs)
if n == 0:
return 0.0
return _sum(xs, 0, n, block)
def _sum(xs: Sequence[float], lo: int, hi: int, block: int) -> float:
n = hi - lo # @line L1
if n <= block: # @line L2
s = xs[lo] # @line L3
for i in range(lo + 1, hi): # @line L4
s += xs[i]
return s # @line L5
mid = lo + n // 2 # @line L6
s1 = _sum(xs, lo, mid, block) # @line L7
s2 = _sum(xs, mid, hi, block) # @line L8
return s1 + s2 # @line L9
def pairwise_sum_iterative(xs: Sequence[float]) -> float:
"""Bottom-up pairwise sum: add neighbours level by level.
Each level adds x[0]+x[1], x[2]+x[3], ...; an odd last term is carried
up unchanged. The tree has depth ceil(log2 n), the same bound as the
recursive split, though the tree shape can differ when n is not a power
of two. Cost: n - 1 additions and n/2 extra floats.
"""
level = list(xs)
if not level:
return 0.0
while len(level) > 1:
nxt = [level[i] + level[i + 1] for i in range(0, len(level) - 1, 2)]
if len(level) % 2:
nxt.append(level[-1])
level = nxt
return level[0]
def depth(n: int, block: int = 1) -> int:
"""The largest number of additions any term goes through in pairwise_sum.
It is the h in the error bound |error| <= gamma_h * sum |x_i|.
"""
if n <= 1:
return 0
if n <= block:
return n - 1
return 1 + depth(n - n // 2, block)
test_pairwise_summation.py 114 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/pairwise-summation/. The authors are listed in prov/QGCJX6.jsonl.
import math
import random
from fractions import Fraction
from pairwise_summation import depth, pairwise_sum, pairwise_sum_iterative, recursive_sum
U = Fraction(1, 2**53)
def gamma(h: int) -> Fraction:
return h * U / (1 - h * U)
def exact(xs) -> Fraction:
return sum((Fraction(x) for x in xs), Fraction(0))
def ceil_log2(n: int) -> int:
return (n - 1).bit_length()
def test_small_cases():
assert pairwise_sum([]) == 0.0
assert pairwise_sum_iterative([]) == 0.0
assert pairwise_sum([2.5]) == 2.5
assert pairwise_sum([0.1, 0.2]) == 0.1 + 0.2
assert pairwise_sum([1.0, 2.0, 3.0, 4.0]) == (1.0 + 2.0) + (3.0 + 4.0)
def test_worked_example():
xs = [1.0] + [2.0**-53] * 7
assert recursive_sum(xs) == 1.0
assert pairwise_sum(xs) == 1.0 + 3 * 2.0**-52
assert Fraction(pairwise_sum(xs)) == exact(xs) - Fraction(1, 2**53)
assert math.fsum(xs) == 1.0 + 2.0**-50
def test_repeated_tenth():
xs = [0.1] * 10**5
assert pairwise_sum(xs) == math.fsum(xs)
assert recursive_sum(xs) != math.fsum(xs)
def test_depth_is_ceil_log2():
for n in range(1, 3000):
assert depth(n) == ceil_log2(n)
for n in (1, 7, 8, 9, 100, 1000):
for b in (1, 2, 8, 128):
d = depth(n, b)
assert d <= (min(n, b) - 1) + max(0, ceil_log2(-(-n // b)))
def test_error_bound_random():
rng = random.Random(20240601)
for trial in range(200):
n = rng.randint(1, 600)
kind = trial % 4
if kind == 0:
xs = [rng.random() for _ in range(n)]
elif kind == 1:
xs = [rng.uniform(-1, 1) for _ in range(n)]
elif kind == 2:
xs = [math.ldexp(rng.uniform(-1, 1), rng.randint(-60, 60)) for _ in range(n)]
else:
xs = [rng.choice([1e16, -1e16, 1.0, 0.5, 3.0]) for _ in range(n)]
ex = exact(xs)
abssum = exact(abs(x) for x in xs)
for block in (1, 4, 32):
err = abs(Fraction(pairwise_sum(xs, block)) - ex)
assert err <= gamma(depth(n, block)) * abssum
err = abs(Fraction(pairwise_sum_iterative(xs)) - ex)
assert err <= gamma(ceil_log2(n)) * abssum
def test_exact_on_integers():
rng = random.Random(5)
xs = [float(rng.randint(-(2**40), 2**40)) for _ in range(1000)]
assert Fraction(pairwise_sum(xs)) == exact(xs)
def test_beats_recursive_on_large_positive_sum():
rng = random.Random(1)
xs = [rng.random() for _ in range(100000)]
ex = exact(xs)
e_rec = abs(Fraction(recursive_sum(xs)) - ex)
e_pw = abs(Fraction(pairwise_sum(xs)) - ex)
assert e_pw < e_rec
assert e_pw <= 2 * U * ex
def test_special_values():
inf, nan = math.inf, math.nan
assert pairwise_sum([1.0, inf, 2.0]) == inf
assert math.isnan(pairwise_sum([inf, -inf]))
assert math.isnan(pairwise_sum([1.0, nan]))
big = 1.7976931348623157e308
# (big + big) + (-big - big) = inf + (-inf): the exact sum 0 is lost to overflow
assert math.isnan(pairwise_sum([big, big, -big, -big]))
assert recursive_sum([big, big, -big, -big]) == inf
assert math.copysign(1.0, pairwise_sum([-0.0, -0.0])) == -1.0
def test_block_argument():
try:
pairwise_sum([1.0], 0)
except ValueError:
pass
else:
raise AssertionError("block 0 accepted")
xs = [0.1 * i for i in range(50)]
assert pairwise_sum(xs, 100) == recursive_sum(xs)
What the next level needs
- the tests did not run in this build