Solution de systèmes polynomiaux par homotopie

Éric Schost, LIX, École polytechnique

`Eric.Schost[at]polytechnique.fr`

Vue d'ensemble

Le but de ce projet est de “résoudre” des systèmes de n équations polynomiales de degré 2 à n inconnues (un tel système a généralement 2ⁿ solutions). Il s'agit de mettre en place un cas particulier simplifié de ce qu'on appelle méthodes par homotopie.

Résoudre un système F=F₁,...,F_n dans k[X₁,...,X_n] (k est un corps pour l'instant quelconque) reviendra à calculer des polynômes à une variable P,V₂,...,V_n tels que l'ensemble des solutions de F=0 soit décrit par

⎧
⎪
⎨
⎪
⎩

P(X₁)	=	0
X₂	=	V₂(X₁)
	⋮
X_n	=	V_n(X₁);

on conçoit que sous cette forme, on accède aisément à la plupart des caractéristiques de l'ensemble des solutions (par exemple, trouver les racines de P dans k nous donne les solutions à coordonnées dans k, etc ...). Tous les systèmes n'admettent pas des solutions de cette forme. Cependant, ce sera le cas pour les systèmes qui nous intéresseront : on choisira pour F₁,...,F_n des polynômes “génériques” (tous les coefficients sont aléatoires) de degré 2.

Pour calculer la résolution, on commence par définir un système G de référence. L'objectif est de pouvoir résoudre G “à vue” ; on prend donc

G₁ = X₁²−X₂, G₂ = X₂² − X₃, ..., G_n−1 = X_n−1²−X_n, G_n = X_n²−1.

Il est très facile de résoudre ce système sur le modèle ci-dessus ; on obtient

⎧
⎪
⎨
⎪
⎩

X₁²ⁿ	=	1
X₂	=	X₁²
	⋮
X_n	=	X₁²ⁿ⁻¹.

Pour utiliser ce système de référence, on définit une nouvelle variable ε, et on transforme le système F, en introduisant F^ε donné par

F₁^ε= ε F₁ + (1− ε) G₁, ... , F_n^ε= ε F_n + (1−ε) G_n,

En faisant ε = 1, on retrouve F ; en faisant ε =0, on retrouve G. Ceci mène à l'algorithme suivant :

On résout le système en ε=0 : on a vu que c'est immédiat.
On en déduit une solution “paramétrée par ε” du système F^ε. Cette solution s'écrit

⎧
⎪
⎨
⎪
⎩

P^ε(X₁) = 0

X₂ = V₂^ε(X₁)

⋮

X_n = V_n^ε(X₁),

où P^ε,V₂^ε,...,V_n^ε sont des polynômes en une variable à coefficients dans k(ε), c'est-à-dire des fractions rationnelles en ε (voir l'exemple à la fin du sujet).
On fait ε=1 dans la solution paramétrée.

Graphiquement, on cherche à déterminer les points de droite, à partir des points de gauche, en “calculant” une description des courbes qui les relient.

C'est bien entendu la seconde étape qui pose toutes les difficultés, et sur laquelle porte l'essentiel du projet.

Calcul de la solution paramétrée

Avant de donner les détails de l'algorithme, on présente rapidement les structures et les outils dont on aura besoin.

Séries formelles sur k.

Une série formelle à coefficients dans k est une “somme formelle” de la forme Σ_{i ≥ 0} a_i εⁱ ; on note l'ensemble de ces séries sous la forme k[[ε]].

Ces séries se comportent comme des polynômes pour les opérations d'addition et de multiplication. Concrètement, on manipule uniquement des séries tronquées, c'est-à-dire avec un nombre fini de termes. On introduit donc une notion de précision : la série a=Σ_i=0ⁿ⁻¹ a_i εⁱ sera dite de précision n. Quand on ajoute ou on multiplie deux séries tronquées, on effectue l'opération comme pour des polynômes, puis on tronque le résultat à la précision courante.

