Finding a solution of a linear Diophantine equation in Python
Uses the extended Euclidean algorithm to find one integer solution of a x plus b y equals c, or reports that none exists.
- File
linear_diophantine_solution.py, 32 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/linear-diophantine-solution/. The authors are listed in prov/747E8M.jsonl.
def solve_linear_diophantine(a: int, b: int, c: int) -> tuple[int, int] | None:
"""Return one integer pair (x, y) with a*x + b*y == c, or None if there is none.
Precondition: a, b and c are any Python integers; negative values and 0
are allowed. If a == b == 0, the answer is (0, 0) when c == 0 and None
otherwise. Otherwise a solution exists exactly when gcd(a, b) divides c,
and the pair returned is (s * c/g, t * c/g) for the Bezout coefficients
s and t that the extended Euclidean algorithm finds.
Cost: O(1 + log min(|a|, |b|)) loop iterations. Python integers have
arbitrary precision; for inputs of at most n bits the loop takes O(n^2)
bit operations and the final products O(n^2) more.
"""
if a == 0 and b == 0: # @line L1
return (0, 0) if c == 0 else None # @line L2
r0, s0, t0, r1, s1, t1 = abs(a), 1, 0, abs(b), 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
g = r0 # @line L9: g == |a|*s0 + |b|*t0 == gcd(a, b) > 0
if c % g != 0: # @line L10
return None
m = c // g
x = s0 * m if a >= 0 else -s0 * m # @line L11
y = t0 * m if b >= 0 else -t0 * m
return (x, y) # @line L12
test_linear_diophantine_solution.py 133 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/linear-diophantine-solution/. The authors are listed in prov/747E8M.jsonl.
import math
import random
from linear_diophantine_solution import solve_linear_diophantine
def brute_solvable(a: int, b: int, c: int) -> bool:
"""Decide solvability by direct search (small inputs only)."""
if a == 0 and b == 0:
return c == 0
if b == 0:
return c % a == 0
# x and x + |b| give the same residue of c - a*x modulo b.
return any((c - a * x) % b == 0 for x in range(abs(b)))
def check(a: int, b: int, c: int) -> tuple[int, int] | None:
r = solve_linear_diophantine(a, b, c)
if r is not None:
x, y = r
assert a * x + b * y == c
return r
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 solve_linear_diophantine(1071, 462, 63) == (-9, 21)
assert solve_linear_diophantine(6, 9, 7) is None
assert solve_linear_diophantine(3, 5, 1) == (2, -1)
def test_zero_and_signs() -> None:
assert solve_linear_diophantine(0, 0, 0) == (0, 0)
assert solve_linear_diophantine(0, 0, 5) is None
assert check(5, 0, 15) == (3, 0)
assert check(0, 5, 15) == (0, 3)
assert check(5, 0, 7) is None
assert check(-5, 0, 15) == (-3, 0)
assert check(-4, 6, 10) is not None
assert check(-4, -6, -10) is not None
assert check(4, 6, 0) == (0, 0)
assert check(1, 1, 1) is not None
assert check(7, 7, 14) is not None
assert check(7, 7, 15) is None
assert check(-7, 7, 0) is not None
def test_powers_of_two() -> None:
for i in range(0, 70, 7):
for j in range(0, 70, 5):
a, b = 2**i, 2**j
g = 2 ** min(i, j)
assert check(a, b, g) is not None
assert check(a, b, 3 * g) is not None
assert (check(a, b, g + 1) is None) == (g > 1)
def test_fibonacci_worst_case() -> None:
assert solve_linear_diophantine(832040, 1346269, 1) == (514229, -317811)
for n in range(2, 300, 7):
a, b = fib(n), fib(n + 1)
assert check(a, b, 1) is not None
assert check(b, a, -12345) is not None
def test_random_small_against_brute_force() -> None:
rng = random.Random(747)
for _ in range(3000):
a = rng.randint(-60, 60)
b = rng.randint(-60, 60)
c = rng.randint(-200, 200)
r = check(a, b, c)
assert (r is not None) == brute_solvable(a, b, c)
def test_random_large() -> None:
rng = random.Random(2024)
for _ in range(500):
digits = rng.choice([1, 5, 20, 100, 300, 500])
a = rng.randint(-(10**digits), 10**digits)
b = rng.randint(-(10**digits), 10**digits)
g = math.gcd(a, b)
if rng.random() < 0.5 or g == 0:
c = rng.randint(-(10**digits), 10**digits)
else:
c = g * rng.randint(-(10**digits), 10**digits)
r = check(a, b, c)
expect = (c == 0) if g == 0 else (c % g == 0)
assert (r is not None) == expect
def test_common_factor_large() -> None:
rng = random.Random(9)
for _ in range(100):
f = rng.randint(2, 10**200)
a = f * rng.randint(1, 10**150)
b = -f * rng.randint(1, 10**150)
g = math.gcd(a, b)
assert check(a, b, 7 * g) is not None
assert check(a, b, g + 1) is None
def test_general_solution_window() -> None:
rng = random.Random(31)
for _ in range(300):
a = rng.randint(-30, 30)
b = rng.randint(-30, 30)
c = rng.randint(-60, 60)
if a == 0 or b == 0:
continue
r = solve_linear_diophantine(a, b, c)
sols = {(x, (c - a * x) // b) for x in range(-100, 101) if (c - a * x) % b == 0}
if r is None:
assert not sols
continue
x0, y0 = r
g = math.gcd(a, b)
dx, dy = b // g, -a // g
expect = set()
for k in range(-300, 301):
x, y = x0 + k * dx, y0 + k * dy
if -100 <= x <= 100:
expect.add((x, y))
assert sols == expect
What the next level needs
- the tests did not run in this build
Read the full page about Finding a solution of a linear Diophantine equation