Dot2 compensated dot product in Python
Computes the dot product as if in twice the working precision using error-free transformations.
- File
dot2_compensated_dot_product.py, 69 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/dot2-compensated-dot-product/. The authors are listed in prov/7CXZPB.jsonl.
"""Dot2 (Ogita, Rump and Oishi, 2005): a dot product as accurate as if computed in twice the precision.
Every product is split exactly into p + e by TwoProduct, the running sum is split exactly into
p + q by TwoSum, and all the small parts are added in ordinary floating point.
"""
import math
from collections.abc import Sequence
SPLITTER = 134217729.0 # 2^27 + 1, Veltkamp's constant for binary64
def two_sum(a: float, b: float) -> tuple[float, float]:
"""Knuth's TwoSum: s = fl(a + b) and e with a + b = s + e exactly (no overflow). Cost: 6 flops."""
s = a + b
z = s - a
return s, (a - (s - z)) + (b - z)
def two_product_fma(a: float, b: float) -> tuple[float, float]:
"""p = fl(a b) and e = a b - p exactly, by one FMA (exact when there is no underflow). Cost: 2 flops."""
p = a * b
return p, math.fma(a, b, -p)
def _split(a: float) -> tuple[float, float]:
c = SPLITTER * a
hi = c - (c - a)
return hi, a - hi
def two_product_dekker(a: float, b: float) -> tuple[float, float]:
"""Dekker's TwoProduct with Veltkamp splitting, without FMA (exact when there is no
underflow and |a|, |b| < 2^996). Cost: 17 flops."""
p = a * b
ah, al = _split(a)
bh, bl = _split(b)
return p, al * bl - (((p - ah * bh) - al * bh) - ah * bl)
def dot2(x: Sequence[float], y: Sequence[float], fma: bool = True) -> float:
"""The dot product sum(x[i] * y[i]) with Dot2.
Precondition: len(x) == len(y) == n, finite entries, n u < 1, no underflow in the products
and no overflow. Then |result - x.y| <= u |x.y| + (1 + u) gamma_n^2 sum |x[i] y[i]| (proved).
Cost: 25n - 7 flops with fma=False, 10n - 8 with fma=True; O(1) extra space.
"""
if len(x) != len(y):
raise ValueError("vectors of different lengths")
if not x:
return 0.0
two_product = two_product_fma if fma else two_product_dekker
p, s = two_product(x[0], y[0]) # @line L1
for i in range(1, len(x)): # @line L2
h, r = two_product(x[i], y[i]) # @line L3
p, q = two_sum(p, h) # @line L4
s = s + (q + r) # @line L5
return p + s # @line L6
def naive_dot(x: Sequence[float], y: Sequence[float]) -> float:
"""Recursive dot product in working precision, for comparison. Cost: 2n flops."""
s = 0.0
for a, b in zip(x, y):
s += a * b
return s
test_dot2_compensated_dot_product.py 132 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/dot2-compensated-dot-product/. The authors are listed in prov/7CXZPB.jsonl.
"""Tests for Dot2: exact error-free transformations, the proved bound, and ill-conditioned cases."""
import math
import random
from fractions import Fraction
import pytest
from dot2_compensated_dot_product import dot2, naive_dot, two_product_dekker, two_product_fma, two_sum
U = Fraction(1, 2 ** 53)
def gamma(n: int) -> Fraction:
return n * U / (1 - n * U)
def exact(x: list[float], y: list[float]) -> Fraction:
return sum((Fraction(a) * Fraction(b) for a, b in zip(x, y)), Fraction(0))
def bound(x: list[float], y: list[float]) -> Fraction:
n = len(x)
return U * abs(exact(x, y)) + (1 + U) * gamma(n) ** 2 * sum(abs(Fraction(a) * Fraction(b)) for a, b in zip(x, y))
def gen_dot(n: int, log10_cond: float, rng: random.Random) -> tuple[list[float], list[float]]:
"""Ill-conditioned dot products in the spirit of Ogita, Rump and Oishi's GenDot (Algorithm 6.1)."""
n2 = n // 2
b = log10_cond * math.log2(10) / 2
e = [round(rng.random() * b) for _ in range(n2)]
e[0], e[-1] = round(b / 2) + 1, 0
x = [(2 * rng.random() - 1) * 2.0 ** ei for ei in e]
y = [(2 * rng.random() - 1) * 2.0 ** ei for ei in e]
for i in range(n2, n):
ei = round(b * (n - 1 - i) / (n - n2))
xi = (2 * rng.random() - 1) * 2.0 ** ei
yi = ((2 * rng.random() - 1) * 2.0 ** ei - float(exact(x, y))) / xi
x.append(xi)
y.append(yi)
return x, y
def test_error_free_transformations_are_exact() -> None:
rng = random.Random(1)
for _ in range(3000):
a = rng.uniform(-1, 1) * 2.0 ** rng.randint(-300, 300)
b = rng.uniform(-1, 1) * 2.0 ** rng.randint(-300, 300)
s, e = two_sum(a, b)
assert Fraction(s) + Fraction(e) == Fraction(a) + Fraction(b) and s == a + b
for tp in (two_product_fma, two_product_dekker):
p, f = tp(a, b)
assert Fraction(p) + Fraction(f) == Fraction(a) * Fraction(b) and p == a * b
def test_small_example() -> None:
x = [1e16, 1.0, -1e16]
y = [1.0, 1.0, 1.0]
assert naive_dot(x, y) == 0.0
assert dot2(x, y) == 1.0 and dot2(x, y, fma=False) == 1.0
def test_random_well_conditioned_is_correctly_rounded_or_near() -> None:
rng = random.Random(2)
for _ in range(300):
n = rng.randint(1, 60)
x = [rng.uniform(0, 1) for _ in range(n)]
y = [rng.uniform(0, 1) for _ in range(n)]
r, t = dot2(x, y), exact(x, y)
assert abs(Fraction(r) - t) <= bound(x, y)
assert abs(Fraction(r) - t) <= Fraction(math.ulp(r)) # nonnegative data: within one ulp
def test_ill_conditioned_bound_holds() -> None:
rng = random.Random(3)
worse = 0
for c in range(5, 40, 3):
for _ in range(4):
x, y = gen_dot(40, c, rng)
t = exact(x, y)
for f in (True, False):
r = dot2(x, y, fma=f)
assert abs(Fraction(r) - t) <= bound(x, y)
if abs(Fraction(naive_dot(x, y)) - t) > abs(Fraction(dot2(x, y)) - t):
worse += 1
assert worse > 20
def test_relative_error_like_twice_the_precision() -> None:
# Condition about 10^20: the naive result has no correct digit, Dot2 is still accurate.
rng = random.Random(4)
x, y = gen_dot(50, 20, rng)
t = exact(x, y)
cond = 2 * sum(abs(Fraction(a) * Fraction(b)) for a, b in zip(x, y)) / abs(t)
rel = abs(Fraction(dot2(x, y)) - t) / abs(t)
assert rel <= U + gamma(50) ** 2 * cond / 2 * (1 + U)
assert rel < Fraction(1, 10 ** 6)
def test_zeros_signs_and_empty() -> None:
assert dot2([], []) == 0.0
assert dot2([0.0, -0.0], [5.0, 3.0]) == 0.0
assert dot2([-2.0], [3.0]) == -6.0
big = 2.0 ** 1023
assert dot2([big, -big], [1.0, 1.0]) == 0.0
def test_nonfinite_and_underflow_limits() -> None:
# Outside the assumptions: an infinite product makes the correction NaN.
# (Python's math.fma raises ValueError on the invalid operation inf - inf.)
assert math.isnan(dot2([math.inf], [1.0], fma=False))
with pytest.raises(ValueError):
dot2([math.inf], [1.0])
# Underflow: TwoProduct is no longer exact, and the bound without an underflow term fails.
x = [2.0 ** -600, 1.0]
y = [2.0 ** -600, 0.0]
assert dot2(x, y) == 0.0 and exact(x, y) == Fraction(1, 2 ** 1200)
assert abs(Fraction(dot2(x, y)) - exact(x, y)) > bound(x, y)
def test_long_vectors_against_fsum_of_exact_parts() -> None:
rng = random.Random(5)
n = 2000
x = [rng.uniform(-1, 1) for _ in range(n)]
y = [rng.uniform(-1, 1) for _ in range(n)]
parts = [v for a, b in zip(x, y) for v in two_product_fma(a, b)]
ref = math.fsum(parts) # the correctly rounded dot product
assert abs(dot2(x, y) - ref) <= 2 * math.ulp(ref) + 1e-300
assert abs(Fraction(dot2(x, y)) - exact(x, y)) <= bound(x, y)
What the next level needs
- the tests did not run in this build