# Diophantine equation ax + by = c

Of course with a, b, c, x, y in integers ∈ Z.

First of all, if p divides a and b, it should also divide c, else equation has no solution. We can then divide everything by p.
We care from now only of equations where a and b coprimes, that is without any common divisor but 1.
Let ax0 + by0 = c one solution. Then if ax + by = c is another one, by substracting comes :
a(x - x0) + b(y - y0) = 0. As a and b coprimes, a divides y - y0, that is y - y0 = ka or :
y = y0 + ka then x = x0 - kb

Consider now equation ax + by = 1, still with a and b coprime. The Bézout theorem says that two numbers a and b are coprime if and only if there are two numbers u and v (in Z) with au + bv = 1. This proves that equation ax + by = c, GCD(a, b) = 1 has always solutions. If (u, v) is solution of ax + by = 1, then (cu, cv) is solution of ax + by = c.
Summary :

 Equation ax + by = c has solutions if and only if p = GCD(a,b) divides c Then let a = pa', b = pb' c = pc', if x0,y0 is some solution of  a'x + b'y = 1, then all solutions of ax + by = c are :      x = c'x0 - kb'     y = c'y0 + ka' for all k ∈ Z

The problem is now to find some specific solution.
In very simple cases, we can just guess or try, but equations like 97x + 29y = 4000 already requires a more efficient method.

## Variable substitution

Let's take as example equation 97x + 29y = 4000.

As 97 = 29×3 + 10 and 4000 = 29×137+27, equation 97x + 29y = 4000 can be written :
10x + 29(y + 3x - 137) = 27 that is substituting z = y + 3x - 137 : 10x + 29z = 27, which is a simpler equation.

We can repeat in writing : 10(x + 2z - 2) + 9z = 7 that is 10t + 9z = 7 then :
t + 9(z + t) = 7, or t + 9u = 7 which has obvious solution t = 7, u = 0

Doing variable substitution backwards, we get : z = u - t = -7, then x = t + 2 - 2z = 23 and finally
y = z + 137 - 3x = 61, hence a specific solution x0 = 23, y0=61,
and all solutions x = 23 - 29k, y = 61 + 97k.

## Congruences

An interresting method is to write ax + by = 1 as a congruence : ax ≡ 1 mod (b).
As a and b coprime, the Euler-Fermat theorem says aφ(b) ≡ 1 mod (b)
with φ(b) being the totient or Euler function of b.
Multiplying the two terms of ax ≡ 1 mod (b) by aφ(b)-1 results into...
The general solution of ax + by = c :

 x ≡caφ(b)-1 mod (b), y = (c - ax)/b

Remains to calculate φ(b) (and for those who don't believe, to prove the Euler-Fermat theorem)
A simple method to calculate φ(b) from the prime decomposition of b :
Let b = ∏piai, then φ(b) = b∏(1 - 1/pi) = ∏(piai - piai-1). Specifically, when p prime, φ(p) = p - 1.

And for the above equation : x ≡ 4000×97φ(29)-1 mod (29), that is 4000×9727
Hopefully, no need of calculate huge numbers as all calculations are modulo 29, The numbers will never be more than 29² = 841.
4000 ≡ 27 mod (29) and 97 ≡ 10 mod (29) already simplfy calculations into :
4000×9727 ≡ 27×1027 mod (29), Pfffiuu !
Remains to effectively calculate 1027 mod (29). We use a method known as "Indian exponentiation"

Algorithm for calculation of V = AB

 Initialize V = 1 while B ≥ 1     - if B odd, multiply V by A  and decrease B by one     - else, square A and divide B by 2

Or a variante in binary :

 Initialize V = 1 For each bit in exposant B, beginning from the most significant bit :      - square V     - if the bit is 1, multiply V by A

which gives here 27 = 11011(2) and 1027 = (((1²×10)²×10)²)²×10)²×10
As every multiplication or squaring is modulo 29 :
102 = 100 ≡ 13 mod (29)
103 ≡ 13×10 = 130 ≡ 14 mod (29)
106 ≡ 142 =196 ≡ 22 mod (29)
1012 ≡ 222 =484 ≡ 20 mod (29)
1013 ≡ 20×10 = 200 ≡ 26 mod (29)
1026 ≡ 262 = 676 ≡ 9 mod (29)
1027 ≡ 9×10 = 90 ≡ 3 mod (29)
And finally x ≡ 27×3 = 81 ≡ 23 mod (29)

