Batch GCD in Python
Finds all nontrivial pairwise gcds of n integers via product and remainder trees in quasi-linear time.
- File
batch_gcd.py, 82 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/batch-gcd/. The authors are listed in prov/FT9662.jsonl.
"""Bernstein's batch GCD.
batch_gcd(values) returns, for each i, g_i = gcd(a_i, product of the a_j
with j != i), using a product tree and a remainder tree. Then g_i > 1 exactly
when a_i shares a factor with another input. Signs are ignored. Zero is
allowed: gcd(0, x) = |x|, so a single zero gets the product of the others,
and every other entry a_j gets |a_j|. The empty product is 1.
shared_factor_pairs finishes the job and lists the pairs with a common factor.
"""
def euclid_gcd(a: int, b: int) -> int:
"""gcd(a, b) >= 0 by the Euclidean algorithm."""
a, b = abs(a), abs(b)
while b:
a, b = b, a % b
return a
def product_tree(values: list[int]) -> list[list[int]]:
"""Levels of the product tree, leaves first; the last level is [product]."""
tree = [list(values)] # @line L1
h = 0 # @line L2
while len(tree[h]) > 1: # @line L3
level = tree[h]
tree.append( # @line L4
[level[j] * level[j + 1] for j in range(0, len(level) - 1, 2)]
+ ([level[-1]] if len(level) % 2 else [])
)
h += 1 # @line L5
return tree
def _batch_positive(values: list[int]) -> list[int]:
tree = product_tree(values)
h = len(tree) - 1
rem = tree[h] # @line L6
for level in range(h - 1, -1, -1): # @line L7
nodes = tree[level]
rem = [rem[j // 2] % (nodes[j] * nodes[j]) for j in range(len(nodes))] # @line L8
g = []
for i, a in enumerate(values): # @line L9
g.append(euclid_gcd(a, rem[i] // a)) # @line L10
return g # @line L11
def batch_gcd(values: list[int]) -> list[int]:
"""[gcd(a_i, prod_{j != i} a_j) for each i]."""
a = [abs(x) for x in values]
zeros = [i for i, x in enumerate(a) if x == 0]
if len(zeros) >= 2:
return [0 if x == 0 else x for x in a]
if len(zeros) == 1:
rest = 1
for x in a:
if x != 0:
rest *= x
return [rest if x == 0 else x for x in a]
if not a:
return []
return _batch_positive(a)
def shared_factor_pairs(values: list[int]) -> list[tuple[int, int, int]]:
"""All (i, j, gcd(a_i, a_j)) with i < j and gcd(a_i, a_j) > 1.
Batch GCD flags the inputs that share a factor with some other input;
pairwise gcds are then taken among the flagged inputs only.
"""
g = batch_gcd(values)
flagged = [i for i in range(len(values)) if g[i] != 1]
pairs = []
for x in range(len(flagged)):
for y in range(x + 1, len(flagged)):
i, j = flagged[x], flagged[y]
d = euclid_gcd(values[i], values[j])
if d != 1:
pairs.append((i, j, d))
return pairs
test_batch_gcd.py 123 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/batch-gcd/. The authors are listed in prov/FT9662.jsonl.
import math
import random
from batch_gcd import batch_gcd, euclid_gcd, product_tree, shared_factor_pairs
def fib(n: int) -> int:
a, b = 0, 1
for _ in range(n):
a, b = b, a + b
return a
def reference(values: list[int]) -> list[int]:
"""gcd(a_i, product of the others) with math.gcd (reference only)."""
out = []
for i, a in enumerate(values):
rest = 1
for j, b in enumerate(values):
if j != i:
rest *= b
out.append(math.gcd(a, rest))
return out
def reference_pairs(values: list[int]) -> list[tuple[int, int, int]]:
out = []
for i in range(len(values)):
for j in range(i + 1, len(values)):
d = math.gcd(values[i], values[j])
if d != 1:
out.append((i, j, d))
return out
def test_example() -> None:
assert batch_gcd([12, 18, 35, 49, 10]) == [12, 6, 35, 7, 10]
moduli = [101 * 103, 101 * 107, 109 * 113, 113 * 127, 131 * 137]
assert batch_gcd(moduli) == [101, 101, 113, 113, 1]
assert shared_factor_pairs(moduli) == [(0, 1, 101), (2, 3, 113)]
def test_product_tree_shape() -> None:
tree = product_tree([2, 3, 5, 7, 11])
assert tree == [[2, 3, 5, 7, 11], [6, 35, 11], [210, 11], [2310]]
def test_edge_inputs() -> None:
assert batch_gcd([]) == []
assert batch_gcd([7]) == [1]
assert batch_gcd([0]) == [1]
assert batch_gcd([0, 6, -4]) == [24, 6, 4]
assert batch_gcd([0, 0, 5]) == [0, 0, 5]
assert batch_gcd([1, 1, 1]) == [1, 1, 1]
assert batch_gcd([-6, 10, 1]) == [2, 2, 1]
assert batch_gcd([9, 9]) == [9, 9]
for values in ([0], [0, 6, -4], [0, 0, 5], [-6, 10, 1], [9, 9], [1, 1, 1]):
assert batch_gcd(values) == reference(values)
def test_euclid_gcd() -> None:
rng = random.Random(3)
for _ in range(500):
a = rng.randrange(-(10**rng.randrange(1, 300)), 10**rng.randrange(1, 300))
b = rng.choice([0, 1, -1, a, 2**rng.randrange(200), rng.randrange(-(10**50), 10**50)])
assert euclid_gcd(a, b) == math.gcd(a, b)
for k in range(2, 400):
assert euclid_gcd(fib(k + 1), fib(k)) == 1
def test_random_small_lists() -> None:
rng = random.Random(2012)
for _ in range(500):
n = rng.randrange(1, 12)
values = [
rng.choice(
[
0,
1,
-1,
rng.randrange(-200, 200),
2**rng.randrange(0, 40),
fib(rng.randrange(1, 60)),
rng.randrange(1, 10**6),
]
)
for _ in range(n)
]
if rng.random() < 0.2 and values:
values.append(values[0]) # equal inputs
assert batch_gcd(values) == reference(values)
if 0 not in values:
assert shared_factor_pairs(values) == reference_pairs(values)
def test_consecutive_fibonacci_are_coprime() -> None:
values = [fib(k) for k in range(3, 60, 2)]
# Odd-index Fibonacci numbers: gcd(F_i, F_j) = F_gcd(i, j).
assert batch_gcd(values) == reference(values)
values = [fib(k) for k in (101, 102, 103)]
assert batch_gcd(values) == [1, 1, 1]
def test_large_rsa_like_moduli() -> None:
rng = random.Random(77)
for _ in range(25):
factors = [rng.randrange(10**149, 10**150) | 1 for _ in range(12)]
values = []
for _ in range(10):
p, q = rng.sample(factors, 2)
values.append(p * q)
values.append(rng.randrange(10**299, 10**300))
assert batch_gcd(values) == reference(values)
assert shared_factor_pairs(values) == reference_pairs(values)
def test_many_inputs() -> None:
rng = random.Random(5)
values = [rng.randrange(2, 10**30) for _ in range(300)]
assert batch_gcd(values) == reference(values)
What the next level needs
- the tests did not run in this build