Binary GCD in Python
Computes gcd using only subtraction, parity tests and shifts, avoiding division.
- File
binary_gcd.py, 52 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/binary-gcd/. The authors are listed in prov/NZK9TA.jsonl.
"""Binary GCD (Stein's algorithm).
Computes gcd(a, b) with parity tests, shifts and subtractions only, with no
division. Inputs may be any integers; the result is non-negative, and
gcd(0, 0) = 0. On inputs of at most n binary digits the loop runs at most
2n + 1 times, and each pass costs O(n) bit operations, so the total is
O(n^2) bit operations. Extra space is O(n) bits for the two working values.
"""
def trailing_zeros(x: int) -> int:
"""The exponent of 2 in x, for x > 0: the number of trailing zero digits."""
return (x & -x).bit_length() - 1
def binary_gcd(a: int, b: int) -> int:
"""Return gcd(a, b) >= 0 using shifts and subtractions."""
a, b = abs(a), abs(b) # @line L1
if a == 0: # @line L2
return b
if b == 0: # @line L3
return a
k = trailing_zeros(a | b) # @line L4
a >>= trailing_zeros(a) # @line L5
while True: # @line L6
b >>= trailing_zeros(b) # @line L7
if a > b: # @line L8
a, b = b, a
b -= a # @line L9
if b == 0: # @line L10
break
return a << k # @line L11
def loop_passes(a: int, b: int) -> int:
"""The number of passes of the loop of lines L6 to L10 on (a, b)."""
a, b = abs(a), abs(b)
if a == 0 or b == 0:
return 0
a >>= trailing_zeros(a)
passes = 0
while True:
passes += 1
b >>= trailing_zeros(b)
if a > b:
a, b = b, a
b -= a
if b == 0:
return passes
test_binary_gcd.py 102 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/binary-gcd/. The authors are listed in prov/NZK9TA.jsonl.
import math
import random
from binary_gcd import binary_gcd, loop_passes, trailing_zeros
def fib(n: int) -> int:
a, b = 0, 1
for _ in range(n):
a, b = b, a + b
return a
def brute_gcd(a: int, b: int) -> int:
"""The largest d that divides both, by direct search (small inputs only)."""
a, b = abs(a), abs(b)
if a == 0 and b == 0:
return 0
for d in range(max(a, b), 0, -1):
if a % d == 0 and b % d == 0:
return d
raise AssertionError("unreachable")
def test_example() -> None:
assert binary_gcd(240, 46) == 2
assert loop_passes(240, 46) == 5
assert binary_gcd(1071, 462) == 21
def test_trailing_zeros() -> None:
for x in range(1, 2000):
t = 0
y = x
while y % 2 == 0:
y //= 2
t += 1
assert trailing_zeros(x) == t
def test_zero_one_negative_equal() -> None:
assert binary_gcd(0, 0) == 0
assert binary_gcd(0, 7) == 7
assert binary_gcd(-7, 0) == 7
assert binary_gcd(0, -8) == 8
assert binary_gcd(1, 10**30) == 1
assert binary_gcd(-12, 18) == 6
assert binary_gcd(12, -18) == 6
assert binary_gcd(-12, -18) == 6
assert binary_gcd(42, 42) == 42
assert binary_gcd(-42, 42) == 42
assert binary_gcd(-(2**63), 2**62) == 2**62
def test_powers_of_two() -> None:
for i in range(0, 200, 7):
for j in range(0, 200, 11):
assert binary_gcd(2**i, 2**j) == 2 ** min(i, j)
assert binary_gcd(3 * 2**i, 5 * 2**j) == 2 ** min(i, j)
assert binary_gcd(6 * 2**i, 9 * 2**j) == 3 * 2 ** min(i + 1, j)
def test_fibonacci() -> None:
for k in range(1, 400):
a, b = fib(k + 2), fib(k + 1)
assert binary_gcd(a, b) == 1
assert binary_gcd(fib(2 * k), fib(k)) == fib(k)
def test_pass_bound() -> None:
rng = random.Random(7)
for _ in range(2000):
a = rng.randrange(1, 2**64)
b = rng.randrange(1, 2**64)
assert loop_passes(a, b) <= 1 + (a * b).bit_length()
for n in range(1, 200):
# (2^n - 1, 1) needs one pass for each binary digit
assert loop_passes(2**n - 1, 1) == n
def test_random_small_against_direct_search() -> None:
rng = random.Random(1)
for _ in range(1500):
a = rng.randint(-300, 300)
b = rng.randint(-300, 300)
assert binary_gcd(a, b) == brute_gcd(a, b)
def test_random_large_against_math_gcd() -> None:
rng = random.Random(2)
for _ in range(500):
g = rng.getrandbits(rng.randint(1, 400)) + 1
a = g * rng.getrandbits(rng.randint(0, 1200)) * rng.choice([-1, 1])
b = g * rng.getrandbits(rng.randint(0, 1200)) * rng.choice([-1, 1])
assert binary_gcd(a, b) == math.gcd(a, b)
big = 10**300 + 7
for m in range(1, 50):
assert binary_gcd(big * m, big * (m + 1)) == big * math.gcd(m, m + 1)
assert binary_gcd(big * 2**m, big * 3**m) == big
What the next level needs
- the tests did not run in this build