Lehmer's GCD in Python
Speeds up Euclid on big integers by running single-precision steps on the leading digits.
- File
lehmers_gcd.py, 44 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/lehmers-gcd/. The authors are listed in prov/RN6MV6.jsonl.
def lehmer_gcd(a: int, b: int, w: int = 64) -> int:
"""Return gcd(a, b) >= 0 by Lehmer's algorithm, with gcd(0, 0) == 0.
Precondition: a and b are any Python integers (negative values and 0 are
allowed; signs are dropped first), and w >= 2 is the word size in bits.
Numbers below 2**w count as single precision.
The outer loop works on the leading w bits of u and v. It runs Euclid's
steps on these short numbers while Knuth's test proves that each quotient
is the true quotient of u by v, then applies all those steps to u and v
at once with a 2x2 matrix. If no step is proved, it does one ordinary
division step on u and v. When v fits in w bits, it ends with Euclid's
algorithm on single-precision numbers.
Cost: for n-bit inputs, at most O(n) outer iterations, each costing
O(n / w) word operations plus O(1) per Euclid step on the leading parts,
so O(n^2) bit operations in the worst case, like Euclid's algorithm; on
typical inputs each outer iteration replaces many multiprecision steps.
"""
u, v = abs(a), abs(b)
if u < v:
u, v = v, u
base = 1 << w
while v >= base: # @line L1
s = u.bit_length() - w # @line L2
uh, vh = u >> s, v >> s
A, B, C, D = 1, 0, 0, 1 # @line L3
while vh + C > 0 and vh + D > 0: # @line L4
q = (uh + A) // (vh + C) # @line L5
if q == 0 or q != (uh + B) // (vh + D): # @line L6
break
A, B, C, D = C, D, A - q * C, B - q * D # @line L7
uh, vh = vh, uh - q * vh # @line L8
if B == 0: # @line L9
u, v = v, u % v
else:
u, v = A * u + B * v, C * u + D * v # @line L10
while v != 0: # @line L11
u, v = v, u % v
return u # @line L12
test_lehmers_gcd.py 79 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/lehmers-gcd/. The authors are listed in prov/RN6MV6.jsonl.
import math
import random
from lehmers_gcd import lehmer_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 fib(n: int) -> int:
x, y = 0, 1
for _ in range(n):
x, y = y, x + y
return x
def test_example() -> None:
assert lehmer_gcd(27182818, 10000000, w=8) == 2
assert lehmer_gcd(1071, 462, w=4) == 21
assert lehmer_gcd(12257478, 400528, w=8) == math.gcd(12257478, 400528)
def test_zero_one_signs_equal() -> None:
for w in (2, 8, 64):
assert lehmer_gcd(0, 0, w) == 0
assert lehmer_gcd(0, 5, w) == 5
assert lehmer_gcd(-5, 0, w) == 5
assert lehmer_gcd(1, 10**50, w) == 1
assert lehmer_gcd(10**50, 1, w) == 1
assert lehmer_gcd(-12, -18, w) == 6
assert lehmer_gcd(-(10**40), 10**30, w) == 10**30
x = 3**200
assert lehmer_gcd(x, x, w) == x
assert lehmer_gcd(-x, x, w) == x
def test_powers_of_two() -> None:
for w in (3, 16, 64):
for i in range(0, 600, 37):
for j in range(0, 600, 41):
assert lehmer_gcd(2**i, 2**j, w) == 2 ** min(i, j)
assert lehmer_gcd(2**i * 3, 2**j * 9, w) == 2 ** min(i, j) * 3
def test_fibonacci() -> None:
for w in (2, 8, 32, 64):
for n in range(1, 1500, 47):
assert lehmer_gcd(fib(n + 1), fib(n), w) == 1
# gcd(F_m, F_n) = F_gcd(m, n)
assert lehmer_gcd(fib(2 * n), fib(3 * n), w) == fib(n)
def test_small_against_brute_force() -> None:
rng = random.Random(1938)
for _ in range(1500):
a = rng.randint(-3000, 3000)
b = rng.randint(-3000, 3000)
w = rng.choice([2, 3, 4, 6, 8])
assert lehmer_gcd(a, b, w) == brute_gcd(a, b)
def test_large_against_reference() -> None:
rng = random.Random(64)
for _ in range(500):
da = rng.choice([5, 30, 100, 300, 600])
db = rng.choice([5, 30, 100, 300, 600])
g = rng.randint(1, 10 ** rng.choice([1, 20, 150]))
a = g * rng.randint(-(10**da), 10**da)
b = g * rng.randint(-(10**db), 10**db)
w = rng.choice([2, 7, 16, 32, 64])
assert lehmer_gcd(a, b, w) == math.gcd(a, b)
What the next level needs
- the tests did not run in this build