Equation de Pell x² - Dy² = k

Cette équation n'a d'intérêt que si on cherche des solutions en nombre entiers.

Considérations générales

Si D est négatif (x² + Ay² = k) le nombre de solutions s'il y en a est fini.
De même si D est un carré x² - A²y² = k s'écrit (x + Ay)(x - Ay) = K où x + Ay et x - Ay sont à choisir parmi les diviseurs de K, en nombre fini.

On ne s'intéressera ici qu'au cas D positif non carré et tout d'abord au cas k = 1 :
x² - Dy² = 1 (équation de Pell proprement dite).

On a alors les résultats suivants :

L'équation x² - Dy² = 1 possède la solution triviale x = 1, y = 0
L'équation x² - Dy² = 1 possède toujours une infinité de solutions non triviales.
Elles peuvent être ordonnées par ordre croissant de x ou de y indifféremment
Pour deux solutions x₁ < x₂ <=> y₁ < y₂
La plus petite solution non triviale est appelée solution fondamentale et notée (x₀, y₀)
Si (x₁, y₁) et (x₂, y₂) sont deux solutions, alors (x, y) = (x₁, y₁)o(x₂, y₂) donnée par :
x = x₁x₂ + Dy₁y₂ et y = x₁y₂ + x₂y₁ aussi, en vertu de l'identité de Brahmagoupta
(X² + AY²)(X'² + AY'²) = (XX' -/+ AYY')² + A(XY' +/- X'Y)²
Toutes les solutions se déduisent de la solution fondamentale en appliquant la formule précédente itérativement ce qui donne :
x_n+1 = x₀x_n + Dy₀y_n, y_n+1 = y₀x_n + x₀y_n
On peut en tirer les formules de récurence donnant indépendemment les suites x_n et y_n :
x_n+2 = 2x₀x_n+1-x_n y_n+2 = 2x₀y_n+1-y_n
On peut aussi écrire : x_n + y_n√D = (x₀ + y₀√ D)ⁿ⁺¹ ou la forme conjuguée en remplaçant + par -
Cela donne les formules donnant directement x_n et y_n en fonction de n. Ces formules sont généralement de peu d'intérêt immédiat à cause des termes en √D irrationnels qui donnent pourtant des x_n et y_n entiers !

Exemple x² - 8y² = 1

Il n'est ici pas trop difficile de trouver (par divination) une solution non triviale x = 3, y = 1. Cette solution est la solution fondamentale car toute solution plus petite aurait y < 1 et serait donc la solution triviale (1, 0).
En appliquant les formules de récurrence, on obtient les solutions :
(3,1) (17,6) (99,35) (577,204) ...

Solution fondamentale

Dans des cas moins évidents, il reste à trouver la solution fondamentale.
Mentionnons la perfidie de Fermat quand il soumit à Frenicle l'équation x² - 61y² = 1.
En effet la solution fondamentale de cette équation est x₀ = 1766319049 y₀ = 226153980. (Rappelons que la solution fondamentale est la plus petite solution autre que x = 0, y = 1 ! )
La taille de cette solution est sans commune mesure avec la valeur de D, il n'y a pas de règle simple reliant ces deux grandeurs.
Par exemple l'équation voisine x² - 63y² = 1 possède la solution fondamentale évidente x = 8, y = 1.
x² - 13y² = 1 nécessite déja une méthode efficace de recherche de la solution fondamentale.
Heureusement il en existe une que nous allons donner sous forme d'un algorithme :

Posons U₀ = 0, V₀ = 1, a₀ = E[√D], P₀ = a₀, Q₀ = 1 et par convention P_-1 = 1 et Q_-1 = 0
avec la notation E[x] = partie entière de x. Puis itérons les formules :

U_n+1 = a_nV_n - U_n, V_n+1 = (D - U_n+1²)/V_n, a_n+1 = E[ (a₀ + U_n+1)/V_n+1 ],
(P_n+1, Q_n+1) = a_n+1(P_n, Q_n) + (P_n-1, Q_n-1)

