Granlund-Montgomery division in Python
Division by a fixed divisor using a precomputed multiplicative inverse and shifts.
- File
granlund_montgomery_division.py, 84 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/granlund-montgomery-division/. The authors are listed in prov/N5CE12.jsonl.
"""Division by an invariant integer with a precomputed multiplier.
Words have w bits and beta = 2^w. Every operation below is one that a w-bit
machine provides: the high word of a w x w product, addition and subtraction
modulo beta, shifts and comparisons. Python ints model the words and every
intermediate result is checked to fit in one word.
"""
from dataclasses import dataclass
@dataclass(frozen=True)
class Divisor:
"""The precomputed data for dividing w-bit words by d."""
d: int
w: int
m: int
s1: int
s2: int
def _word(x: int, w: int) -> int:
if not 0 <= x < (1 << w):
raise AssertionError("intermediate value does not fit in a word")
return x
def precompute(d: int, w: int) -> Divisor:
"""Granlund-Montgomery data for 1 <= d < 2^w (their Figure 4.1)."""
if w < 1 or not 1 <= d < (1 << w):
raise ValueError("need 1 <= d < 2^w")
ell = (d - 1).bit_length() # @line L1
m = ((1 << w) * ((1 << ell) - d)) // d + 1 # @line L2
s1, s2 = min(ell, 1), max(ell - 1, 0) # @line L3
return Divisor(d, w, _word(m, w), s1, s2)
def divide(n: int, dv: Divisor) -> tuple[int, int]:
"""Return (n // d, n % d) for 0 <= n < 2^w with one multiply-high."""
w = dv.w
_word(n, w)
t = (dv.m * n) >> w # @line L4
if t > n:
raise AssertionError("t must not exceed n")
q = _word(t + ((n - t) >> dv.s1), w) >> dv.s2 # @line L5
return q, _word(n - q * dv.d, w) # @line L6
def reciprocal(d: int, w: int) -> int:
"""The Moller-Granlund reciprocal v = floor((beta^2 - 1) / d) - beta.
Precondition: d is normalized, 2^(w-1) <= d < 2^w.
"""
if w < 1 or not (1 << (w - 1)) <= d < (1 << w):
raise ValueError("d must be normalized")
return _word(((1 << (2 * w)) - 1) // d - (1 << w), w)
def div_2by1(u1: int, u0: int, d: int, v: int, w: int) -> tuple[int, int]:
"""Divide the two-word number (u1, u0) by a normalized d with u1 < d.
Moller and Granlund, Algorithm 4: one full product, one low product,
and at most two corrections. Returns (q, r) with u1 beta + u0 = q d + r.
"""
beta = 1 << w
mask = beta - 1
if not (0 <= u1 < d and 0 <= u0 < beta):
raise ValueError("need u1 < d and u0 < beta")
prod = v * u1 + (u1 << w) + u0 # @line L7
q1, q0 = prod >> w, prod & mask
_word(q1, w)
q1 = (q1 + 1) & mask # @line L8
r = (u0 - q1 * d) & mask
if r > q0: # @line L9
q1 = (q1 - 1) & mask
r = (r + d) & mask
if r >= d: # @line L10
q1 += 1
r -= d
return q1, r # @line L11
test_granlund_montgomery_division.py 71 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/granlund-montgomery-division/. The authors are listed in prov/N5CE12.jsonl.
import random
import pytest
from granlund_montgomery_division import div_2by1, divide, precompute, reciprocal
def test_example_w8_d7():
dv = precompute(7, 8)
assert (dv.m, dv.s1, dv.s2) == (37, 1, 2)
assert divide(200, dv) == (28, 4)
def test_exhaustive_w8():
for d in range(1, 256):
dv = precompute(d, 8)
for n in range(256):
assert divide(n, dv) == divmod(n, d)
def test_random_wide():
rng = random.Random(1994)
for w in (16, 32, 64, 128):
for _ in range(250):
d = rng.randint(1, (1 << w) - 1)
if rng.random() < 0.3:
d = (1 << rng.randint(0, w - 1)) + rng.randint(-1, 1)
d = min(max(d, 1), (1 << w) - 1)
dv = precompute(d, w)
for n in (rng.getrandbits(w), (1 << w) - 1, 0, d, d - 1):
assert divide(n, dv) == divmod(n, d)
def test_2by1_exhaustive_w6():
w = 6
for d in range(1 << (w - 1), 1 << w):
v = reciprocal(d, w)
for u1 in range(d):
for u0 in range(1 << w):
assert div_2by1(u1, u0, d, v, w) == divmod((u1 << w) + u0, d)
def test_2by1_example_w8():
v = reciprocal(200, 8)
assert v == 71
assert div_2by1(150, 77, 200, v, 8) == divmod(150 * 256 + 77, 200)
def test_2by1_random():
rng = random.Random(2011)
for w in (8, 16, 32, 64):
for _ in range(400):
d = rng.randint(1 << (w - 1), (1 << w) - 1)
v = reciprocal(d, w)
u1, u0 = rng.randint(0, d - 1), rng.getrandbits(w)
assert div_2by1(u1, u0, d, v, w) == divmod((u1 << w) + u0, d)
assert div_2by1(d - 1, (1 << w) - 1, d, v, w) == divmod(((d - 1) << w) + (1 << w) - 1, d)
def test_errors():
with pytest.raises(ValueError):
precompute(0, 8)
with pytest.raises(ValueError):
precompute(256, 8)
with pytest.raises(ValueError):
reciprocal(100, 8)
with pytest.raises(ValueError):
div_2by1(200, 0, 200, 71, 8)
What the next level needs
- the tests did not run in this build