Garner's algorithm in Python
Reconstructs the CRT solution in mixed-radix form using only small modular operations.
- File
garner_algorithm.py, 81 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/garner-algorithm/. The authors are listed in prov/PHSQZC.jsonl.
"""Garner's algorithm: Chinese remaindering through mixed-radix digits.
Given pairwise coprime moduli m_1, ..., m_t and residues v_i modulo m_i, it
finds the unique x with 0 <= x < M = m_1 ... m_t and x = v_i (mod m_i). The
digits u_i of x in the mixed radix (m_1, ..., m_t),
x = u_1 + u_2 m_1 + u_3 m_1 m_2 + ... + u_t m_1 ... m_(t-1),
are computed with operations modulo m_i only. The setup stores the t - 1
constants C_i = (m_1 ... m_(i-1))^(-1) mod m_i. Only the last step, which
evaluates x from its digits, works with numbers as large as M.
"""
from typing import List, NamedTuple, Sequence, Tuple
class Garner(NamedTuple):
m: Tuple[int, ...] # the moduli
c: Tuple[int, ...] # c[i] = (m[0] ... m[i-1])^(-1) mod m[i], c[0] = 0
def _inverse(a: int, n: int) -> Tuple[int, int]:
"""Return (g, u) with g = gcd(a, n) and u * a = g (mod n), by the
extended Euclidean algorithm with a loop."""
r0, r1 = n, a % n
u0, u1 = 0, 1
while r1 != 0:
t = r0 // r1
r0, r1 = r1, r0 - t * r1
u0, u1 = u1, u0 - t * u1
return r0, u0 % n
def garner_setup(moduli: Sequence[int]) -> Garner:
"""Precompute the constants C_i; raise ValueError if the moduli are not
positive and pairwise coprime."""
m = tuple(moduli)
if any(mi < 1 for mi in m):
raise ValueError("the moduli must be positive integers")
c = [0] * len(m)
for i in range(1, len(m)): # @line L1
p = 1
for j in range(i): # @line L2
p = p * (m[j] % m[i]) % m[i]
g, c[i] = _inverse(p, m[i]) # @line L3
if g != 1: # @line L4
raise ValueError(f"modulus {m[i]} is not coprime to the earlier moduli")
return Garner(m, tuple(c))
def garner_digits(residues: Sequence[int], ctx: Garner) -> List[int]:
"""Return the mixed-radix digits u_i, with 0 <= u_i < m_i."""
m, c = ctx.m, ctx.c
t = len(m)
if len(residues) != t:
raise ValueError("one residue per modulus is needed")
if t == 0:
return []
u = [0] * t
u[0] = residues[0] % m[0] # @line L5
for i in range(1, t): # @line L6
s = u[i - 1] % m[i] # @line L7
for j in range(i - 2, -1, -1): # @line L8
s = (s * m[j] + u[j]) % m[i] # @line L9
u[i] = (residues[i] - s) * c[i] % m[i] # @line L10
return u
def garner(residues: Sequence[int], ctx: Garner) -> int:
"""Return the unique x in [0, M) with x = residues[i] (mod m[i])."""
u = garner_digits(residues, ctx)
if not u:
return 0 # M = 1, and x = 0 is the only value
t = len(u)
x = u[t - 1] # @line L11
for i in range(t - 2, -1, -1): # @line L12
x = x * ctx.m[i] + u[i]
return x # @line L13
test_garner_algorithm.py 114 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/garner-algorithm/. The authors are listed in prov/PHSQZC.jsonl.
import math
import random
import pytest
from garner_algorithm import garner, garner_digits, garner_setup
def reference(residues, moduli):
"""The classical CRT sum with Python's pow, used only as the reference."""
big = math.prod(moduli)
x = 0
for v, mi in zip(residues, moduli):
w = big // mi
x += v * w * pow(w, -1, mi)
return x % big
def coprime_moduli(rng, pool, count):
chosen = []
for mi in rng.sample(pool, len(pool)):
if all(math.gcd(mi, mj) == 1 for mj in chosen):
chosen.append(mi)
if len(chosen) == count:
break
return chosen
POOL = [
1, 2, 3, 4, 5, 7, 8, 9, 11, 13, 16, 25, 27, 49, 64, 101, 121, 1 << 64, 15, 77, 221,
10**9 + 7, 998244353, 2**61 - 1, 2**127 - 1, 10**300 + 1, 3**700, 10**400 - 3, 2**521 - 1, 17**250,
]
def test_hac_example():
ctx = garner_setup([5, 7, 11, 13])
assert ctx.c == (0, 3, 6, 5)
assert garner_digits([2, 1, 3, 8], ctx) == [2, 4, 7, 5]
assert garner([2, 1, 3, 8], ctx) == 2192
def test_rsa_example():
ctx = garner_setup([61, 53])
assert ctx.c[1] == 20
assert garner([pow(2790, 53, 61), pow(2790, 49, 53)], ctx) == 65 == pow(2790, 2753, 3233)
def test_small_cases():
assert garner([], garner_setup([])) == 0
assert garner([0], garner_setup([1])) == 0
assert garner([0, 1], garner_setup([1, 2])) == 1
assert garner([1, 0], garner_setup([2, 1])) == 1
ctx = garner_setup([2, 3])
assert [garner([x % 2, x % 3], ctx) for x in range(6)] == list(range(6))
def test_every_value_small_systems():
for moduli in ([3, 4, 5], [8, 9, 7], [2, 1, 3], [16, 5, 3, 7], [11, 13, 17]):
ctx = garner_setup(moduli)
big = math.prod(moduli)
for x in range(big):
assert garner([x % mi for mi in moduli], ctx) == x
def test_random_against_reference():
rng = random.Random(5)
count = 0
for _ in range(600):
moduli = coprime_moduli(rng, POOL, rng.randrange(1, 9))
ctx = garner_setup(moduli)
big = math.prod(moduli)
for x in (0, 1 % big, big - 1, rng.randrange(big)):
res = [x % mi for mi in moduli]
got = garner(res, ctx)
assert got == x == reference(res, moduli)
digits = garner_digits(res, ctx)
assert all(0 <= d < mi for d, mi in zip(digits, moduli))
count += 1
assert count >= 2000
def test_random_residues_and_word_primes():
rng = random.Random(6)
primes = [p for p in range(10**6, 10**6 + 3000) if all(p % d for d in range(2, 1001))]
for _ in range(100):
moduli = rng.sample(primes, rng.randrange(1, 40))
ctx = garner_setup(moduli)
res = [rng.randrange(-10**20, 10**20) for _ in moduli]
x = garner(res, ctx)
assert 0 <= x < math.prod(moduli)
assert x == reference([v % mi for v, mi in zip(res, moduli)], moduli)
def test_not_coprime():
for moduli in ([4, 6], [3, 5, 9], [2, 2], [10**300 + 1, 3, 10**300 + 1], [0, 3], [5, -7]):
with pytest.raises(ValueError):
garner_setup(moduli)
rng = random.Random(7)
for _ in range(200):
a, b = rng.randrange(2, 10**6), rng.randrange(2, 10**6)
moduli = [a, b]
if math.gcd(a, b) == 1:
assert garner([1, 1], garner_setup(moduli)) == 1 % (a * b)
else:
with pytest.raises(ValueError):
garner_setup(moduli)
def test_length_mismatch():
with pytest.raises(ValueError):
garner([1], garner_setup([3, 5]))
What the next level needs
- the tests did not run in this build