En s'arrêtant quand a_k > a₀ (ou V_k = 1).
Alors (P_k-1, Q_k-1) est solution de l'équation x² - Dy² = (-1)^k.

Si k est pair c'est la solution cherchée.
Sinon, la solution de x² - Dy² = 1 s'obtient à partir de la solution (P, Q) de x² - Dy² = -1 par :
(P² + DQ², 2PQ)

Cet algorithme est un cas particulier de l'algorithme PQa donnant le développement en fraction continue d'un nombre quadratique (u + √D)/v avec ici u = 0 et v = 1.

Appliquons ces formules à la résolution de x² - 13y² = 1 :

U₀ = 0	V₀ = 1	a₀ = 3	P₀ = 3	Q₀ = 1
U₁ = 3	V₁ = (13 - 9)/1 = 4	a₁ = E[(3 + 3)/4] = 1	P₁ = 4	Q₁ = 1
U₂ = 4 - 3 = 1	V₂ = (13 - 1)/4 = 3	a₂ = E[(3 + 1)/3] = 1	P₂ = 7	Q₂ = 2
U₃ = 3 - 1 = 2	V₃ = (13 - 4)/3 = 3	a₃ = E[(3 + 2)/3] = 1	P₃ = 11	Q₃ = 3
U₄ = 3 - 2 = 1	V₄ = (13 - 1)/1 = 4	a₄ = E[(3 + 1)/4] = 1	P₄ = 18	Q₄ = 5
U₅ = 4 - 1 = 3	V₅ = (13 - 9)/4 = 1	a₅ = E[(3 + 3)/1] = 6 > 3 => stop

5 étant impair, (18, 5) est solution de x² - 13y² = -1.
La solution fondamentale de x² - 13y² = 1 est (18² + 13×5², 2×18×5) = (649, 180) ce qui n'était pas évident à trouver sans disposer de cet algorithme.

Passons maintenant au cas général où k ≠ 1 et tout d'abord

Equation x² - Dy² = -1

Cette équation peut ne pas avoir de solution, mais si elle en a, elle en a une infinité.
Là encore, elles peuvent être rangées en ordre croissant et la plus petite est appelée "solution fondamentale" ou "primitive" et notée (x₀, y₀)
Toutes les solutions sont données par x_n + y_n√D = (x₀ + y₀√ D)²ⁿ⁺¹

ou la forme conjuguée en remplaçant + par -

Nota : les solutions de l'équation x² - Dy² = 1 sont obtenues à partir de la solution fondamentale
(x₀, y₀) de x² - Dy² = -1 par : x_n + y_n√D = (x₀ + y₀√D)²ⁿ
La solution fondamentale de x² - Dy²=1 est ainsi (x₀² + Dy₀², 2x₀y₀).
On a aussi les relations de récurence :

x_n+1 = (x₀² + Dy₀²)x_n + 2Dx₀y₀y_n, y_n+1 = 2x₀y₀x_n + (x₀² + Dy₀²)y_n

Ainsi que des formules en x_n+2 = ax_n+1 + bx_n que je vous laisse chercher à titre d'exercice.

Quant à la recherche de la solution fondamentale, on se reportera à l'algorithme précédent en notant que si k est pair, c'est que x² - Dy² = -1 n'a pas de solution. Il n'y a pas en général de critère simple pour le savoir à priori. Citons toutefois le cas ou D est un nombre premier p :
Une condition nécessaire et suffisante pour que l'équation x² - py² = -1 soit soluble est que p = 4k + 1.

Exemple :

Soit l'équation x² - 2y² = -1
x = 1, y = 1 est une solution évidente, qui est fondamentale puisqu'une solution plus petite aurait x et y < 1.
La solution fondamentale de x² - 2y² = +1 est alors x = 1² + 2*1² = 3 et y = 1² + 1² = 2 .
Et toutes les solutions sont données par x_n+1 = 3x_n + 4y_n, y_n+1 = 2x_n + 3y_n, c'est à dire :
(1,1) (7,5) (41,29) (239,169) (1393,985) ...

