Heuristic GCD in Python
Computes integer polynomial GCDs by evaluating at a large integer and reading the GCD from the result.
- File
heuristic_gcd.py, 122 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/heuristic-gcd/. The authors are listed in prov/N1KXF6.jsonl.
"""Heuristic gcd of two polynomials in Z[x] (Char, Geddes and Gonnet).
A polynomial is a list of integers [a_0, a_1, ..., a_d], lowest degree first,
with a nonzero last entry. The zero polynomial is []. The gcd is returned with
a positive leading coefficient.
"""
from math import gcd
from typing import List, Optional, Tuple
Poly = List[int]
def trim(f: Poly) -> Poly:
"""Remove zero coefficients at the top."""
f = list(f)
while f and f[-1] == 0:
f.pop()
return f
def content(f: Poly) -> int:
"""The gcd of the coefficients, 0 for the zero polynomial."""
c = 0
for a in f:
c = gcd(c, a)
return c
def primitive_part(f: Poly) -> Poly:
"""f divided by its content, with a positive leading coefficient."""
c = content(f)
if c == 0:
return []
if f[-1] < 0:
c = -c
return [a // c for a in f]
def norm(f: Poly) -> int:
"""The max norm |f|_inf."""
return max((abs(a) for a in f), default=0)
def evaluate(f: Poly, xi: int) -> int:
"""f(xi) by Horner's rule."""
v = 0
for a in reversed(f):
v = v * xi + a
return v
def genpoly(gamma: int, xi: int) -> Poly:
"""The polynomial whose coefficients are the symmetric xi-adic digits of gamma.
Each digit lies in (-xi/2, xi/2], and the result h satisfies h(xi) = gamma.
"""
h: Poly = []
while gamma != 0:
r = gamma % xi
if 2 * r > xi:
r -= xi
h.append(r)
gamma = (gamma - r) // xi
return h
def divides(h: Poly, f: Poly) -> bool:
"""True when h divides f in Z[x] (h nonzero), by exact long division."""
r = list(f)
dh = len(h) - 1
lc = h[-1]
while len(r) - 1 >= dh and r:
q, m = divmod(r[-1], lc)
if m != 0:
return False
shift = len(r) - 1 - dh
for i, a in enumerate(h):
r[shift + i] -= q * a
r = trim(r)
return not r
def heuristic_gcd(f: Poly, g: Poly, tries: int = 6) -> Optional[Poly]:
"""Return gcd(f, g) in Z[x], or None if every evaluation point fails.
None means: fall back to another method, such as a modular gcd.
"""
f, g = trim(f), trim(g)
if not f or not g:
h = f or g
return [-a for a 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)
xi = 2 * min(norm(f), norm(g)) + 2 # @line L2
for _ in range(tries): # @line L3
gamma = gcd(evaluate(f, xi), evaluate(g, xi)) # @line L4
h = genpoly(gamma, xi) # @line L5
h = primitive_part(h) # @line L6
if divides(h, f) and divides(h, g): # @line L7
return [c * a for a in h]
xi = xi * 73794 // 27011 # @line L8
return None # @line L9
def heuristic_gcd_trace(f: Poly, g: Poly, tries: int = 6) -> List[Tuple[int, int, Poly, bool]]:
"""The list of (xi, gamma, h, accepted) for each evaluation point tried."""
f, g = primitive_part(trim(f)), primitive_part(trim(g))
xi = 2 * min(norm(f), norm(g)) + 2
out = []
for _ in range(tries):
gamma = gcd(evaluate(f, xi), evaluate(g, xi))
h = primitive_part(genpoly(gamma, xi))
ok = divides(h, f) and divides(h, g)
out.append((xi, gamma, h, ok))
if ok:
break
xi = xi * 73794 // 27011
return out
test_heuristic_gcd.py 122 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/heuristic-gcd/. The authors are listed in prov/N1KXF6.jsonl.
import random
from fractions import Fraction
from math import gcd, lcm
from heuristic_gcd import (
content,
divides,
evaluate,
genpoly,
heuristic_gcd,
heuristic_gcd_trace,
primitive_part,
trim,
)
def mul(a, b):
if not a or not b:
return []
r = [0] * (len(a) + len(b) - 1)
for i, x in enumerate(a):
for j, y in enumerate(b):
r[i + j] += x * y
return r
def oracle_gcd(f, g):
"""gcd in Z[x] from the Euclidean algorithm over Q, then made primitive."""
f, g = trim(f), trim(g)
if not f or not g:
h = f or g
return [-a for a in h] if h and h[-1] < 0 else h
a = [Fraction(x) for x in f]
b = [Fraction(x) for x in g]
while b:
r = list(a)
while len(r) >= len(b) and r:
q = r[-1] / b[-1]
s = len(r) - len(b)
for i, y in enumerate(b):
r[s + i] -= q * y
while r and r[-1] == 0:
r.pop()
a, b = b, r
den = 1
for x in a:
den = lcm(den, x.denominator)
h = primitive_part([int(x * den) for x in a])
c = gcd(content(f), content(g))
return [c * x for x in h]
def rand_poly(rng, deg, size):
p = [rng.randint(-size, size) for _ in range(deg + 1)]
if p[-1] == 0:
p[-1] = rng.choice([-1, 1]) * rng.randint(1, size)
return p
def test_worked_example():
f, g = [0, 2, 1], [-4, 0, 1] # x^2 + 2x and x^2 - 4
tr = heuristic_gcd_trace(f, g)
assert tr == [(6, 16, [-2, 3], False), (16, 36, [2, 1], True)]
assert heuristic_gcd(f, g) == [2, 1]
def test_coprime_example():
tr = heuristic_gcd_trace([1, -1, 1], [1, 3])
assert tr == [(4, 13, [1, -1, 1], False), (10, 1, [1], True)]
assert heuristic_gcd([1, -1, 1], [1, 3]) == [1]
def test_genpoly_inverts_evaluation():
rng = random.Random(11)
for _ in range(500):
xi = rng.randint(2, 50)
h = [rng.randint(-((xi - 1) // 2), xi // 2) for _ in range(rng.randint(1, 6))]
h = trim(h)
assert genpoly(evaluate(h, xi), xi) == h
def test_random_against_euclid():
rng = random.Random(2024)
answered = 0
for _ in range(600):
h = rand_poly(rng, rng.randint(0, 3), 6)
a = rand_poly(rng, rng.randint(0, 3), 6)
b = rand_poly(rng, rng.randint(0, 3), 6)
f, g = mul(h, a), mul(h, b)
want = oracle_gcd(f, g)
got = heuristic_gcd(f, g)
if got is not None:
answered += 1
assert got == want
assert divides(primitive_part(got), primitive_part(f))
assert answered >= 590
def test_large_coefficients():
rng = random.Random(7)
for _ in range(100):
h = rand_poly(rng, 4, 10 ** 6)
f = mul(h, rand_poly(rng, 3, 10 ** 6))
g = mul(h, rand_poly(rng, 5, 10 ** 6))
got = heuristic_gcd(f, g)
assert got is None or got == oracle_gcd(f, g)
def test_edge_cases():
assert heuristic_gcd([], []) == []
assert heuristic_gcd([], [0, -2]) == [0, 2]
assert heuristic_gcd([3, 6], []) == [3, 6]
assert heuristic_gcd([6], [4]) == [2]
assert heuristic_gcd([6], [1, 1]) == [1]
assert heuristic_gcd([-2, 0, -2], [0, 0, 4]) == [2]
f = [-1, 0, 0, 0, 0, 0, 1] # x^6 - 1
assert heuristic_gcd(f, f) == f
assert heuristic_gcd([1, -1, 0, 0, 1], [-1, 1, 0, 0, -1]) == [1, -1, 0, 0, 1]
assert heuristic_gcd(mul([1, 1], [2, 3]), mul([1, 1], [5, -7])) == [1, 1]
What the next level needs
- the tests did not run in this build