Une différence importante entre séries et polynômes est la possibilité d'inverser une série formelle dont le terme constant n'est pas nul. Par exemple, l'inverse de la série 1−ε est la série Σ_{i ≥ 0} εⁱ. La formule suivante donne l'inverse de la série a = Σ_{i ≥ 0} a_i εⁱ : on pose b₀ = 1/a₀ et

b_i+1 = 2 b_i − a b_i² mod ε²ⁱ⁺¹

pour i ≥ 0. Alors b_i est l'inverse de a modulo ε²ⁱ.

Approximation de Padé.

Le second outil dont on a besoin est le calcul d'approximants de Padé. Si f=p/q est une fraction rationnelle dans k(ε), le calcul d'inverse ci-dessus permet de calculer le développement de Taylor de f en 0. L'approximation de Padé permet de faire l'inverse : si on connaît le développement de Taylor de f avec suffisamment de termes (degp +degq+1, en l'occurence), on peut retrouver p et q (à une constante près). L'algorithme est simple à implanter, c'est une variante de l'algorithme d'Euclide, il est donné au chapitre 5 de Modern Computer Algebra, J. von zur Gathen et J. Gerhard, Cambridge University Press, 1999.

Solutions en séries.

On fait désormais l'hypothèse que le polynôme X₁²ⁿ −1 a toutes ses racines dans k (on verra plus bas des situations où c'est le cas). Dans ce cas, il est possible d'expliciter toutes les solutions du système G. À partir de ces solutions, on va calculer des solutions sous forme de séries formelles en ε pour le système F^ε.

Comme pour l'inverse, on travaille itérativement. On se donne une solution initiale x⁰=(x₁⁰,...,x_n⁰) pour F^εmod ε, c'est à dire G. On peut alors calculer la suite xⁱ=(x₁ⁱ,...,x_nⁱ), où chaque x_jⁱ est une série à précision 2ⁱ, et telle que F^ε(xⁱ)=0 modulo ε²ⁱ. On vérifie facilement que la suite définie par

xⁱ⁺¹ = xⁱ − J(xⁱ)⁻¹ F^ε(xⁱ) mod ε²ⁱ⁺¹

convient. Ici, J est la matrice jacobienne du sysème F^ε par rapport à X₁,...,X_n, J(xⁱ) est son évaluation en xⁱ, et J(xⁱ)⁻¹ l'inverse de celle-ci. Remarquer qu'il s'agit ici d'une version formelle de l'opérateur de Newton, qui généralise ce qu'on a écrit pour l'inverse d'une série.

Solutions paramétrées.

On a maintenant tous les outils nécessaires. La première chose à faire est de calculer toutes les solutions de G ; c'est immédiat, en utilisant l'hypothèse qu'on connaît toutes les racines de X₁²ⁿ −1. Il y a 2ⁿ solutions.

En utilisant l'opérateur de Newton, on calcule 2ⁿ solutions en séries formelles de F^ε, en allant jusqu'à la précision 2ⁿ⁺¹+1. On va noter x⁽¹⁾,...,x⁽²ⁿ⁾ ces solutions, où x⁽ⁱ⁾=(x₁⁽ⁱ⁾,...,x_n⁽ⁱ⁾).

On peut d'ores et déjà en déduire le polynôme P^ε. On calcule pour cela le polynôme dont les racines sont les x₁⁽ⁱ⁾. Il est facile de voir que ce polynôme coïncide “presque” avec P^ε : ses coefficients sont les développements de Taylor à précision 2ⁿ⁺¹+1 de ceux de P^ε. Par ailleurs, on peut montrer que les coefficients de P^ε ont des numérateurs et dénominateurs de degrés au plus 2ⁿ. Donc pour retrouver P^ε, on calcule un approximant de Padé de chacun des coefficients, en prenant comme degré maximal pour les numérateurs et dénominateurs 2ⁿ.

