Equation de Diophante ax + by = c

Bien entendu avec a, b, c, x, y en nombres entiers ∈ Z.

Tout d'abord si p est un diviseur de a et b, il doit diviser c sinon l'équation est insoluble. On peut alors tout diviser par p.
On ne s'intéressera dorénavant qu'aux équations où a et b sont premiers entre eux, c'est à dire sans diviseurs communs (autre que 1).
Soit ax₀ + by₀ = c une solution particulière. Alors si ax + by = c est une solution quelconque, par soustraction membre à membre, il vient :
a(x - x₀) + b(y - y₀) = 0. Comme a et b sont premiers entre eux, a divise y - y₀, soit y - y₀ = ka ou encore :
y = y₀ + ka et x = x₀ - kb

Soit maintenant l'équation ax + by = 1, toujours avec a et b premiers entre eux. Le théorème de Bézout affirme qu'une condition nécessaire et suffisante pour que deux nombres a et b soient premiers entre eux est qu'il existe deux nombres u et v (dans Z) tels que au + bv = 1. Ceci prouve que l'équation ax + by = c, PGCD(a, b) = 1 a toujours des solutions. Si (u, v) est solution de ax + by = 1, alors (cu, cv) est solution de ax + by = c. En résumé :

L'équation ax + by = c est soluble si et seulement si p = PGCD(a,b) divise c En posant alors a = pa', b = pb' c = pc', si x0,y0 est une solution particulière de  a'x + b'y = 1, alors toutes les solutions de ax + by = c sont données par :      x = c'x0 - kb'     y = c'y0 + ka' où k parcourt Z

Le problème est de trouver une solution particulière.
Si on peut dans les cas simples opérer par tâtonnements, une équation comme 97x + 29y = 4000 nécessite une méthode plus efficace.

Changements de variables

Prenons comme exemple l'équation 97x + 29y = 4000.

Comme 97 = 29×3 + 10 et 4000 = 29×137+27, l'équation 97x + 29y = 4000 peut s'écrire :
10x + 29(y + 3x - 137) = 27 soit en posant z = y + 3x - 137 : 10x + 29z = 27, qui est une équation plus simple.

On peut recommencer en écrivant : 10(x + 2z - 2) + 9z = 7 où encore 10t + 9z = 7 puis :
t + 9(z + t) = 7 soit t + 9u = 7 qui possède la solution évidente t = 7, u = 0

En remontant les changements de variables, il vient : z = u - t = -7, puis x = t + 2 - 2z = 23 et enfin
y = z + 137 - 3x = 61, d'où une solution particulière x₀ = 23, y₀=61,
et la solution générale x = 23 - 29k, y = 61 + 97k.

Congruences

Une méthode séduisante consiste à écrire ax + by = 1 comme une congruence : ax ≡ 1 mod (b).
Comme a et b premiers entre eux, le théorème d'Euler-Fermat affirme que a^φ(b) ≡ 1 mod (b)
avec φ(b) l'indicateur d'Euler de b.
En multipliant les deux membres de la congruence ax ≡ 1 mod (b) par a^φ(b)-1 on obtient...
la solution générale de ax + by = c :

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

Reste qu'il faut calculer φ(b) (et pour les incrédules à démontrer le théorème d'Euler-Fermat)
Une méthode simple pour calculer φ(b) à partir de la décomposition en facteurs premiers de b :
Si b = ∏p_i^ai, alors φ(b) = b∏(1 - 1/p_i) = ∏(p_i^ai - p_i^ai-1). En particulier si p est premier, φ(p) = p - 1.

Pour notre équation : x ≡ 4000.97^φ(29)-1 mod (29), soit 4000×97²⁷
Heureusement, point n'est besoin de calculer des nombres astronomiques car tous les calculs se font modulo 29, les nombres mis en oeuvre ne dépasseront pas 29² = 841.
4000 ≡ 27 mod (29) et 97 ≡ 10 mod (29) permettent déjà de simplifier les calculs en :
4000×97²⁷ ≡ 27×10²⁷ mod (29), ouf !
Reste à calculer efficacement 10²⁷ mod (29). On utilise une méthode connue sous le nom "d'exponentiation indienne"

Algorithme de calcul de V = A^B

Initialiser V = 1 Tant que B ≥ 1     - si B est impair, multiplier V par A  et retrancher 1 à B     - sinon, élever A au carré et diviser B par 2

Ou une variante en binaire :