## Euclidean Algorithm

The previous calculation, although using an elegant relation, are not very efficient, at least because they use the prime factorization of b (to get φ(b)).
The Euclidean algorithm, used to calculate the GCD of a and b, gives also as a result the representation of a/b as continued fraction. The study of such fractions gives a variant of Euclidean algorithm resulting into the solution of ax + by = 1.
First of all, note that by variable substitution (exchange of x and y and sign changes) we could assert a > b > 0. In such conditions :

 Let : r0 = a, r1 = b (a > b > 0) P0 = 0, P1 = 1 and Q0 = 1, Q1 = 0 Then repeat from i = 2     - Integer division of ri-2 by ri-1 : ri-2 = ri-1qi + ri      - Calculate Pi = -qiPi-1 + Pi-2 and Qi = -qiQi-1 + Qi-2      - until rn = 0 p = rn-1 is GCD of a and b (Euclidean algorithm)  Let a' = a/p and b' = b/p, x0 = Qn-1, y0 = Pn-1 is solution of a'x + b'y = 1
Proof of this algorithm.

With 97x + 29y = 4000 this results into :
97 = 29×3 + 10, r2 = 10, q2 = 3, P2 = -3, Q2 = 1
29 = 10×2 + 9, r3 = 9, q3 = 2, P3 = 7, Q3 = -2
10 = 9×1 + 1, r4 = 1, q4 = 1, P4 = -10, Q4 = 3
9 = 1×9 + 0, r5 = 0, that's all and x = 3, y = -10 is solution of 97x + 29y = 1
Multiplying by c = 4000, which is a multiple of p = r4 = 1, we get x = 12000 + 29k, y = -40000 - 97k
If we let m = k + [12000/29] = k + 413, we get back x = 23 + 29m, y = 61 - 97m as above.

Script to solve ax + by = c with this algorithm.

# Equations with more than 2 unknowns

Take as example 15x + 21y + 35z = 200.
To solve such an equation, just consider all unknowns but two as parameters, and put them on the other side.
That is : 15x + 21y = 200 - 35z = c and we have to solve 15x + 21y = c as a function of parameter c.
Of course a little problem as 15 and 21 are not coprime. Hence c should be a multiple of GCD(15,21) = 3.
That is 200 - 35z ≡ 0 mod (3), or z ≡ 1 mod(3). Which could be written as z = 3m + 1.
The real unknown is m instead of z.
The equation then becomes 15x + 21y = 165 - 3×35m and dividing everything by 3 : 5x + 7y = 55 - 35m.
As x0 = 3, y0 = -2 is solution of 5x + 7y = 1 (exercice...), the general solution of this equation is
x = 3(55 - 35m) + 7k, y = -2(55 - 35m) - 5k and finally

 x = 165 - 105m - 7k y = -110 + 70m + 5k z = 3m + 1 k and m independant in Z

Let's give general formulas for solving equation ax + by + cz = D

 Let p = GCD(a,b), a' = a/p, b' = b/p Let u0 and v0 some solution of a'u + b'v = c  z0, t0 some solution of cz + pt = D x0, y0 some solution of a'x + b'y = t0 The general solution of ax + by + cz = D is : x = x0 + b'k - u0m y = y0 - a'k - v0m z = z0 + pm with k and m ∈ Z independant
A Script with this method.

Direct application to 15x + 21y + 35z = 200 really gives :
p = 3, a' = 5, b' = 7 then 5u + 7v = 35 and u0 = 0 v0 = 5 (for instance)
35z + 3t = 200 z0 = 1, t0 = 55
5x + 7y = 55 : x0 = 11, y0 = 0 and finally :

 x = 11 + 7k y = 0 - 5k - 5m z = 1 + 3m