Un léger problème se pose pour les polynômes V^ε : leurs coefficients (dans k(ε)) ont des degrés de l'ordre de 4ⁿ, ce qui est beaucoup trop. Pour contourner cette difficulté, on définit W_j^ε = (V_j^ε × (P^ε)') mod P^ε, qui vit aussi dans k(ε)[X₁], et on montre que ce choix est le bon : les coefficients de W_j^ε ont un degré au plus en 2ⁿ. Par ailleurs, on voit facilement que W_j^ε(X₁) vaut

i≤ 2ⁿ

x_j⁽ⁱ⁾

ℓ ≠ i

(X₁ − x_j^(ℓ)), (1)

si on calcule “à précision infinie”. On commence donc par calculer cette expression avec des coefficients dans k[[ε]] à précision 2ⁿ⁺¹+1, et on en déduit les coefficients de W_j^ε par approximation de Padé en degré 2ⁿ. On fait alors ε=0 dans W_j^ε pour obtenir un W_j, et on retrouve V_j par la formule W_j = (V_j × P') mod P (ce dernier calcul nécessite un pgcd étendu).

Racines de l'unité dans k.

Ce qui précède demande d'avoir certaines racines de l'unité dans k. Un choix naturel est alors de travailler dans un corps fini de la forme

, avec p un premier bien choisi. En effet, si p est un premier de la forme ℓ 2ⁿ +1, on peut montrer que X²ⁿ−1 a toutes ses racines dans k.

Par exemple, on pourra prendre p=29 ⋅2⁵⁷+1 et ω = 21. Alors, ω est une racine 2⁵⁷-ième de 1 dans

. Dans ce cas, pour tout n ≤ 57, on pose τ = ω²⁵⁷⁻ⁿ ; les racines de X²ⁿ−1 sont les puissances 1,τ,...,τ²ⁿ⁻¹, et elles sont distinctes.

Exemple

On choisit n=2, on travaille modulo 9001, et le système à résoudre est

F₁=8996X₂² + 66X₂X₁ + 99X₂ + 77X₁² + 54X₁ + 8940 = 0,

F₂=31X₂² + 8989X₂X₁ + 8975X₂ + 8951X₁² + 8983X₁ + 8939=0.

Modulo 9001, on les a 4 racines quatrièmes de l'unité 1, 9000,1237,7764. On en déduit que les solutions de G: {X₁²−X₂,X₂²−1} sont

( 1, 1 ), ( 9000, 1 ), ( 1237, 9000 ), ( 7764, 9000 ).

En partant de (1,1), on obtient successivement par l'opérateur de Newton les solutions de F^ε en “séries tronquées”, jusqu'à la précision 9 :

( 1 + 6670ε , 1 + 4569ε )

( 1 + 6670ε + 4674ε² + 8487ε³ , 1 + 4569ε + 4028ε² + 1959ε³ )

( 1 + 6670ε + 4674ε² + 8487ε³ + 6058ε⁴ + 4716ε⁵ + 1854ε⁶ + 1415ε⁷,

1 + 4569ε + 4028ε² + 1959ε³ + 5533ε⁴ + 4596ε⁵ + 5107ε⁶ + 1568ε⁷)

( 1 + 6670ε + 4674ε² + 8487ε³ + 6058ε⁴ + 4716ε⁵ + 1854ε⁶ + 1415ε⁷ + 8903ε⁸,

1 + 4569ε + 4028ε² + 1959ε³ + 5533ε⁴ + 4596ε⁵ + 5107ε⁶ + 1568ε⁷+ 840ε⁸)

On fait de même avec les 3 autres solutions. Ceci permet de reconstituer le polynôme ayant pour racines les coordonnées X₁ de ces solutions, à coefficients à précision 9 :

		X₁⁴ + (96ε + 2400ε² + 5521ε³ + 82ε⁴ + 3518ε⁵ + 3262ε⁶ + 8605ε⁷+ + 5763ε⁸ )X₁³
		+ (8793ε + 7463ε² + 7042ε³ + 6137ε⁴ + 789ε⁵ + 3829ε⁶ + 4951ε⁷ + 5255ε⁸)X₁²
		+ (114ε + 5848ε² + 3466ε³ + 8386ε⁴ + 2077ε⁵ + 4952ε⁶ + 364ε⁷ + 4994ε⁸)X₁
		+ 9000 + 321ε + 5223ε² + 1843ε³ + 5697ε⁴ + 4214ε⁵ + 7588ε⁶ + 3976ε⁷ + 2096ε⁸.

On applique alors une approximation de Padé sur chacun des coefficients, pour obtenir P^ε :

X₁⁴

2346ε⁴ + 2660ε³ + 3938ε² + 4726ε

ε⁴ + 976ε³ + 1411ε² + 5186ε + 2487

X₁³

2916ε⁴ + 8872ε³ + 1851ε² + 4762ε

ε x⁴ + 976ε³ + 1411ε² + 5186ε + 2487

X₁²

1207ε⁴ + 8180ε³ + 4499ε² + 4487ε

ε⁴ + 976ε³ + 1411ε² + 5186ε + 2487

X₁

4696ε⁴ + 6406ε³ + 8269ε² + 1053ε + 6514

ε⁴ + 976ε³ + 1411ε² + 5186ε + 2487

En faisant ε=1, on retrouve P=X₁⁴ + 7910X₁³ + 1877X₁² + 5851X₁ + 3439. Pour obtenir W₂, on commence par calculer le développement de Taylor de W₂^ε, en utilisant l'équation (1) ; on trouve

(152ε + 3151ε² + 7932ε³ + 2801ε⁴ + 6656ε⁵ + 8728ε⁶ + 667ε⁷ + 5977ε⁸ )X₁³

+ (8923ε + 1644ε² + 2127ε³ + 3934ε⁴ + 8030ε⁵ + 4864ε⁶ + 170ε⁷ + 5966ε⁸)X₁²

+ (4 + 8421ε + 2873ε² + 1813ε³ + 1603ε⁴ + 998ε⁵ + 1190ε⁶ + 1118ε⁷ + 3333ε⁸)X₁

+ 120ε + 1936ε² + 2120ε³ + 5755ε⁴ + 8628ε⁵ + 1888ε⁶ + 4008ε⁷ + 2380ε⁸.

On applique une approximation de Padé aux coefficients pour retrouver W₂^ε :

3998ε⁴ + 8412ε³ + 1851ε² + 8983ε

ε⁴ + 976ε³ + 1411ε² + 5186ε + 2487

X₁³ +

6475ε⁴ + 6053ε³ + 2711ε² + 4036ε

ε⁴ + 976ε³ + 1411ε² + 5186ε + 2487

X₁² +

8732ε⁴ + 6768ε³ + 2455ε² + 442ε + 947

ε⁴ + 976ε³ + 1411ε² + 5186ε + 2487

X₁+

693ε⁴ + 136ε³ + 548ε² + 1407ε

ε⁴ + 976ε³ + 1411ε² + 5186ε + 2487

En faisant ε=1, on trouve W₂=6951 X₁³ + 7703X₁² + 6268X₁ + 8630 ; par un calcul de pgcd étendu, on en déduit V₂ = 4359X₁³ + 2251X₁² + 2523X₁ + 5498. Donc, les solutions du système F₁,F₂ sont décrites par :

⎧
⎨
⎩

X₁⁴ + 7910X₁³ + 1877X₁² + 5851X₁ + 3439 =0,

X₂ = 4359X₁³ + 2251X₁² + 2523X₁ + 5498.

En factorisant le premier polynôme, on trouve par exemple que (2353,5333) et (2334,4888) sont solutions du système (F₁,F₂).

Ce document a été traduit de L^AT_EX par H^EV^EA