Goldschmidt division in Python
Iteratively scales numerator and denominator by correction factors until the denominator reaches one.
- File
goldschmidt_division.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/goldschmidt-division/. The authors are listed in prov/CRTBFS.jsonl.
"""Goldschmidt division of binary64 numbers in W-bit fixed point, with a final exact rounding step."""
import math
P = 53 # precision of binary64
W = 64 # fractional bits of the fixed-point datapath
TABLE_BITS = 7 # the seed table is indexed by the 7 leading fraction bits of the divisor
SEED_BITS = 10 # fractional bits of each table entry
ONE = 1 << W
def seed(big_b: int) -> int:
"""The table entry f0 ~ 1/B for a divisor B in [1, 2) given as an integer scaled by 2^W.
Entry i is 1/(1 + (i + 1/2)/2^7) rounded to SEED_BITS fractional bits, so that
|1 - B*f0| < 2^(-7). Returned scaled by 2^W. Cost: O(1).
"""
i = (big_b - ONE) >> (W - TABLE_BITS)
num = 1 << (SEED_BITS + TABLE_BITS + 1) # 2^(SEED_BITS) * 2^(TABLE_BITS + 1)
den = (1 << (TABLE_BITS + 1)) + 2 * i + 1 # 2^(TABLE_BITS + 1) * (1 + (i + 1/2)/2^7)
f = (2 * num + den) // (2 * den) # nearest integer to num/den
return f << (W - SEED_BITS)
def split(x: float) -> tuple[int, int]:
"""For a positive normal float x, the integer M with x = M * 2^(e - 52), 2^52 <= M < 2^53, and e. Cost: O(1)."""
m, e = math.frexp(x) # x = m * 2^e with 1/2 <= m < 1
return int(m * (1 << P)), e - 1
def goldschmidt_divide(a: float, b: float, iterations: int = 3, trace: list | None = None) -> float:
"""The correctly rounded (to nearest) binary64 quotient a / b by Goldschmidt's iteration.
Precondition: a and b are finite, nonzero and normal, the quotient is a normal
binary64 number (no underflow or overflow), and 3 <= iterations <= 64. The iteration runs
in W-bit fixed point with truncation; after 3 iterations the approximation is within
2^(-54) of A/B, and lines L8 and L9 fix the last bit with exact integer comparisons.
With fewer iterations the result can be wrong. Cost: O(iterations) multiplications
of (W+2)-bit integers, O(1) for fixed W.
If trace is a list, the tuples (f, n, d) of each step are appended to it.
"""
if not (math.isfinite(a) and math.isfinite(b)) or a == 0.0 or b == 0.0:
raise ValueError("a and b must be finite and nonzero")
sign = -1.0 if (a < 0.0) != (b < 0.0) else 1.0
ma, ea = split(abs(a)) # @line L1
mb, eb = split(abs(b))
big_a = ma << (W - P + 1) # A in [1, 2), scaled by 2^W
big_b = mb << (W - P + 1)
f = seed(big_b) # @line L2
n = (big_a * f) >> W # @line L3
d = (big_b * f) >> W
if trace is not None:
trace.append((f, n, d))
for _ in range(iterations): # @line L4
f = 2 * ONE - d # @line L5
n, d = (n * f) >> W, (d * f) >> W # @line L6
if trace is not None:
trace.append((f, n, d))
# The quotient A/B lies in (1/2, 2); its unit in the last place is 2^(-52) or 2^(-53).
shift = P - 1 if ma >= mb else P # unit = 2^(-shift)
q = (n + (1 << (W - shift - 1))) >> (W - shift) # @line L7
# Exact tests: A/B < (q - 1/2) * unit <=> 2 * ma * 2^shift < (2q - 1) * mb.
if 2 * ma * (1 << shift) < (2 * q - 1) * mb: # @line L8
q -= 1
if 2 * ma * (1 << shift) > (2 * q + 1) * mb: # @line L9
q += 1
return sign * math.ldexp(q, ea - eb - shift) # @line L10
test_goldschmidt_division.py 121 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/goldschmidt-division/. The authors are listed in prov/CRTBFS.jsonl.
"""Tests for Goldschmidt division: the result must equal the exactly rounded quotient."""
import math
import random
from fractions import Fraction
import pytest
from goldschmidt_division import ONE, W, goldschmidt_divide, seed
def rn(q: Fraction) -> float:
"""Round a positive or negative rational to the nearest binary64 (ties to even).
Fraction.__float__ divides the integers numerator and denominator, which CPython
rounds correctly; this is an oracle independent of the float division a / b.
"""
return q.numerator / q.denominator
def exact_quotient(a: float, b: float) -> float:
return rn(Fraction(a) / Fraction(b))
def test_examples() -> None:
assert goldschmidt_divide(1.0, 3.0) == 0.3333333333333333
assert goldschmidt_divide(1.0, 10.0) == 0.1
assert goldschmidt_divide(4195835.0, 3145727.0) == 1.333820449136241
assert goldschmidt_divide(3.0, 5.0) == 0.6
assert goldschmidt_divide(6.0, 3.0) == 2.0
assert goldschmidt_divide(1.0, 1.0) == 1.0
def test_signs() -> None:
assert goldschmidt_divide(-1.0, 3.0) == -exact_quotient(1.0, 3.0)
assert goldschmidt_divide(1.0, -3.0) == -exact_quotient(1.0, 3.0)
assert goldschmidt_divide(-7.0, -2.0) == 3.5
def test_rejects_zero_and_nonfinite() -> None:
for a, b in [(0.0, 1.0), (1.0, 0.0), (math.inf, 1.0), (1.0, math.nan)]:
with pytest.raises(ValueError):
goldschmidt_divide(a, b)
def test_seed_accuracy_proved_bound() -> None:
# |1 - B f0| < 2^(-7) for every table interval; check both ends and the middle of each.
for i in range(128):
for t in (0, 1 << 30, (1 << (W - 7)) - 1):
big_b = ONE + (i << (W - 7)) + t
f = seed(big_b)
eps = abs(Fraction(ONE * ONE - big_b * f, ONE * ONE))
assert eps < Fraction(1, 128)
def test_trace_converges_quadratically() -> None:
trace: list = []
goldschmidt_divide(1.0, 3.0, 3, trace)
errors = [abs(ONE - d) for _, _, d in trace]
assert errors[0] < ONE >> 7
for before, after in zip(errors, errors[1:]):
assert after <= (before * before >> W) + 1 # eta' <= eta^2 + 2^(-W)
def test_extremes_of_the_significands() -> None:
big = math.nextafter(2.0, 0.0) # 2 - 2^(-52)
cases = [(big, 1.0), (1.0, big), (big, big), (1.0, math.nextafter(1.0, 2.0)),
(math.nextafter(1.0, 2.0), 1.0), (1.5, 1.0), (1.0, 1.5)]
for a, b in cases:
assert goldschmidt_divide(a, b) == exact_quotient(a, b)
def test_large_and_small_normal_results() -> None:
tiny = 2.0 ** -1022
cases = [(1.0, 2.0 ** 1000), (2.0 ** 1000, 3.0), (tiny * 3.0, 3.0), (1e308, 1.0),
(math.ulp(0.0) * 2.0 ** 52 * 1.5, 1.5)]
for a, b in cases:
q = goldschmidt_divide(a, b)
assert q == exact_quotient(a, b)
def test_random_against_rational_oracle() -> None:
rng = random.Random(20240611)
for _ in range(20000):
a = rng.uniform(1.0, 2.0) * 2.0 ** rng.randint(-400, 400) * rng.choice((1.0, -1.0))
b = rng.uniform(1.0, 2.0) * 2.0 ** rng.randint(-400, 400) * rng.choice((1.0, -1.0))
assert goldschmidt_divide(a, b) == exact_quotient(a, b)
def test_random_bit_patterns_near_ties() -> None:
# Quotients of nearby significands and of integers stress the rounding decision.
rng = random.Random(7)
for _ in range(5000):
m = rng.randrange(1 << 52, 1 << 53)
k = rng.randrange(1, 1 << 20)
a = float(m)
b = float(m + k) if m + k < (1 << 53) else float(m - k)
assert goldschmidt_divide(a, b) == exact_quotient(a, b)
x, y = float(rng.randrange(1, 1 << 30)), float(rng.randrange(1, 1 << 30))
assert goldschmidt_divide(x, y) == exact_quotient(x, y)
def test_more_iterations_also_correct_and_fewer_can_fail() -> None:
rng = random.Random(3)
wrong_with_two = 0
for _ in range(2000):
a, b = rng.uniform(1.0, 2.0), rng.uniform(1.0, 2.0)
assert goldschmidt_divide(a, b, 4) == exact_quotient(a, b)
wrong_with_two += goldschmidt_divide(a, b, 2) != exact_quotient(a, b)
assert wrong_with_two > 0
def test_agrees_with_hardware_division() -> None:
rng = random.Random(11)
for _ in range(5000):
a, b = rng.uniform(-1e6, 1e6), rng.uniform(-1e6, 1e6)
if a != 0.0 and b != 0.0:
assert goldschmidt_divide(a, b) == a / b
What the next level needs
- the tests did not run in this build