Initialiser V = 1 Pour chaque bit de l'exposant B, en commençant par les poids forts :      - élever V au carré     - si ce bit vaut 1, multiplier V par A

Ce qui donne ici 27 = 11011₍₂₎ et 10²⁷ = (((1²×10)²×10)²)²×10)²×10
Comme on calcule à chaque multiplication ou élévation au carré modulo 29 :
10² = 100 ≡ 13 mod (29)
10³ ≡ 13×10 = 130 ≡ 14 mod (29)
10⁶ ≡ 14² =196 ≡ 22 mod (29)
10¹² ≡ 22² =484 ≡ 20 mod (29)
10¹³ ≡ 20×10 = 200 ≡ 26 mod (29)
10²⁶ ≡ 26² = 676 ≡ 9 mod (29)
10²⁷ ≡ 9×10 = 90 ≡ 3 mod (29)
et enfin x ≡ 27×3 = 81 ≡ 23 mod (29)

Algorithme d'Euclide

Les calculs précédents, même s'ils partent d'une formule élégante, ne sont toutefois pas très efficaces, ne serait-ce que parce qu'ils nécessitent de décomposer b en facteurs premiers (pour calculer φ(b)).
L'algorithme d'Euclide, qui permet d'obtenir le PGCD de a et b, fournit comme sous produit la représentation en fraction continue de a/b. L'étude des fractions continues conduit alors à une variante de l'algorithme d'Euclide qui nous donne en prime la solution de ax + by = 1.
Tout d'abord notons que par changement de variables (échange de x et y et changement de signes) on peut imposer a > b > 0. Dans ces conditions :

Posons : r0 = a, r1 = b (a > b > 0) P0 = 0, P1 = 1 et Q0 = 1, Q1 = 0 Puis répéter depuis i = 2     - Division entière de ri-2 par ri-1 : ri-2 = ri-1qi + ri      - Calculer Pi = -qiPi-1 + Pi-2 et Qi = -qiQi-1 + Qi-2      - jusqu'à rn = 0 p = rn-1 est le PGCD de a et b (algorithme d'Euclide)  En posant a' = a/p et b' = b/p, x0 = Qn-1, y0 = Pn-1 est solution de a'x + b'y = 1

Démonstration de cet algorithme.

Avec 97x + 29y = 4000 cela donne :
97 = 29×3 + 10, r₂ = 10, q₂ = 3, P₂ = -3, Q₂ = 1
29 = 10×2 + 9, r₃ = 9, q₃ = 2, P₃ = 7, Q₃ = -2
10 = 9×1 + 1, r₄ = 1, q₄ = 1, P₄ = -10, Q₄ = 3
9 = 1×9 + 0, r₅ = 0, c'est fini et x = 3 y = -10 est solution de 97x + 29y = 1
En multipliant par c = 4000, qui est bien un multiple de p = r₄ = 1, on obtient x = 12000 + 29k, y = -40000 - 97k
Si on pose m = k + [12000/29] = k + 413, on retrouve x = 23 + 29m, y = 61 - 97m comme ci-dessus.

Script de résolution de ax + by = c avec cet algorithme.

Equations à plus de 2 inconnues

Prenons comme exemple 15x + 21y + 35z = 200.
Pour résoudre une telle équation, il suffit de considérer toutes les inconnues sauf deux comme des paramètres et de les faire passer dans le second membre.
Ce qui donne : 15x + 21y = 200 - 35z = c et on est amené à résoudre 15x + 21y = c en fonction de c.
Petit problème bien sûr car 15 et 21 ne sont pas premiers entre eux. c doit donc être un multiple du PGCD 3 de 15 et 21.
C'est à dire que 200 - 35z ≡ 0 mod (3) ou encore z ≡ 1 mod(3). Ce qui s'écrit z = 3m + 1. La véritable inconnue est m et non pas z.
Notre équation s'écrit 15x + 21y = 165 - 3×35m et en divisant tout par 3 : 5x + 7y = 55 - 35m.
Comme x₀ = 3, y₀ = -2 est solution de 5x + 7y = 1 (exercice de révision...), la solution générale de cette équation s'écrit
x = 3(55 - 35m) + 7k, y = -2(55 - 35m) - 5k et finalement

x = 165 - 105m - 7k y = -110 + 70m + 5k z = 3m + 1 k et m parcourant indépendamment Z

Donnons les formules générales pour l'équation ax + by + cz = D

