Veltkamp splitting in Python
Splits a floating-point number into two halves with fewer significant bits, enabling exact products.
- File
veltkamp_splitting.py, 70 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/veltkamp-splitting/. The authors are listed in prov/Z6SW6X.jsonl.
"""Veltkamp splitting of a binary64 float into a high part and a low part.
With s = 27 the high part has at most 26 significant bits and the low part
at most 26, so products of parts are exact in binary64 (Dekker's product).
"""
from fractions import Fraction
P = 53
def veltkamp_split(a: float, s: int = 27) -> tuple[float, float]:
"""Return (a_hi, a_lo) with a = a_hi + a_lo exactly.
a_hi is a rounded to a multiple of 2^s ulp(a), so it has at most P - s
significant bits; |a_lo| <= 2^(s-1) ulp(a). Valid for 1 <= s <= P - 1,
finite a, and (2^s + 1) * a not overflowing. Cost: 1 multiplication and
3 additions.
"""
c = float((1 << s) + 1) # @line L1
gamma = c * a # @line L2
delta = a - gamma # @line L3
hi = gamma + delta # @line L4
lo = a - hi # @line L5
return hi, lo # @line L6
def significant_bits(x: float) -> int:
"""The number of bits from the leading to the trailing one bit of |x| (0 for 0)."""
if x == 0.0:
return 0
f = abs(Fraction(x))
n, d = f.numerator, f.denominator
# d is a power of two; strip trailing zero bits of n
while n % 2 == 0:
n //= 2
return n.bit_length()
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 veltkamp_split_emulated(a: Fraction, s: int, p: int, emin: int) -> tuple[Fraction, Fraction]:
"""Veltkamp splitting in a simulated binary format of precision p."""
def rn(v: Fraction) -> Fraction:
return round_to_precision(v, p, emin)
gamma = rn((2**s + 1) * a)
delta = rn(a - gamma)
hi = rn(gamma + delta)
return hi, rn(a - hi)
test_veltkamp_splitting.py 97 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/veltkamp-splitting/. The authors are listed in prov/Z6SW6X.jsonl.
import math
import random
from fractions import Fraction
from veltkamp_splitting import significant_bits, veltkamp_split, veltkamp_split_emulated
def ulp_exp(a: float) -> int:
"""The exponent of ulp(a) for a normal float a."""
return math.frexp(a)[1] - 53
def check(a: float, s: int) -> None:
hi, lo = veltkamp_split(a, s)
assert Fraction(hi) + Fraction(lo) == Fraction(a)
assert significant_bits(hi) <= 53 - s
assert significant_bits(lo) <= max(s - 1, 1)
if a != 0.0:
quantum = Fraction(2) ** (ulp_exp(a) + s)
assert (Fraction(hi) / quantum).denominator == 1
assert abs(Fraction(lo)) <= quantum / 2
def test_examples():
hi, lo = veltkamp_split(0.1)
assert hi == 0.09999999962747097 and lo == 3.7252903539730653e-10
hi, lo = veltkamp_split(math.pi)
assert (hi, lo) == (3.1415926814079285, -2.781813535079891e-08)
assert veltkamp_split(1.0) == (1.0, 0.0)
assert veltkamp_split(0.0) == (0.0, 0.0)
def test_random_s27():
rng = random.Random(2718)
for _ in range(10000):
a = math.ldexp(rng.uniform(-1, 1), rng.randint(-1000, 990))
check(a, 27)
def test_all_split_points():
rng = random.Random(31)
for s in range(1, 53):
for _ in range(300):
check(math.ldexp(rng.uniform(-1, 1), rng.randint(-900, 900)), s)
def test_hard_significands():
# all-ones significands and halfway patterns force rounding ties and carries
for s in (2, 26, 27, 28):
for a in (2.0 - 2.0**-52, 1.0 + 2.0**-52, 1.5, 1.0 + 2.0 ** (-53 + s), 1.0 + 2.0 ** (-54 + s)):
check(a, s)
check(-a * 2.0**100, s)
def test_products_of_parts_are_exact():
rng = random.Random(4)
for _ in range(5000):
a = rng.uniform(-1e5, 1e5)
b = rng.uniform(-1e5, 1e5)
ah, al = veltkamp_split(a)
bh, bl = veltkamp_split(b)
for x, y in ((ah, bh), (ah, bl), (al, bh), (al, bl)):
assert Fraction(x * y) == Fraction(x) * Fraction(y)
def test_exhaustive_small_precision():
for p in (4, 5, 6, 7):
emin = -4
for s in range(1, p):
for m in range(2 ** (p - 1), 2**p):
for e in (-4, 0, 3):
a = Fraction(m) * Fraction(2) ** (e - p + 1)
hi, lo = veltkamp_split_emulated(a, s, p, emin)
assert hi + lo == a
quantum = Fraction(2) ** (e - p + 1 + s)
assert (hi / quantum).denominator == 1, (p, s, a)
assert abs(lo) <= quantum / 2
def test_overflow_limit():
big = 1.7976931348623157e308
hi, lo = veltkamp_split(big)
assert math.isinf(hi) or math.isnan(hi) or math.isnan(lo)
hi, lo = veltkamp_split(big / 2.0**28)
assert Fraction(hi) + Fraction(lo) == Fraction(big / 2.0**28)
def test_subnormal_inputs():
rng = random.Random(5)
for _ in range(5000):
a = rng.uniform(-1, 1) * 2.0**-1022 * rng.choice([1, 2.0**-10, 2.0**-30, 2.0**-45])
hi, lo = veltkamp_split(a)
assert Fraction(hi) + Fraction(lo) == Fraction(a)
assert significant_bits(hi) <= 26 and significant_bits(lo) <= 26
What the next level needs
- the tests did not run in this build