Linear Diophantine equations

A linear Diophantine equation is an equation of the form $ax + by = c$ where $a, b, c \in \mathbb{N}$. They require a specific set of solutions where all the coefficients are integers. When can a linear Diophantine equation be solved? The equation only has solutions if $\textrm{GCD}(a, b)$ is a divisor of $c$, otherwise there are no solutions. How many solutions are there? Well, if they can be solved, there are always infinitely many solutions.

Diophantine equations are notoriously difficult to solve, however this linear version is rather simple (in relative terms). A single solution can be found by using the extended Euclid's algorithm. If we have found one solution $(x, y)$, then all solutions can be found with $(x+kv, y+ku)$, where $k$ is any integer number, and $u$ and $v$ are the quotients of $a$ and $b$ (respectively) by the greatest common divisor of $a$ and $b$. The steps to solve a linear Diophantine equation are as follows:

  1. First check if the equation has any solutions by checking if $\textrm{GCD}(a, b)$ is a divisor of $c$.
  2. If it has solutions, we divide the equation by $\textrm{GCD}(a, b)$ to get the equivalent form $Ax + By = C$.
  3. Find one solution of $Ax + By = C$ by using the extended Euclid's algorithm.
  4. Write all the solution by using the single solution which is a variation of all solutions.

Example: Find all solutions of $3900x + 744y = 36$.

First we find $\textrm{GCD}(3900, 744) = 12\ |\ 36$, which means that there are solutions. Now we divide $a, b, c$ by $12$ to find the general form $325x + 62y = 3$. We will solve the general form with the extended Euclid's algorithm.

$325$ $62$
$325$ $1$ $0$
$62$ $0$ $1$
$15$ $1$ $-5$ $R_1 - 5R_2$
$2$ $-4$ $21$ $R_2-4R_3$
$1$ $29$ $-152$ $R_3-7R_4$
$(\textrm{blow-up})\quad 3$ $87$ $-456$ $3R_5$

From the table we can read the first solution, which is $(87, -456)$. However, notice that $325\cdot87 + 62(-456) = 3$. To get back to our initial equation, we have to multiply it by $12$. This gives the equation $3900\cdot87 + 744(-456) = 36$. Using this single solution, we can write all the solutions as $3900(87 - 62k) + 744(-456 + 325k) = 36$. In other words: $$ \begin{align} \begin{cases} x_k &= 87 & - 62k \\ y_k &= -456 &+ 325k \end{cases} \end{align} $$

Note that if the equation has different signs, the $\textrm{GCD}$ is calculated with the positive values of $a, b$. Then the it is solved in the same way. However, the final solution needs to have the signs adjusted. Just something to keep in mind!

The following Python code will solve linear Diophantine equations.

def gcd(a, b):
    if b > a:
        a, b = b, a
    while b != 0:
        a, b = b, a - a//b*b
    return a


def extended_euclid(a, b):
    sign_a = 1 if a >= 0 else -1
    sign_b = 1 if b >= 0 else -1
    a = abs(a)
    b = abs(b)
    x1, y1 = 1, 0
    x2, y2 = 0, 1
    swap = b > a
    if swap:
        a, b = b, a
    while b != 0:
        q = a // b
        a, b = b, a - q*b
        x1, y1, x2, y2 = x2, y2, x1 - q*x2, y1 - q*y2
    if swap:
        x1, y1 = y1, x1
    return sign_a * x1, sign_b * y1


sign = lambda x: 1 if x >= 0 else -1 
def solve_linear_diophantine(a, b, c):
    d = gcd(a, b)
    if c % d != 0:
        print("Warning: GCD(a, b) is not a divisor of c. The equation has no solutions.")
        return -1
    A, B, C = a // d, b // d, c // d
    x, y = extended_euclid(A, B)
    x, y = x * c/d, y * c/d
    print('a={}, b={}, c={}, A={}, B={}, C={}, x={}, y={}'.format(a,b,c,A,B,C,x,y))
    return lambda t: (x + abs(b/d*t) * sign(x), y + abs(a/d*t) * sign(y))

The solve_linear_diophantine(a, b, c) function, returns a function with a generator for all the solutions. This is a function f(t). If we test this function, with the coefficients from the example, we will find that:

>>> f = lambda x, y: 3900*x + 744*y
>>> g = solve_linear_diophantine(3900, 744, 36)
a=3900, b=744, c=36, A=325, B=62, C=3, x=87.0, y=-456.0
>>> f(*g(1))
36.0
>>> f(*g(2))
36.0
>>> f(*g(3))
36.0