Half-GCD in Python
A divide-and-conquer gcd for big integers and polynomials that runs in quasi-linear time.
- File
half_gcd.py, 116 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/half-gcd/. The authors are listed in prov/NZPF3Q.jsonl.
"""Half-GCD: a divide and conquer gcd for arbitrary precision integers.
``half_gcd(a, b)`` runs the Euclidean algorithm on a > b >= 0 up to the
first remainder below 2^m, with m = ceil(bitlength(a) / 2), and returns the
pair reached together with the 2 x 2 matrix of the quotients used. It finds
most quotients from the leading halves of the numbers, by two recursive calls
on numbers of about half the length. ``fast_gcd`` calls it in a loop.
Cost: O(M(n) log n) bit operations for n-bit inputs, where M(n) is the cost
of one n-bit multiplication. Python multiplies with Karatsuba's method, so
M(n) = O(n^1.59) here. The recursion depth is about log2(n / threshold).
"""
Matrix = tuple[int, int, int, int] # (m00, m01, m10, m11), row by row
IDENTITY: Matrix = (1, 0, 0, 1)
def _mul(m: Matrix, n: Matrix) -> Matrix:
a, b, c, d = m
e, f, g, h = n
return (a * e + b * g, a * f + b * h, c * e + d * g, c * f + d * h)
def _times_q(m: Matrix, q: int) -> Matrix:
"""Return m times the quotient matrix [[q, 1], [1, 0]]."""
a, b, c, d = m
return (a * q + b, a, c * q + d, c)
def _fixup(
r: Matrix, qs: list[int], a: int, b: int, limit: int
) -> tuple[Matrix, list[int], int, int]:
"""Apply r^-1 to (a, b), then undo trailing quotients until valid.
Unless every quotient is undone (then r is the identity and the result is
(a, b) itself), the result satisfies alpha > beta > 0 and alpha >= limit,
which makes it a pair of consecutive Euclidean remainders of (a, b).
"""
m00, m01, m10, m11 = r
det = m00 * m11 - m01 * m10 # +1 or -1
alpha = det * (m11 * a - m01 * b) # @line L11
beta = det * (m00 * b - m10 * a)
while qs and not (alpha > beta > 0 and alpha >= limit): # @line L12
q = qs.pop()
a0, a1, a2, a3 = r
r = (a1, a0 - q * a1, a3, a2 - q * a3) # r times [[0, 1], [1, -q]]
alpha, beta = q * alpha + beta, alpha
return r, qs, alpha, beta # @line L13
def half_gcd(
a: int, b: int, threshold: int = 64
) -> tuple[Matrix, list[int], int, int]:
"""Run Euclid on a > b >= 0 until the first remainder below 2^m.
Precondition: a > b >= 0 (a == b > 0 is also accepted), threshold >= 2.
Returns (M, qs, alpha, beta): qs are the Euclidean quotients used,
M = [[q1, 1], [1, 0]] ... [[qk, 1], [1, 0]], (a, b) = M (alpha, beta),
and alpha >= 2^m > beta with m = ceil(bitlength(a) / 2).
"""
if a < b or b < 0:
raise ValueError("half_gcd needs a >= b >= 0")
n = a.bit_length()
m = (n + 1) // 2 # @line L1
limit = 1 << m
mat: Matrix = IDENTITY
qs: list[int] = []
if b < limit: # @line L2
return mat, qs, a, b
if n <= threshold: # @line L3
while b >= limit:
q, rem = divmod(a, b)
a, b = b, rem
mat = _times_q(mat, q)
qs.append(q)
return mat, qs, a, b
r, rq, _, _ = half_gcd(a >> m, b >> m, threshold) # @line L4
mat, qs, a, b = _fixup(r, rq, a, b, limit) # @line L5
if b >= limit: # @line L6
q, rem = divmod(a, b)
a, b = b, rem
mat = _times_q(mat, q)
qs.append(q)
if b >= limit and 2 * (a.bit_length() - m) < n: # @line L7
k = 2 * m - a.bit_length()
s, sq, _, _ = half_gcd(a >> k, b >> k, threshold)
s, sq, a, b = _fixup(s, sq, a, b, limit) # @line L8
mat = _mul(mat, s)
qs.extend(sq)
while b >= limit: # @line L9
q, rem = divmod(a, b)
a, b = b, rem
mat = _times_q(mat, q)
qs.append(q)
return mat, qs, a, b # @line L10
def fast_gcd(a: int, b: int, threshold: int = 64) -> int:
"""Return gcd(a, b) >= 0 for any integers, with gcd(0, 0) = 0.
Signs are dropped first. Cost O(M(n) log n) bit operations.
"""
a, b = abs(a), abs(b) # @line L14
if a < b:
a, b = b, a
while b.bit_length() > threshold: # @line L15
_, _, a, b = half_gcd(a, b, threshold)
if b != 0:
a, b = b, a % b
while b != 0: # @line L16
a, b = b, a % b
return a # @line L17
test_half_gcd.py 139 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/half-gcd/. The authors are listed in prov/NZPF3Q.jsonl.
import math
import random
from half_gcd import fast_gcd, half_gcd
def brute_gcd(a: int, b: int) -> int:
"""Largest common divisor by direct search (small inputs only)."""
a, b = abs(a), abs(b)
if a == 0 or b == 0:
return a + b
return max(d for d in range(1, min(a, b) + 1) if a % d == 0 and b % d == 0)
def euclid_state(a: int, b: int) -> tuple[list[int], int, int]:
"""Reference: plain Euclid on a >= b >= 0 until the first remainder below 2^m."""
limit = 1 << ((a.bit_length() + 1) // 2)
qs = []
while b >= limit:
qs.append(a // b)
a, b = b, a % b
return qs, a, b
def fib(n: int) -> int:
x, y = 0, 1
for _ in range(n):
x, y = y, x + y
return x
def check_half(a: int, b: int, threshold: int) -> None:
mat, qs, alpha, beta = half_gcd(a, b, threshold)
assert (qs, alpha, beta) == euclid_state(a, b)
m00, m01, m10, m11 = mat
assert (m00 * alpha + m01 * beta, m10 * alpha + m11 * beta) == (a, b)
assert m00 * m11 - m01 * m10 == (-1) ** len(qs)
def test_example() -> None:
a, b = 3141592653, 2718281828
mat, qs, alpha, beta = half_gcd(a, b, threshold=16)
assert qs == [1, 6, 2, 2, 1, 2, 6, 8, 2, 1]
assert (alpha, beta) == (99065, 56321)
assert mat == (22873, 15548, 19791, 13453)
assert fast_gcd(a, b, threshold=16) == 1
# leading halves of the example and of the second call
assert half_gcd(47936, 41477, 16)[1:] == ([1, 6, 2, 2, 1, 2], 316, 65)
assert half_gcd(3353, 400, 16)[1:] == ([8, 2, 1], 94, 59)
def test_fixup_example() -> None:
# the leading 16 bits suggest a 7th quotient 2, the true one is 1
a, b = 2723649695, 1569894979
assert half_gcd(a >> 16, b >> 16, 16)[1] == [1, 1, 2, 1, 3, 2, 2]
check_half(a, b, 16)
for t in (2, 3, 4, 8):
check_half(a, b, t)
def test_zero_one_negative() -> None:
assert fast_gcd(0, 0) == 0
assert fast_gcd(0, 5) == 5
assert fast_gcd(-5, 0) == 5
assert fast_gcd(1, 0) == 1
assert fast_gcd(1, 1) == 1
assert fast_gcd(-12, 18) == 6
assert fast_gcd(12, -18) == 6
assert fast_gcd(-12, -18) == 6
big = 3**400
assert fast_gcd(-big, 0) == big
assert fast_gcd(0, -big) == big
def test_equal_inputs() -> None:
rng = random.Random(1)
for bits in (1, 10, 65, 300, 2000):
x = rng.getrandbits(bits) | 1
assert fast_gcd(x, x, 8) == x
assert fast_gcd(x, -x) == x
check_half(x, x, 4)
def test_powers_of_two() -> None:
for i in range(0, 700, 37):
for j in range(0, 700, 53):
assert fast_gcd(1 << i, 1 << j, 8) == 1 << min(i, j)
x, y = max(1 << i, 1 << j), min(1 << i, 1 << j)
check_half(x, y, 8)
def test_fibonacci() -> None:
for k in (10, 50, 200, 1000, 3000):
a, b = fib(k + 1), fib(k)
assert fast_gcd(a, b, 16) == 1
check_half(a, b, 16)
check_half(a, b, 64)
# gcd(F_m, F_n) = F_gcd(m, n)
assert fast_gcd(fib(1200), fib(900), 16) == fib(300)
def test_large_common_factor() -> None:
rng = random.Random(2)
for _ in range(40):
g = rng.getrandbits(rng.randint(1, 1500)) | 1
x = rng.getrandbits(rng.randint(1, 1500))
y = rng.getrandbits(rng.randint(1, 1500))
assert fast_gcd(g * x, g * y, 32) == g * math.gcd(x, y)
# Euclid reaches 0 while the remainder is still above 2^m
c = rng.randint(1, 60)
d = rng.randint(0, 60)
hi, lo = max(g * c, g * d), min(g * c, g * d)
check_half(hi, lo, 4)
def test_small_brute_force() -> None:
rng = random.Random(3)
for _ in range(600):
a = rng.randint(-300, 300)
b = rng.randint(-300, 300)
assert fast_gcd(a, b, rng.choice([2, 3, 4, 8])) == brute_gcd(a, b)
for a in range(0, 70):
for b in range(0, a + 1):
check_half(a, b, 2)
def test_random_large() -> None:
rng = random.Random(4)
for _ in range(600):
digits = rng.choice([1, 20, 100, 300, 600])
a = rng.randrange(10**digits) * rng.choice([1, -1])
b = rng.randrange(10 ** rng.randint(1, digits)) * rng.choice([1, -1])
threshold = rng.choice([2, 4, 16, 64])
assert fast_gcd(a, b, threshold) == math.gcd(a, b)
x, y = max(abs(a), abs(b)), min(abs(a), abs(b))
check_half(x, y, threshold)
What the next level needs
- the tests did not run in this build