Brown's modular GCD algorithm in Python
Computes integer polynomial gcds by gcds modulo several primes followed by Chinese remaindering.
- File
brown_modular_gcd.py, 185 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/brown-modular-gcd/. The authors are listed in prov/02TR5G.jsonl.
"""Brown's modular gcd of two polynomials in Z[x], with the small primes method.
A polynomial a_0 + a_1 x + ... + a_n x^n is the list [a_0, ..., a_n] of
Python ints, lowest degree first; the zero polynomial is []. The gcd is
returned with a positive leading coefficient, and gcd(0, 0) = [].
For each prime p that does not divide b = gcd(lc f, lc g), the algorithm
computes the monic gcd v of f and g modulo p and the cofactors f / v and
g / v modulo p. Images of the smallest degree are combined by the Chinese
remainder theorem. A norm test, not a trial division, proves that the
combined polynomial is the gcd. All loops are iterative.
"""
from __future__ import annotations
from math import gcd
from typing import Iterable, Iterator, List, Optional, Sequence, Tuple
Poly = List[int]
def trim(a: Sequence[int]) -> Poly:
"""Return a copy of a without trailing zero coefficients."""
out = list(a)
while out and out[-1] == 0:
out.pop()
return out
def content(a: Sequence[int]) -> int:
"""The gcd of the coefficients, 0 for the zero polynomial."""
c = 0
for x in a:
c = gcd(c, x)
return c
def primitive_part(a: Sequence[int]) -> Poly:
"""a divided by its content, with a positive leading coefficient."""
a = trim(a)
if not a:
return []
c = content(a)
if a[-1] < 0:
c = -c
return [x // c for x in a]
def _mod_p(a: Sequence[int], p: int) -> Poly:
return trim([x % p for x in a])
def _divmod_p(a: Sequence[int], b: Sequence[int], p: int) -> Tuple[Poly, Poly]:
"""Quotient and remainder of a by the nonzero b in (Z/pZ)[x]."""
r = _mod_p(a, p)
inv = pow(b[-1], -1, p)
db = len(b) - 1
if len(r) <= db:
return [], r
q = [0] * (len(r) - db)
for i in range(len(r) - 1, db - 1, -1):
c = r[i] * inv % p
if c:
q[i - db] = c
for j in range(db + 1):
r[i - db + j] = (r[i - db + j] - c * b[j]) % p
return trim(q), trim(r[:db])
def gcd_mod_p(a: Sequence[int], b: Sequence[int], p: int) -> Poly:
"""The monic gcd of a and b in (Z/pZ)[x], by the Euclidean algorithm."""
a, b = _mod_p(a, p), _mod_p(b, p)
while b:
a, b = b, _divmod_p(a, b, p)[1]
if not a:
return []
inv = pow(a[-1], -1, p)
return [x * inv % p for x in a]
def _symmetric(a: Sequence[int], m: int) -> Poly:
return trim([x - m if 2 * x > m else x for x in a])
def _norm1(a: Sequence[int]) -> int:
return sum(abs(x) for x in a)
def _crt(old: Sequence[int], m: int, new: Sequence[int], p: int) -> Poly:
"""Combine residues modulo m and modulo p into residues modulo m p (0 .. m p - 1)."""
n = max(len(old), len(new))
minv = pow(m, -1, p)
out = []
for i in range(n):
x = old[i] if i < len(old) else 0
y = new[i] if i < len(new) else 0
out.append(x + m * ((y - x) * minv % p))
return out
def _is_prime(n: int) -> bool:
if n < 2:
return False
for q in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37):
if n % q == 0:
return n == q
d, s = n - 1, 0
while d % 2 == 0:
d, s = d // 2, s + 1
for a in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37):
x = pow(a, d, n)
if x in (1, n - 1):
continue
for _ in range(s - 1):
x = x * x % n
if x == n - 1:
break
else:
return False
return True
def word_primes(start: int = 2 ** 31) -> Iterator[int]:
"""Primes below start, in decreasing order."""
n = start - 1
while n > 2:
if _is_prime(n):
yield n
n -= 1
def brown_gcd(f: Sequence[int], g: Sequence[int], primes: Optional[Iterable[int]] = None,
log: Optional[list] = None) -> Poly:
"""Return gcd(f, g) in Z[x] with a positive leading coefficient.
primes is an iterable of distinct primes to use (default: primes below
2^31). If log is a list, one entry (p, degree, action) is appended for
every prime tried.
"""
f, g = trim(f), trim(g)
if not f or not g:
h = f or g
return [-x for x in h] if h and h[-1] < 0 else h
c = gcd(content(f), content(g)) # @line L1
f, g = primitive_part(f), primitive_part(g) # @line L1
if len(f) == 1 or len(g) == 1:
return [c]
b = gcd(f[-1], g[-1]) # @line L2
d = min(len(f), len(g)) - 1 # @line L2
m = 1 # @line L2
bf = [b * x for x in f]
bg = [b * x for x in g]
big = 2 * max(max(abs(x) for x in bf), max(abs(x) for x in bg))
H: Poly = []
U: Poly = []
W: Poly = []
for p in (primes if primes is not None else word_primes()): # @line L3
if b % p == 0: # @line L3
continue
v = gcd_mod_p(f, g, p) # @line L4
u = _divmod_p(f, v, p)[0] # @line L4
w = _divmod_p(g, v, p)[0] # @line L4
v = [b * x % p for x in v] # @line L4
if len(v) - 1 > d: # @line L5
if log is not None:
log.append((p, len(v) - 1, "skipped"))
continue
if len(v) - 1 < d or m == 1: # @line L6
H, U, W, m, d = v, u, w, p, len(v) - 1 # @line L6
step = "restarted"
else:
step = "combined"
H, U, W = _crt(H, m, v, p), _crt(U, m, u, p), _crt(W, m, w, p) # @line L7
m = m * p # @line L7
hs, us, ws = _symmetric(H, m), _symmetric(U, m), _symmetric(W, m) # @line L8
n1 = _norm1(hs)
ok = 2 * n1 * _norm1(us) < m and 2 * n1 * _norm1(ws) < m and big < m # @line L9
if log is not None:
log.append((p, d, "accepted" if ok else step))
if ok:
return [c * x for x in primitive_part(hs)] # @line L10
raise ValueError("ran out of primes")
test_brown_modular_gcd.py 239 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/brown-modular-gcd/. The authors are listed in prov/02TR5G.jsonl.
import random
from fractions import Fraction
from math import gcd
import pytest
from brown_modular_gcd import brown_gcd, gcd_mod_p, word_primes
class ModP:
"""An element of Z/pZ, enough for the tests."""
def __init__(self, v: int, p: int) -> None:
self.p = p
self.v = v % p
def _c(self, o):
return o.v if isinstance(o, ModP) else o
def __add__(self, o):
return ModP(self.v + self._c(o), self.p)
__radd__ = __add__
def __sub__(self, o):
return ModP(self.v - self._c(o), self.p)
def __rsub__(self, o):
return ModP(self._c(o) - self.v, self.p)
def __mul__(self, o):
return ModP(self.v * self._c(o), self.p)
__rmul__ = __mul__
def __neg__(self):
return ModP(-self.v, self.p)
def __truediv__(self, o):
d = self._c(o) % self.p
if d == 0:
raise ZeroDivisionError
return ModP(self.v * pow(d, -1, self.p), self.p)
def __rtruediv__(self, o):
return ModP(o, self.p) / self
def __eq__(self, o):
return self.v == self._c(o) % self.p
def __hash__(self):
return hash(self.v)
def __repr__(self):
return f"{self.v} mod {self.p}"
# Independent oracle: the Euclidean algorithm over a field (Fraction or ModP).
def ref_trim(a):
a = list(a)
while a and a[-1] == 0:
a.pop()
return a
def ref_rem(a, b):
a = ref_trim(a)
while len(a) >= len(b):
c = a[-1] / b[-1]
s = len(a) - len(b)
a = ref_trim([a[i] - (c * b[i - s] if i >= s else 0) for i in range(len(a))])
return a
def ref_monic_gcd(a, b):
a, b = ref_trim(a), ref_trim(b)
while b:
a, b = b, ref_rem(a, b)
if not a:
return []
lead = a[-1]
return [x / lead for x in a]
def ref_mul(a, b):
if not a or not b:
return []
out = [0] * (len(a) + len(b) - 1)
for i in range(len(a)):
for j in range(len(b)):
out[i + j] += a[i] * b[j]
return ref_trim(out)
def ref_int_gcd(f, g):
"""gcd in Z[x] = gcd of contents times the primitive part of the gcd over Q."""
f, g = ref_trim(f), ref_trim(g)
if not f or not g:
h = f or g
return [-x for x in h] if h and h[-1] < 0 else h
cf = cg = 0
for x in f:
cf = gcd(cf, x)
for x in g:
cg = gcd(cg, x)
q = ref_monic_gcd([Fraction(x) for x in f], [Fraction(x) for x in g])
den = 1
for x in q:
den = den * x.denominator // gcd(den, x.denominator)
ints = [int(x * den) for x in q]
c = 0
for x in ints:
c = gcd(c, x)
return [gcd(cf, cg) * x // c for x in ints]
def small_primes():
n = 2
while True:
if all(n % q for q in range(2, int(n ** 0.5) + 1)):
yield n
n += 1
def rand_poly(rng, deg, bound, monic=False):
a = [rng.randint(-bound, bound) for _ in range(deg + 1)]
if monic:
a[-1] = 1
while a[-1] == 0:
a[-1] = rng.randint(-bound, bound)
return a
def test_worked_example():
log = []
f = [3, 1, 6, 2]
g = [-1, 3, -2, 6]
assert brown_gcd(f, g, [5, 7, 11, 13], log) == [1, 0, 2]
assert log == [(5, 3, "restarted"), (7, 2, "restarted"), (11, 2, "accepted")]
assert ref_int_gcd(f, g) == [1, 0, 2]
@pytest.mark.parametrize("which", ["word", "small", "ntt"])
def test_random_against_rational_euclid(which):
rng = random.Random({"word": 1, "small": 2, "ntt": 3}[which])
for _ in range(60):
G = rand_poly(rng, rng.randint(0, 4), 9, monic=rng.random() < 0.3)
A = rand_poly(rng, rng.randint(0, 5), 9)
B = rand_poly(rng, rng.randint(0, 5), 9)
f = ref_mul([rng.choice([1, 2, -3, 6])], ref_mul(G, A))
g = ref_mul([rng.choice([1, 4, -9])], ref_mul(G, B))
if which == "word":
primes = word_primes()
elif which == "small":
primes = small_primes()
else:
primes = iter([998244353, 3, 7, 2] + list(_take(word_primes(2 ** 20), 40)))
assert brown_gcd(f, g, primes) == ref_int_gcd(f, g)
def _take(it, k):
out = []
for x in it:
out.append(x)
if len(out) == k:
break
return out
def test_unlucky_primes_are_skipped():
# x + 3 and x - 2 agree modulo 5, so 5 is unlucky for this pair
f = ref_mul([1, 0, 1], [3, 1])
g = ref_mul([1, 0, 1], [-2, 1])
log = []
assert brown_gcd(f, g, [5, 2, 3, 7, 11, 13, 17], log) == [1, 0, 1]
assert log[0] == (5, 3, "restarted")
# many primes in a row are unlucky: x and x + 30 agree modulo 2, 3 and 5
f = ref_mul([1, 1, 1], [0, 1])
g = ref_mul([1, 1, 1], [30, 1])
log = []
assert brown_gcd(f, g, small_primes(), log) == [1, 1, 1]
assert [e[1] for e in log[:3]] == [3, 3, 3]
@pytest.mark.parametrize("p", [2, 3, 7, 998244353])
def test_gcd_mod_p(p):
rng = random.Random(p % 1000)
for _ in range(80):
a = [rng.randrange(p) for _ in range(rng.randint(0, 7))]
b = [rng.randrange(p) for _ in range(rng.randint(0, 7))]
if rng.random() < 0.5:
c = [rng.randrange(p) for _ in range(rng.randint(1, 3))]
a, b = ref_mul(a, c), ref_mul(b, c)
want = ref_monic_gcd([ModP(x, p) for x in a], [ModP(x, p) for x in b])
assert gcd_mod_p(a, b, p) == [x.v for x in want]
def test_edge_cases():
assert brown_gcd([], []) == []
assert brown_gcd([0, 0], [0]) == []
assert brown_gcd([], [-2, -4]) == [2, 4]
assert brown_gcd([3, 6], []) == [3, 6]
assert brown_gcd([6], [4]) == [2]
assert brown_gcd([6, 12], [4, 0, 8]) == [2]
assert brown_gcd([5], [1, 1]) == [1]
# leading zero coefficients are stripped
assert brown_gcd([1, 0, -1, 0, 0], [1, 1, 0]) == [1, 1]
# coprime, equal, and dividing inputs
assert brown_gcd([1, 1], [-1, 1]) == [1]
assert brown_gcd([-2, 0, -4], [-2, 0, -4]) == [2, 0, 4]
assert brown_gcd([-1, 0, 1], [1, 1]) == [1, 1]
# negative leading coefficients give a positive gcd
assert brown_gcd([1, -1], [-1, 0, 1]) == [-1, 1]
# contents: gcd(6 (x + 1), 4 (x + 1) (x - 1)) = 2 (x + 1)
assert brown_gcd([6, 6], [-4, 0, 4]) == [2, 2]
# very different degrees
big = ref_mul([1, 1], [1] * 30)
assert brown_gcd(big, [3, 3]) == [1, 1]
# degree 1 non-monic gcd
assert brown_gcd(ref_mul([1, 3], [5, 7]), ref_mul([1, 3], [-2, 0, 9])) == [1, 3]
def test_large_coefficients_need_many_primes():
rng = random.Random(9)
for _ in range(10):
G = rand_poly(rng, 3, 10 ** 30)
f = ref_mul(G, rand_poly(rng, 3, 10 ** 20))
g = ref_mul(G, rand_poly(rng, 2, 10 ** 20))
log = []
assert brown_gcd(f, g, small_primes(), log) == ref_int_gcd(f, g)
assert len(log) > 10
def test_runs_out_of_primes():
with pytest.raises(ValueError):
brown_gcd([1, 0, 1], [1, 0, 1], [2])
What the next level needs
- the tests did not run in this build