which seems to be very different of the previously found solution.
But just substitute k = 22 - k' - 15m gives exactly the same solution.

This method is quite tedious as it requires to successively solve several equations in two unknowns, that is to perform the Euclide's algorithm several times.
It therefore becomes impracticable with more than 3 unknowns.
We then could efficiently solve this equation in considering it as a "system" (hum) of one équation (!) in n unknowns...

#### Generalization of the Euclide's algorithm

Firsts of all, note that the GCD of n numbers could be obtained in "crunching" the successive Euclide's algorithms.
We could consider the Euclide's algorithm (for GCD only) as successive linear combinations of the numbers :
For instance to get the GCD(15, 21, 35)
The smallest of the 3 numbers is 15, we then do an Euclidean division of all others by 15, keeping just the remainders.
That is GCD(15, 21, 35) = GCD(15, 21 - 15 , 35 - 2×15) = GCD(15, 6, 5)
We go on with the new smallest number 5 : GCD(15, 6, 5) = GCD(15 - 3×5, 6 - 5, 5) = GCD(0, 1, 5)
And one more time (with the smallest non null) : GCD(0, 1, 5) = GCD(0, 1, 5 - 5×1) = GCD(0, 1, 0) = 1
This process allways ends as for each step, all the remainders are < smallest non null value.
This one is then strictly decreasing.
The extended Euclide's algorithm (which also gives the Bézout relation ax + by = 1) just "operates" simultaneously the same reduction to the lines of a unit matrix, that is we just consider the matrix obtained by concatenating the column a, b and a unit matrix :
```(a | 1 0)
(b | 0 1)```
Generalizing this reduction process to our equation 15x + 21y + 35z = 200, we perform a "row Euclidean reduction" of the matrix :
```(15 | 1  0  0)
(21 | 0  1  0)
(35 | 0  0  1) resulting into :

(15 | 1  0  0)
(6  |-1  1  0) = minus [21/15]=1 times the 1st line
(5  |-2  0  1) = minus [35/15]=2 times the 1st line

(0  | 7  0 -3) = minus [15/5]=3 times the 3rd line
(1  | 1  1 -1) = minus [6/5]=1 times the 3rd line
(5  |-2  0  1)

and finally :

(0  | 7  0 -3)
(1  | 1  1 -1)
(0  |-7 -5  6) = minus [5/5]=5 times the 2nd line```
The only line with a non null value in the 1st column (= GCD = 1) gives directly the Bézout relation : 15×1 + 21×1 + 35×(-1) = PGCD = 1
That is a solution of our equation : x = 200×1, y = 200×1, z = 200×(-1)

The lines beginning with 0 in the above matrix give solutions of the homogenous equation 15x + 21y + 35z = 0
That is (x = 7, y = 0, z = -3) and (x = -7, y = -5, z = 6), which are linearily independant, and we may as well replace them by two independant linear combinations, for instance u = (7, -5, 0) and v = (0, -5, 3)
We can then freely add any linear combination of u and v : ku + mv.

This gives the general solution of our equation :

 x = 200 + 7k + 0m   y = 200 - 5k - 5m   z = -200 + 0k + 3m
Which is, at a shift (m = 67 + m', k = -27 + k'), the values found previously.

This method may be used for any number of unknowns.

# Systems of equations

And even more : we could generalize to a system of equations, without first reducing it to a single equation.