Soit p = PGCD(a,b), a' = a/p, b' = b/p Soient u0 et v0 une solution particulière de a'u + b'v = c  z0, t0 une solution particulière de cz + pt = D x0, y0 une solution particulière de a'x + b'y = t0 La solution générale de ax + by + cz = D est : x = x0 + b'k - u0m y = y0 - a'k - v0m z = z0 + pm où k et m parcourent indépendamment Z

Un Script utilisant cette méthode.

L'application directe de ces formules à 15x + 21y + 35z = 200 donne en fait :
p = 3, a' = 5, b' = 7 puis 5u + 7v = 35 et u₀ = 0 v₀ = 5 (par exemple)
35z + 3t = 200 z₀ = 1, t₀ = 55
5x + 7y = 55 : x₀ = 11, y₀ = 0 et finalement :

x = 11 + 7k y = 0 - 5k - 5m z = 1 + 3m

qui semble différente de la solution trouvée ci-dessus. Mais le changement de variable k = 22 - k' - 15m remet tout en place.

Cette méthode est un peu pénible car elle nécessite de résoudre successivement plusieurs équations à deux inconnues, donc d'exécuter plusieurs fois l'algorithme d'Euclide.
Du coup, au dela de 3 inconnues elle devient impraticable.
On peut alors efficacement résoudre cette équation en considérant un système (hum) de une équation (!) à n inconnues...

Généralisation de l'algorithme d'Euclide

Une première remarque est que le PGCD de n nombres peut s'obtenir en "téléscopant" les algorithmes d'Euclide.
En fait l'algorithme d'Euclide (juste le PGCD) peut se décrire par des combinaisons linéaires successives des valeurs :
Par exemple soit à calculer le PGCD(15, 21, 35)
Le plus petit des 3 nombres est 15, on peut donc effectuer une division Euclidienne de tous les autres par 15 en ne gardant que les restes.
Soit PGCD(15, 21, 35) = PGCD(15, 21 - 15 , 35 - 2×15) = PGCD(15, 6, 5)
On recommence alors avec le plus petit 5 : PGCD(15, 6, 5) = PGCD(15 - 3×5, 6 - 5, 5) = PGCD(0, 1, 5)
Et une dernière fois (avec le plus petit non nul) : PGCD(0, 1, 5) = PGCD(0, 1, 5 - 5×1) = PGCD(0, 1, 0) = 1
Cette méthode se termine forcément puisque à chaque étape, les restes sont tous < la plus petite valeur non nulle.
Celle-ci est donc strictement décroissante.
L'algorithme d'Euclide étendu (qui donne en plus la relation de Bézout ax + by = 1) fait en fait "subir" simultanément cette réduction aux lignes de la matrice unité, c'est à dire que l'on considère la matrice obtenue en juxtaposant la colonne a, b et une matrice unité :

(a | 1 0)
(b | 0 1)

Généralisant cette réduction à notre équation 15x + 21y + 35z = 200, on effectue la "réduction Euclidienne" de la matrice :

(15 | 1  0  0)
(21 | 0  1  0)
(35 | 0  0  1) ce qui donne :

(15 | 1  0  0)
(6  |-1  1  0) = moins [21/15]=1 fois la 1ère ligne
(5  |-2  0  1) = moins [35/15]=2 fois la 1ère ligne

(0  | 7  0 -3) = moins [15/5]=3 fois la 3ème ligne
(1  | 1  1 -1) = moins [6/5]=1 fois la 3ème ligne
(5  |-2  0  1)

et finalement :

(0  | 7  0 -3) 
(1  | 1  1 -1) 
(0  |-7 -5  6) = moins [5/5]=5 fois la 2ème ligne

La seule ligne avec une valeur non nulle en 1ère colonne (= PGCD = 1) nous donne directement la relation de Bézout : 15×1 + 21×1 + 35×(-1) = PGCD = 1
C'est à dire une solution particulière de notre équation : x = 200×1, y = 200×1, z = 200×(-1)

Les lignes commençant par 0 de la matrice nous donnent des solutions particulières de l'équation 15x + 21y + 35z = 0
soit (x = 7, y = 0, z = -3) et (x = -7, y = -5, z = 6), qui sont indépendantes, et qu'on peut d'ailleurs remplacer par deux combinaisons linéaires indépendantes, par exemple u = (7, -5, 0) et v = (0, -5, 3)
On peut donc ajouter une combinaison linéaire quelconque de u et v : ku + mv.

