TwoProduct (Dekker) in Python
Computes a product and its exact rounding error as an unevaluated sum of two floats.
- File
two_product.py, 87 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/two-product/. The authors are listed in prov/Q0TV83.jsonl.
"""Dekker's TwoProduct: x = fl(a*b) and the exact error y = a*b - x.
Exact in binary64 with round to nearest when no operation overflows and
e_a + e_b >= e_min + p - 1 = -970 (so the error is representable).
"""
import math
from fractions import Fraction
SPLITTER = 134217729.0 # 2^27 + 1 for binary64
def split(a: float) -> tuple[float, float]:
"""Veltkamp splitting: a = hi + lo with 26-bit parts."""
gamma = SPLITTER * a
delta = a - gamma
hi = gamma + delta
return hi, a - hi
def two_product(a: float, b: float) -> tuple[float, float]:
"""Dekker's product in 17 floating-point operations, no FMA needed."""
x = a * b # @line L1
ah, al = split(a) # @line L2
bh, bl = split(b) # @line L3
t1 = ah * bh - x # @line L4
t2 = t1 + ah * bl # @line L5
t3 = t2 + al * bh # @line L6
y = t3 + al * bl # @line L7
return x, y # @line L8
def _fma_fallback(a: float, b: float, c: float) -> float:
"""Correctly rounded a*b + c for finite inputs, via exact rationals."""
return float(Fraction(a) * Fraction(b) + Fraction(c))
fma = getattr(math, "fma", _fma_fallback)
def two_product_fma(a: float, b: float) -> tuple[float, float]:
"""The FMA variant (2MultFMA): two operations."""
x = a * b
return x, fma(a, b, -x)
def round_to_precision(x: Fraction, p: int, emin: int) -> Fraction:
"""Round x to nearest, ties to even, in binary with precision p, subnormals below 2^emin."""
if x == 0:
return Fraction(0)
sign = -1 if x < 0 else 1
m = abs(x)
e = m.numerator.bit_length() - m.denominator.bit_length()
if Fraction(2) ** e > m:
e -= 1
e = max(e, emin)
quantum = Fraction(2) ** (e - p + 1)
q, r = divmod(m, quantum)
q = int(q)
if r > quantum / 2 or (r == quantum / 2 and q % 2 == 1):
q += 1
return sign * q * quantum
def two_product_emulated(a: Fraction, b: Fraction, p: int, emin: int) -> tuple[Fraction, Fraction]:
"""Dekker's product in a simulated binary format of precision p."""
def rn(v: Fraction) -> Fraction:
return round_to_precision(v, p, emin)
c = 2 ** ((p + 1) // 2) + 1
def sp(v: Fraction) -> tuple[Fraction, Fraction]:
g = rn(c * v)
h = rn(g + rn(v - g))
return h, rn(v - h)
x = rn(a * b)
ah, al = sp(a)
bh, bl = sp(b)
t1 = rn(rn(ah * bh) - x)
t2 = rn(t1 + rn(ah * bl))
t3 = rn(t2 + rn(al * bh))
return x, rn(t3 + rn(al * bl))
test_two_product.py 80 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/two-product/. The authors are listed in prov/Q0TV83.jsonl.
import math
import random
from fractions import Fraction
from two_product import _fma_fallback, two_product, two_product_emulated, two_product_fma
def exact(a: float, b: float, x: float, y: float) -> bool:
return Fraction(x) + Fraction(y) == Fraction(a) * Fraction(b)
def test_examples():
assert two_product(0.1, 0.1) == (0.010000000000000002, -8.326672684688674e-19)
a = 1.0 + 2.0**-52
assert two_product(a, a) == (1.0 + 2.0**-51, 2.0**-104)
assert two_product(3.0, 1.0 / 3.0) == (1.0, -(2.0**-54))
assert two_product(2.0, 3.0) == (6.0, 0.0)
def test_random_exact_and_matches_fma():
rng = random.Random(1971)
for _ in range(20000):
a = math.ldexp(rng.uniform(-1, 1), rng.randint(-480, 480))
b = math.ldexp(rng.uniform(-1, 1), rng.randint(-480, 480))
x, y = two_product(a, b)
assert x == a * b
assert exact(a, b, x, y)
assert (x, y) == two_product_fma(a, b)
assert abs(Fraction(y)) <= abs(Fraction(x)) * Fraction(1, 2**53)
def test_fma_fallback_agrees():
rng = random.Random(7)
for _ in range(3000):
a, b, c = (rng.uniform(-1e3, 1e3) for _ in range(3))
assert _fma_fallback(a, b, c) == float(Fraction(a) * Fraction(b) + Fraction(c))
if hasattr(math, "fma"):
assert math.fma(a, b, c) == _fma_fallback(a, b, c)
def test_hard_significands():
vals = [1.0 + 2.0**-52, 2.0 - 2.0**-52, 1.5 + 2.0**-52, 1.0 + 2.0**-26, 1.0 + 2.0**-27, 1.0 + 2.0**-28]
for a in vals:
for b in vals:
for sa in (1.0, -1.0):
x, y = two_product(sa * a, b)
assert exact(sa * a, b, x, y)
def test_underflow_loses_error():
# e_a + e_b = -1000 < -970: the error 2^-1104 is below the smallest subnormal
a = 2.0**-500 * (1.0 + 2.0**-52)
x, y = two_product(a, a)
assert y == 0.0 and not exact(a, a, x, y)
# e_a + e_b = -960 >= -970: exact again, with a subnormal error
a = 2.0**-480 * (1.0 + 2.0**-52)
x, y = two_product(a, a)
assert exact(a, a, x, y) and y == 2.0**-1064
def test_overflow_of_split():
a = 2.0**1000
x, y = two_product(a, 2.0**-1000)
assert x == 1.0 and math.isnan(y)
x, y = two_product_fma(a, 2.0**-1000)
assert (x, y) == (1.0, 0.0)
def test_exhaustive_small_precision():
for p in (4, 5, 6):
emin = -30
sig = [Fraction(m, 2 ** (p - 1)) for m in range(2 ** (p - 1), 2**p)]
for a in sig:
for b in sig:
for sb in (1, -1):
x, y = two_product_emulated(a, sb * b, p, emin)
assert x + y == a * sb * b, (p, a, b)
What the next level needs
- the tests did not run in this build