18.16 Diophantine Equations
A Diophantine equation is an equation whose solutions must be integers.
18.16 Diophantine Equations
A Diophantine equation is an equation whose solutions must be integers. In algorithmic number theory, the most common form is the linear Diophantine equation:
ax + by = c
where a, b, and c are given integers, and we need integer values of x and y.
This problem is closely related to the Extended Euclidean Algorithm, modular inverses, and linear congruences. Once you know how to solve ax + by = gcd(a,b), solving ax + by = c becomes a matter of scaling and parameterizing the result.
Problem
Given integers:
a
b
c
find integers:
x
y
such that:
ax + by = c
Example:
15x + 25y = 10
One solution is:
x = 4
y = -2
because:
15×4 + 25×(-2) = 60 - 50 = 10
Not every equation has an integer solution. The GCD determines when solutions exist.
Solvability Criterion
The equation:
ax + by = c
has an integer solution if and only if:
gcd(a,b) | c
Let:
g = gcd(a,b)
If g does not divide c, no integer solution exists.
Example with no solution:
6x + 10y = 7
Here:
gcd(6,10) = 2
but:
2 ∤ 7
So the equation has no integer solution.
Example with solutions:
6x + 10y = 14
Here:
gcd(6,10) = 2
and:
2 | 14
So solutions exist.
Using Extended Euclid
The Extended Euclidean Algorithm gives integers x0 and y0 such that:
ax0 + by0 = gcd(a,b)
If:
g = gcd(a,b)
and:
c = gq
then multiply both sides by q:
a(x0q) + b(y0q) = c
So:
x = x0q
y = y0q
is one solution.
Example
Solve:
15x + 25y = 10
First compute:
gcd(15,25) = 5
Extended Euclid gives:
15(2) + 25(-1) = 5
Since:
10 / 5 = 2
multiply the coefficients by 2:
x = 4
y = -2
Check:
15(4) + 25(-2) = 60 - 50 = 10
This gives one solution.
General Solution
Once one solution is known, all solutions can be generated.
Suppose:
ax0 + by0 = c
Let:
g = gcd(a,b)
Then all integer solutions are:
x = x0 + (b/g)t
y = y0 - (a/g)t
for any integer t.
Why this works:
a[x0 + (b/g)t] + b[y0 - (a/g)t]
=
ax0 + by0 + abt/g - abt/g
=
c
The added terms cancel.
Implementation
from math import gcd
def extended_gcd(a, b):
if b == 0:
return abs(a), 1 if a >= 0 else -1, 0
g, x1, y1 = extended_gcd(b, a % b)
x = y1
y = x1 - (a // b) * y1
return g, x, y
def solve_diophantine(a, b, c):
g, x0, y0 = extended_gcd(a, b)
if c % g != 0:
return None
scale = c // g
x = x0 * scale
y = y0 * scale
return x, y, g
Example:
print(solve_diophantine(15, 25, 10))
Output:
(4, -2, 5)
The return value gives one solution and the GCD.
Generating All Solutions
def diophantine_solution(a, b, c, t):
solved = solve_diophantine(a, b, c)
if solved is None:
return None
x0, y0, g = solved
x = x0 + (b // g) * t
y = y0 - (a // g) * t
return x, y
Example:
for t in range(-2, 3):
print(t, diophantine_solution(15, 25, 10, t))
Possible output:
-2 (-6, 4)
-1 (-1, 1)
0 (4, -2)
1 (9, -5)
2 (14, -8)
All satisfy:
15x + 25y = 10
Handling Negative Coefficients
The formula still works when a, b, or c is negative, but implementation must handle signs carefully.
Example:
-6x + 10y = 14
This has solutions because:
gcd(-6,10) = 2
and:
2 | 14
A robust extended_gcd should return a positive gcd while allowing coefficients to satisfy the original equation.
For production use, verify:
g, x, y = extended_gcd(a, b)
assert a * x + b * y == g
This assertion catches sign bugs.
Finding Nonnegative Solutions
Many applications require:
x ≥ 0
y ≥ 0
Given the general solution:
x = x0 + (b/g)t
y = y0 - (a/g)t
we need values of t satisfying both inequalities.
Assume a > 0 and b > 0.
Then:
x0 + (b/g)t ≥ 0
and:
y0 - (a/g)t ≥ 0
So:
t ≥ -x0 / (b/g)
and:
t ≤ y0 / (a/g)
Using integer ceiling and floor division:
def ceil_div(a, b):
return -((-a) // b)
def floor_div(a, b):
return a // b
def nonnegative_solutions(a, b, c):
solved = solve_diophantine(a, b, c)
if solved is None:
return []
x0, y0, g = solved
step_x = b // g
step_y = a // g
t_min = ceil_div(-x0, step_x)
t_max = floor_div(y0, step_y)
result = []
for t in range(t_min, t_max + 1):
x = x0 + step_x * t
y = y0 - step_y * t
result.append((x, y))
return result
Example:
print(nonnegative_solutions(6, 10, 42))
This returns all nonnegative integer solutions to:
6x + 10y = 42
Counting Nonnegative Solutions
If there may be many solutions, return only the count.
def count_nonnegative_solutions(a, b, c):
solved = solve_diophantine(a, b, c)
if solved is None:
return 0
x0, y0, g = solved
step_x = b // g
step_y = a // g
t_min = ceil_div(-x0, step_x)
t_max = floor_div(y0, step_y)
return max(0, t_max - t_min + 1)
This runs in constant time after Extended Euclid.
Example: Coin Problem
Suppose coins have values:
6
10
How many ways can we form:
42
using nonnegative counts?
We solve:
6x + 10y = 42
Divide by 2:
3x + 5y = 21
Try values of y:
y = 0 → x = 7
y = 3 → x = 2
So the solutions are:
(x,y) = (7,0), (2,3)
There are two ways.
Relation to Linear Congruences
The equation:
ax + by = c
can be viewed modulo b:
ax ≡ c (mod b)
So solving a two-variable Diophantine equation is equivalent to solving a linear congruence for one variable, then recovering the other.
This relationship is useful when one variable has range constraints.
Application: Lattice Points on a Line
A line:
ax + by = c
has integer lattice points exactly when:
gcd(a,b) | c
The general solution describes all lattice points on that line.
This appears in computational geometry, grid algorithms, and integer optimization.
Application: Scheduling
Suppose two periodic processes produce offsets in multiples of a and b. To determine whether a target offset c is reachable, solve:
ax + by = c
Reachability again depends on:
gcd(a,b) | c
Common Mistakes
Forgetting the Divisibility Check
If gcd(a,b) does not divide c, no solution exists.
Returning Only One Solution
Extended Euclid gives one solution. Most Diophantine equations have infinitely many integer solutions.
Getting the Sign Wrong
The general solution must use opposite signs:
x = x0 + (b/g)t
y = y0 - (a/g)t
Mishandling Nonnegative Bounds
Do not brute force x or y unless the range is tiny. Convert the inequalities into bounds on t.
Assuming Positive Coefficients
Negative coefficients change inequality directions. Normalize or handle signs explicitly.
Testing
Basic tests:
def test_diophantine():
solved = solve_diophantine(15, 25, 10)
assert solved is not None
x, y, _ = solved
assert 15 * x + 25 * y == 10
assert solve_diophantine(6, 10, 7) is None
Check generated solutions:
def test_general_solution():
a, b, c = 15, 25, 10
for t in range(-10, 11):
x, y = diophantine_solution(a, b, c, t)
assert a * x + b * y == c
Check nonnegative counting with brute force for small values:
def brute_count(a, b, c):
count = 0
for x in range(c + 1):
for y in range(c + 1):
if a * x + b * y == c:
count += 1
return count
def test_nonnegative_count():
for a in range(1, 20):
for b in range(1, 20):
for c in range(0, 100):
assert count_nonnegative_solutions(a, b, c) == brute_count(a, b, c)
This catches errors in the interval arithmetic.
Takeaway
A linear Diophantine equation ax + by = c has integer solutions exactly when gcd(a,b) divides c. The Extended Euclidean Algorithm gives one solution to ax + by = gcd(a,b), and scaling gives one solution to the target equation. All other solutions are obtained by shifting along the line using x = x0 + (b/g)t and y = y0 - (a/g)t. With this parameterization, nonnegative solution counting and range-constrained variants reduce to simple integer interval calculations.