Ceci donne les solutions générales de notre équation :

x = 200 + 7k + 0m   y = 200 - 5k - 5m   z = -200 + 0k + 3m

Ce qui à un décalage près (m = 67 + m', k = -27 + k') redonne les valeurs trouvées précédemment.

Cette méthode est applicable à une équation avec un nombre quelconque d'inconnues.

Systèmes d'équations

Et même plus : on peut généraliser à un système d'équations, sans avoir à effectuer une réduction préalable du système en une seule équation.

Considérons donc notre système d'équations sous la forme matricielle (A)(X) = (B).
Effectuons les opérations précédentes de réduction Euclidienne sur la matrice (A^t|I), concaténation de la transposée de (A) et d'une matrice unité I.
On peut alors obtenir une matrice (R|T) où T est la matrice de passage : (R|T) = (T)(A^t|I)
(T) est inversible, puisqu'on n'a effectué que des combinaisons linéaires de lignes.
Ceci donne aussi (T)(A^t) = (R), et par conséquent :
(A)(T^t) = (R^t), et finalement puisque T est inversible : (A) = (R^t)(T^t)^-1
Le déterminant de T étant ±1 (comme produit de matrices de déterminant ±1, matrices unimodulaires), l'inverse (T^t)^-1 est une matrice de nombres entiers.
L'équation (A)(X) = (B) s'écrit (R^t)(T^t)^-1(X) = (B), et en posant la variable auxilliaire (Y) = (T^t)^-1(X), elle devient (R^t)(Y) = (B), qui est plus "simple",à cause de la réduction Euclidienne de (R) en une matrice finalement au pire triangulaire.
Résoudre (A)(X) = (B) est équivallent à résoudre  (R^t)(Y) = (B) 
puis (Y) = (T^t)^-1(X) nous donne la solution générale de (A)(X) = (B) :  (X) = (T^t)(Y) 

Après cette justification un peu abstraite de la méthode, un petit exemple :
Soit le système
5x + 6y + 8z - 2t = 3
6x - 10y + 7z - 8t = 11
ou encore sous forme matricielle :

(5,  6,  8, -2)×(x) = (3 )
(6, -10, 7, -8) (y)   (11)
                (z)
                (t)

Effectuons la réduction Euclidienne de (A^t|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) Le mini en valeur absolue est 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) On peut aussi changer de signe une ligne

(1  -10 | 1  0  0  2)
(0  -34 | 0  1  0  3)
(0  -25 | 0  0  1  4)
(0   28 |-2  0  0 -5)

On obtient déja le PGCD(5,6,8,-2) = 1. Poursuivons sur les colonnes suivantes de A^t,
en ne choisissant comme pivots que des éléments à partir de la deuxième ligne, le min est ici |-25| et non |-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) Nouveau mini = 3

(1 -1  | -5  0  3 -1) On réduit aussi cette ligne pour obtenir une matrice R plus simple (diagonale)
(0  0  | -6  1  2 -4)
(0  1  | 16  0 -9  4) Nouveau mini = 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 ) Ou en échangeant les lignes :
(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)  et (T) = (11  0 -6  3 )
      (0  1)           (16  0 -9  4 )
      (0  0)           (-6  1  2 -4 )
      (0  0)           (-50 0 28 -13)

Il faut alors résoudre (R^t)(Y) = (B) :

(1 0 0 0)×(y1) = (3 )
(0 1 0 0) (y2)   (11)
          (y3)
          (y4)

Qui donne de façon triviale y1 = 3 et y2 = 11, car R est "diagonale", si l'on peut dire en parlant d'une matrice non carrée.
C'est à cette étape que l'on détermine si le système a des solutions ou non en nombres entiers.
y3 et y4 sont indéterminées, et on les laisse sous la forme de y3 = k, y4 = m. On obtient ainsi

(Y) = (3 )
      (11)
      (k )
      (m )

et (X) = (T^t)(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  )

Bien entendu, tous ces calculs "sordides" sur des matrices sont à faire par un programme,
dès que le nombre de variables ou d'équations devient un tant soit peu conséquent.
Comme pour les exemples précédents, l'écriture la plus "simple" des solutions ne s'obtient pas de façon triviale.
Ici on peut poser m = 4 + m' :

x = 9 - 6k - 50m'   y = k   z = -5 + 2k + 28m'   t = 1 - 4k - 13m'

Où k et m' parcourent Z indépendamment.