Shewchuk expansion arithmetic in Python
Represents results exactly as nonoverlapping expansions of floats, with adaptive precision.
- File
shewchuk_expansion_sum.py, 101 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/shewchuk-expansion-sum/. The authors are listed in prov/A734Z5.jsonl.
"""Shewchuk's expansion arithmetic for exact summation.
An expansion is a list of nonzero binary64 numbers, smallest magnitude first,
whose bits do not overlap. Its value is the exact sum of its components.
grow_expansion adds one float exactly; shewchuk_sum adds a sequence exactly
and rounds the result correctly, which is the method behind math.fsum.
"""
import math
def two_sum(a, b):
"""s = fl(a + b) and t = a + b - s exactly (Knuth, Moller)."""
s = a + b
bv = s - a
av = s - bv
t = (a - av) + (b - bv)
return s, t
def fast_two_sum(a, b):
"""Requires |a| >= |b| (or a == 0); then s + t = a + b exactly."""
s = a + b
t = b - (s - a)
return s, t
def grow_expansion(e, x):
"""Return the expansion of value(e) + x (Shewchuk, Grow-Expansion)."""
Q = x # @line L3
h = []
for y in e: # @line L4
Q, t = two_sum(Q, y) # @line L5
if t != 0.0: # @line L6
h.append(t)
if Q != 0.0: # @line L7
h.append(Q)
return h
def expansion_sum(e, f):
"""Exact sum of two expansions, by growing e with each component of f."""
for y in f:
e = grow_expansion(e, y)
return e
def round_expansion(e):
"""Correctly rounded (to nearest, ties to even) value of an expansion."""
if not e:
return 0.0
j = len(e) - 1
s = e[j] # @line L8
t = 0.0
while j > 0: # @line L9
j -= 1
s, t = fast_two_sum(s, e[j])
if t != 0.0:
break
if t != 0.0 and j > 0 and (t < 0) == (e[j - 1] < 0): # @line L10
y = 2.0 * t
x = s + y
if x - s == y and math.isfinite(x):
s = x
return s # @line L11
def shewchuk_sum(xs):
"""Correctly rounded sum of finite or special floats, like math.fsum."""
e = [] # @line L1
special = 0.0
inf_seen = False
for x in xs: # @line L2
if not math.isfinite(x):
inf_seen = True
special += x
continue
e = grow_expansion(e, x)
if e and not math.isfinite(e[-1]):
raise OverflowError("intermediate overflow in shewchuk_sum")
if inf_seen:
return special
return round_expansion(e)
def value(e):
"""Exact value of an expansion as a Fraction (for checking)."""
from fractions import Fraction
return sum((Fraction(c) for c in e), Fraction(0))
if __name__ == "__main__":
for xs in ([0.1] * 10, [1e100, 1.0, -1e100], [1.0, 2.0 ** -53, 2.0 ** -106]):
print(xs, sum(xs), shewchuk_sum(xs), math.fsum(xs))
e = []
for x in [0.1] * 10:
e = grow_expansion(e, x)
print(e)
test_shewchuk_expansion_sum.py 117 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/shewchuk-expansion-sum/. The authors are listed in prov/A734Z5.jsonl.
import math
import random
from fractions import Fraction
from shewchuk_expansion_sum import (expansion_sum, grow_expansion, round_expansion,
shewchuk_sum, two_sum, value)
def lowest_bit(x):
n, d = abs(x).as_integer_ratio()
return Fraction(n & -n, d)
def nonoverlapping_increasing(e):
for a, b in zip(e, e[1:]):
if a == 0.0 or b == 0.0:
return False
if not abs(a) < lowest_bit(b):
return False
return True
def rand_float(rng, lo=-60, hi=60):
return rng.choice((1, -1)) * math.ldexp(rng.random() + 1.0, rng.randint(lo, hi))
def test_two_sum_exact():
rng = random.Random(1)
for _ in range(3000):
a, b = rand_float(rng), rand_float(rng)
s, t = two_sum(a, b)
assert Fraction(s) + Fraction(t) == Fraction(a) + Fraction(b)
def test_grow_expansion_invariants():
rng = random.Random(2)
for trial in range(200):
e, total = [], Fraction(0)
for _ in range(rng.randint(1, 60)):
x = rand_float(rng, -80, 80)
e = grow_expansion(e, x)
total += Fraction(x)
assert value(e) == total
assert nonoverlapping_increasing(e)
def test_expansion_sum_exact():
rng = random.Random(3)
for trial in range(200):
e = f = []
for _ in range(10):
e = grow_expansion(e, rand_float(rng))
for _ in range(10):
f = grow_expansion(f, rand_float(rng))
g = expansion_sum(e, f)
assert value(g) == value(e) + value(f)
assert nonoverlapping_increasing(g)
def test_correct_rounding_random():
rng = random.Random(4)
for trial in range(1500):
xs = [rand_float(rng, -40, 40) for _ in range(rng.randint(1, 30))]
exact = sum((Fraction(x) for x in xs), Fraction(0))
assert shewchuk_sum(xs) == float(exact)
assert shewchuk_sum(xs) == math.fsum(xs)
def test_ties_and_near_ties():
u = 2.0 ** -53
cases = [
[1.0, u], # tie, stays at 1 (even)
[1.0, u, u * u], # just above the tie
[1.0, u, -u * u], # just below the tie
[1.0 + 2 * u, u], # tie, goes up to the even neighbour
[1.0, -u / 2], # tie below a power of two
[1.0, -u / 2, -2.0 ** -200],
[1e308, -1e308, 1e308, 2.0 ** -1074],
]
for xs in cases:
exact = sum((Fraction(x) for x in xs), Fraction(0))
assert shewchuk_sum(xs) == float(exact), xs
def test_catastrophic_cancellation():
assert shewchuk_sum([1e100, 1.0, -1e100]) == 1.0
assert shewchuk_sum([0.1, 0.2, 0.3]) == 0.6
rng = random.Random(5)
for trial in range(100):
xs = [rand_float(rng, -30, 30) for _ in range(50)]
ys = xs + [-x for x in xs] + [1.5]
rng.shuffle(ys)
assert shewchuk_sum(ys) == 1.5
def test_subnormal_inputs():
rng = random.Random(6)
for trial in range(300):
xs = [rng.randint(-10 ** 6, 10 ** 6) * 5e-324 for _ in range(10)] + [rand_float(rng, -1074, -1000)]
exact = sum((Fraction(x) for x in xs), Fraction(0))
assert shewchuk_sum(xs) == float(exact)
def test_specials():
assert shewchuk_sum([]) == 0.0
assert shewchuk_sum([1.0, math.inf]) == math.inf
assert math.isnan(shewchuk_sum([math.inf, -math.inf]))
try:
shewchuk_sum([1.7e308, 1.7e308])
except OverflowError:
pass
else:
raise AssertionError("expected OverflowError")
assert round_expansion([]) == 0.0
What the next level needs
- the tests did not run in this build