Equation de Pell généralisée

Soit x² - Dy² = K avec |K| ≠ 1
Là encore, cette équation peut ne pas avoir de solution mais si elle en a, elle en a une infinité. Par contre, il n'y a plus de solution fondamentale mais un ensemble de solutions primitives
Continuons à appeler (x₀, y₀) la solution fondamentale de x² - Dy² = 1.
Alors toutes les solutions de x² - Dy² = K sont obtenues à partir de l'ensemble des solutions primitives indépendantes (ξ, η) par l'ensemble des formules

x_n + y_n√D = (x₀ + y₀√D)ⁿ(ξ + η√D) ainsi que les formules de récurence correspondantes.

Le problème de trouver les solutions primitives indépendantes n'est toutefois pas évident. Une première méthode est de balayer systématiquement les valeurs de x et y, la théorie se contentant de fournir une borne supérieure à ces essais (rappelons que l'équation peut ne pas avoir de solution). C'est la méthode la plus efficace si la borne supérieure n'est pas trop grande.

Tout d'abord explicitons ce que sont des solutions indépendantes et plus précisément celles qui ne le sont pas et qui sont alors dites "équivalentes".
Deux solutions (x,y) et (x',y') sont équivalentes s'il existe une solution (X,Y) de x² - Dy² = 1 telle que :

x + y√D = (X + Y√D)(x' + y'√D)

Ce critère étant peu pratique, ajoutons :

Pour que deux solutions (x,y) et (x',y') de x² - Dy² = K soient équivalentes, il faut et il suffit que :
xx' - Dyy' et x'y - xy' soient multiples de K.

Reste à donner les bornes des solutions primitives :

si K > 0, les solutions primitives sont bornées par
0 < x ≤ √( (x₀ + 1)K/2 ) ou 0 ≤ y ≤ y₀√ ( K/2(x₀ + 1) )

si K < 0, les solutions primitives sont bornées par
0 ≤ x ≤ √( (x₀ - 1)|K|/2 ) ou 0 < y ≤ y₀√( |K|/2(x₀ - 1) )

S'il n'y a pas de solution entre ces bornes, c'est qu'il n'y a pas de solution du tout.

On pourrait être tenté de chercher un critère pour que x² - Dy² = K soit soluble en écrivant la congruence
x² ≡ K mod (D), c'est à dire que K doit être un résidu quadratique de D.
La seule chose que l'on peut dire alors est que si K n'est pas un résidu quadratique de D, l'équation de Pell est insoluble. Par exemple l'équation x² - 97y² = 51 est insoluble car le symbole de Legendre (51/97) = -1 montre que 51 n'est pas un résidu quadratique de 97.
Mais cette condition nécessaire n'est pas suffisante...

Exemple :

Soit l'équation x² - 7y² = 9
x² - 7y² = 1 a pour solution fondamentale x₀ = 8, y₀ = 3 (exercice de révision !)
il faut chercher 0 ≤ y ≤ 3√(9/18) soit 0 ≤ y ≤ 2.
y = 0 donne x = ±3
y = 1 donne x = ±4
y = 2 donne x² = 37 impossible.

Il y a donc 4 solutions primitives, mais qui ne sont pas toutes indépendantes
(±3, 0) sont équivalentes, de façon générale (x,y) ~ (-x,-y)
(4,±1) sont indépendantes : 16 + 7 = 26 premier avec 9.
(3,0) et (4,1) sont indépendantes
(3,0) et (4,-1) sont indépendantes

Finalement trois solutions primitives indépendantes : (3, 0), (4, 1) et (4,-1)
Si on s'intéresse aux x,y > 0, (4,-1) ~ (11,4), plus petite solution >0 de cette classe
Les solutions générales sont formées des trois familles :