Thus consider to write the system as matrix form (A)(X) = (B).
Perform the previous row Euclidean reduction on matrix (At|I), obtained by concatenation of the transposed of (A) and a unit matrix I.
We then get a matrix (R|T), in which T is the transforming matrix : (R|T) = (T)(At|I)
(T) is invertible, as we just made linear combinations of rows.
This also results into (T)(At) = (R), hence :
(A)(Tt) = (Rt), and finally, as T is invertible : (A) = (Rt)(Tt)-1
The determinant of T being ±1 (as product of matrices with determinant ±1, unimodular matrices), the inverse (Tt)-1 is a matrix in integers.
Equation (A)(X) = (B) is then (Rt)(Tt)-1(X) = (B). Let the auxiliary unknown (Y) = (Tt)-1(X), this becomes (Rt)(Y) = (B), which is "simpler", because of the row Euclidean reduction of (R) into a matrix which is at worse triangular.
Solving (A)(X) = (B) is then equivallent to solve  (Rt)(Y) = (B)
Then (Y) = (Tt)-1(X) gives the general solution of (A)(X) = (B) :  (X) = (Tt)(Y)

After this a little abstract proof of the method, just some example :
To solve the system
5x + 6y + 8z - 2t = 3
6x - 10y + 7z - 8t = 11
That is in matrix form :

```(5,  6,  8, -2)×(x) = (3 )
(6, -10, 7, -8) (y)   (11)
(z)
(t)```
Perform the row Euclidean reduction of (At|I) :
```(5   6  | 1  0  0  0)
(6  -10 | 0  1  0  0)
(8   7  | 0  0  1  0)
(-2 -8  | 0  0  0  1) The lowest absolute value is 2

(1  -10 | 1  0  0  2)
(0  -34 | 0  1  0  3)
(0  -25 | 0  0  1  4)
(2   8  | 0  0  0 -1) We may also change the signs of a row

(1  -10 | 1  0  0  2)
(0  -34 | 0  1  0  3)
(0  -25 | 0  0  1  4)
(0   28 |-2  0  0 -5)```
We thus get the GCD(5,6,8,-2) = 1. Let's go on with the row reduction from the next columns in At,
choosing as min only elements after the 2nd row, the min is then |-25| and not |-10| :
```(1  -10 | 1  0  0  2)
(0  -9  | 0  1 -1 -1)
(0  25  | 0  0 -1 -4)
(0   3  |-2  0  1 -1) New min = 3

(1 -1  | -5  0  3 -1) We may also reduce this row to get a simpler R (diagonal)
(0  0  | -6  1  2 -4)
(0  1  | 16  0 -9  4) New min = 1
(0  3  | -2  0  1 -1)

(1  0 | 11  0 -6  3 )
(0  0 | -6  1  2 -4 )
(0  1 | 16  0 -9  4 ) Or by exchanging rows :
(0  0 | -50 0 28 -13)

(1  0 | 11  0 -6  3 )
(0  1 | 16  0 -9  4 )
(0  0 | -6  1  2 -4 )
(0  0 | -50 0 28 -13)

(R) = (1  0) and (T) = (11  0 -6  3 )
(0  1)           (16  0 -9  4 )
(0  0)           (-6  1  2 -4 )
(0  0)           (-50 0 28 -13)```
We have then to solve (Rt)(Y) = (B) :
```(1 0 0 0)×(y1) = (3 )
(0 1 0 0) (y2)   (11)
(y3)
(y4)```
Which obviously gives y1 = 3 and y2 = 11, because R is "diagonal", if we can say that of a non square matrix.
It is during this step that we know if there exist integer solutions.
y3 and y4 are undefined, so we keep them as y3 = k, y4 = m. And we get then
```(Y) = (3 )
(11)
(k )
(m )```
and (X) = (Tt)(Y) :
```(x) = (11 16 -6 -50)×(3 ) = (209 - 6k - 50m )
(y)   (0   0  1  0 ) (11)   (0  + k         )
(z)   (-6 -9  2  28) (k )   (-117 + 2k + 28m)
(t)   (3   4 -4 -13) (m )   (53 - 4k - 13m  )```

Of course, all these "filthy" calculations upon matrices are to be done by a program,
as soon as the number of unknowns or of equations become a little large.
As in the above examples, the "simplest" writing for the solutions is not obvious.
Here, we may let m = 4 + m' :

 x = 9 - 6k - 50m'   y = k   z = -5 + 2k + 28m'   t = 1 - 4k - 13m'
With k and m' running independantly in Z.