Extended Euclidean algorithm in Python
Computes gcd(a, b) together with integers x, y satisfying ax + by = gcd(a, b).
- File
extended_euclidean_algorithm.py, 37 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/extended-euclidean-algorithm/. The authors are listed in prov/QQCABS.jsonl.
"""Extended Euclidean algorithm.
Returns (g, x, y) with g = gcd(a, b) >= 0 and a*x + b*y = g, for any
integers a and b. gcd(0, 0) is 0, returned as (0, 1, 0). For a, b >= 1 the
coefficients satisfy |x| <= b/g and |y| <= a/g, and |x| <= b/(2g),
|y| <= a/(2g) unless b divides a. The loop runs O(log min(|a|, |b|)) times;
on inputs of n binary digits it uses O(n^2) bit operations. Extra space is
a constant number of integers no larger than the inputs.
"""
def extended_euclidean_algorithm(a: int, b: int) -> tuple[int, int, int]:
"""Return (g, x, y) with g = gcd(a, b) >= 0 and a*x + b*y = g."""
r0, r1 = abs(a), abs(b) # @line L1
s0, s1 = 1, 0 # @line L2
t0, t1 = 0, 1 # @line L3
while r1 != 0: # @line L4
q = r0 // r1 # @line L5
r0, r1 = r1, r0 - q * r1 # @line L6
s0, s1 = s1, s0 - q * s1 # @line L7
t0, t1 = t1, t0 - q * t1 # @line L8
x = -s0 if a < 0 else s0 # @line L9
y = -t0 if b < 0 else t0
return r0, x, y # @line L10
def modular_inverse(a: int, m: int) -> int:
"""Return the inverse of a modulo m >= 1, in 0..m-1; raise if gcd(a, m) != 1."""
if m < 1:
raise ValueError("modulus must be positive")
g, x, _ = extended_euclidean_algorithm(a, m)
if g != 1:
raise ValueError("a is not invertible modulo m")
return x % m
test_extended_euclidean_algorithm.py 111 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/extended-euclidean-algorithm/. The authors are listed in prov/QQCABS.jsonl.
import math
import random
import pytest
from extended_euclidean_algorithm import extended_euclidean_algorithm as egcd
from extended_euclidean_algorithm import modular_inverse
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 check(a: int, b: int, g_ref: int) -> None:
g, x, y = egcd(a, b)
assert g == g_ref
assert a * x + b * y == g
if a != 0 and b != 0:
assert abs(x) <= abs(b) // g
assert abs(y) <= abs(a) // g
if abs(a) % abs(b) != 0:
assert 2 * g * abs(x) <= abs(b)
assert 2 * g * abs(y) <= abs(a)
else:
assert x == 0 and abs(y) == 1
def test_example() -> None:
assert egcd(240, 46) == (2, -9, 47)
assert egcd(1071, 462) == (21, -3, 7)
def test_zero_one_negative_equal() -> None:
assert egcd(0, 0) == (0, 1, 0)
assert egcd(0, 5) == (5, 0, 1)
assert egcd(5, 0) == (5, 1, 0)
assert egcd(0, -5) == (5, 0, -1)
assert egcd(-5, 0) == (5, -1, 0)
assert egcd(7, 7) == (7, 0, 1)
assert egcd(1, 1) == (1, 0, 1)
assert egcd(-240, 46) == (2, 9, 47)
assert egcd(240, -46) == (2, -9, -47)
assert egcd(-240, -46) == (2, 9, -47)
for a in range(-30, 31):
for b in range(-30, 31):
check(a, b, brute_gcd(a, b))
def test_powers_of_two() -> None:
for i in range(0, 150, 7):
for j in range(0, 150, 11):
check(2**i, 2**j, 2 ** min(i, j))
check(3 * 2**i, -5 * 2**j, 2 ** min(i, j))
def test_fibonacci() -> None:
for k in range(3, 300):
a, b = fib(k + 1), fib(k)
g, x, y = egcd(a, b)
assert g == 1 and a * x + b * y == 1
# by Cassini's identity the coefficients are consecutive Fibonacci numbers
assert {abs(x), abs(y)} == {fib(k - 2), fib(k - 1)}
def test_random_small_against_direct_search() -> None:
rng = random.Random(1)
for _ in range(1500):
a = rng.randint(-500, 500)
b = rng.randint(-500, 500)
check(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])
check(a, b, math.gcd(a, b))
big = 10**300 + 7
for m in range(1, 40):
check(big * m, big * (m + 1), big)
def test_modular_inverse_against_pow() -> None:
rng = random.Random(3)
for _ in range(800):
m = rng.choice([rng.randint(1, 1000), rng.getrandbits(900) + 1])
a = rng.randint(-(10**300), 10**300)
if math.gcd(a, m) == 1:
assert modular_inverse(a, m) == pow(a, -1, m)
else:
with pytest.raises(ValueError):
modular_inverse(a, m)
What the next level needs
- the tests did not run in this build