x_n + y_n√7 = (8 + 3√7)ⁿ(3 + 0√7)
x_n + y_n√7 = (8 + 3√7)ⁿ(4 + √7)
x_n + y_n√7 = (8 + 3√7)ⁿ(4 - √7)

Ou les formules de récurence :

x_n+1 = 8x_n + 21y_n, y_n+1 = 3x_n + 8y_n, en partant de (x, y) = (3,0), (4,1) et (4,-1)

soit les trois familles :
(3,0) (24,9) (381,144) ...
(11,4) (172,65) (2741,1036) ...
(4,1) (53,20) (844,319) ...

La recherche par balayage des solutions devient fastidieuse si la borne est trop grande. Dans ce cas on dispose de deux méthodes plus efficaces.

Cas où K² < D

Si on reprend l'algorithme de recherche des solutions de x² - Dy² = ±1, on a l'invariant de boucle

P_i² - DQ_i² = (-1)ⁱ⁺¹V_i+1V₀

Ceci permet de trouver dans le même algorithme les solutions primitives de x² - Dy² = K
Notons tout d'abord que si (x,y) est une solution de x² - Dy² = K, (mx,my) est une solution de x² - Dy² = m²K.
Si donc on obtient un (-1)ⁱ⁺¹V_i+1 = K/m², alors x = mP_i, y = mQ_i est solution de x² - Dy² = K
Si K² < D, ceci donne toutes les solutions primitives.

Exemple x² - 19y² = -3

U	V	a	P	Q
0	1	4	4	1	1
4	3	2	9	2	-3	4	1	est solution de x² - 19y² = -3
2	5	1	13	3	5
3	2	3	48	11	-2
3	5	1	61	14	5
2	3	2	170	39	-3	61	14	est solution de x² - 19y² = -3
4	1	8			1	170	39	est solution de x² - 19y² = 1

Les solutions primitives de x² - 19y² = -3 sont (4,1) et (61,14)

Notas :
1) l'équation x² - 19y² = K n'a donc pas de solutions pour K = -1, +2, +3, et -4.
2) x² - 19y² = 5 admet les solutions primitives (9,2) et (48,11) mais on ne sait pas si ce sont les seules car 5² > 19.

Algorithme LMM

L'algorithme de Lagrange-Matthews-Mollin permet de s'en sortir si K² > D

Lister les f  tels que f² divise K, m=K/f²
Pour chaque f, trouver tous les z : -|m|/2 < z≤|m|/2 et z² ≡ D modulo |m|
Pour chaque z appliquer l'algorithme PQa avec U₀ = z, V₀ = |m| jusqu'à la fin de la période ou un V_i = ±1
  - Si V_i = 1, alors (P_i-1,Q_i-1) est solution de x² - Dy² = m, ajouter cette solution à la liste.
  - Si V_i = -1, alors (P_i-1,Q_i-1) est solution de x² - Dy² = -m. Si x² - Dy² = -1 a une solution (r,s) (primitive)
    alors (P_i-1,Q_i-1)o(r,s) est solution de x² - Dy²=m. Ajouter (P_i-1,Q_i-1)o(r,s) = (P_i-1r + Q_i-1sD, P_i-1s + Q_i-1r)

Exemple :

x² - 13y² = 12
f = 1, m = 12, z² ≡ 13 modulo 12, z = ±1 et z = ±5
   z = -1 donne la solution (83,23)
   z = 1 donne la solution (47,13)
   z = -5 donne la solution (5,1)
   z = 5 donne la solution (-5,1)
f = 2, m = 3, z² ≡ 13 modulo 3, z = ±1
   z = -1 donne la solution (8,2)
   z = 1 donne la solution (512,142)

toutes ces solutions étant indépendantes et donnant donc 6 familles de solutions par combinaison avec la solution (649,180)ⁿ de x² - 13y² = 1

Equation ax² + bxy + cy² = K

L'algorithme LMM peut être étendu au cas plus général des équations ax² + bxy + cy² = K avec D = b² - 4ac > 0 non carré.