**Entrée** F, G de degrés *m* et *n*.

**Sortie** Q et R tels que $F = QG + R$ avec $\deg R < \deg G$.

1. Calculer $T^m F(1/T)$ et $T^n G(1/T)$ (aucune opération arithmétique n'est nécessaire, il suffit d'inverser l'ordre des coefficients).

2. Calculer le quotient $(T^m F(1/T))/(T^n G(1/T))$ mod $T^{m-n+1}$ par une inversion de série formelle suivie d'un produit.

3. En déduire Q en inversant l'ordre des coefficients.

4. En déduire R, qui est donné par $R = F - QG \bmod X^n$.

### Algorithme 4.1 – Division euclidienne rapide.

Si on pense que $\mathbb{A} = \mathbb{R}$, on voit donc que le développement asymptotique de $F/G$ à l'infini a ses premiers termes donnés par Q, puisque $R/G$ tend vers 0 à l'infini (à cause des contraintes de degré sur R et G). Ceci suggère que l'on va obtenir Q par calcul du développement de Taylor de $F/G$ au voisinage de l'infini.

Concrètement, il suffit de ramener d'abord l'infini en zéro (par le changement de variable $X = 1/T$), pour réduire le calcul à la division de séries formelles. Plus précisément, le point de départ est l'identité

$$ \frac{T^m F(1/T)}{T^n G(1/T)} = T^{m-n} Q(1/T) + \frac{T^m R(1/T)}{T^n G(1/T)}. $$

Dans cette identité les numérateurs et dénominateurs des fractions sont des polynômes, et le polynôme $T^n G(1/T)$ a 1 pour terme constant. En outre, la valuation du second sommant du membre droit est supérieure à $m-n$. En découle donc l'Algorithmme 4.1 de division rapide.

**Théorème 4.1** Soient G un polynôme unitaire de $\mathbb{A}[X]$ de degré $n$ et $F \in \mathbb{A}[X]$ de degré $m \ge n$. Alors on peut calculer le quotient Q et le reste R de la division euclidienne de F par G en $5M(m-n) + M(n) + O(m)$ opérations $(+, -, \times)$ de $\mathbb{A}$.

Ainsi, la division euclidienne ne coûte pas plus cher que la multiplication (à une constante près). En particulier, si on utilise une multiplication à base de FFT, le coût de la division est linéaire en le degré, à des facteurs logarithmiques près.

*Démonstration.* D'après la Proposition 3.4 et la remarque qui la suit, l'inverse du dénominateur à l'étape (2) s'obtient en

$$ 4M(m-n) + O(m-n) $$

opérations dans $\mathbb{A}$, puis une multiplication supplémentaire donne Q. Pour retrouver R, le dernier produit est effectué modulo $X^n$, c'est-à-dire pour un coût de $M(n)$. Toutes les additions sont prises en compte dans le terme $O(m)$. ■
R

Si Q a un degré beaucoup plus petit que G, on peut accélérer le calcul de QG, en découpant G en « tranches » de taille $m-n$; de la sorte, le dernier produit peut se faire en $\frac{n}{m-n}M(m-n)$ opérations. Enfin, si de nombreuses réductions sont à faire modulo le même polynôme G, il est évidemment recommandé de stocker l'inverse du réciproque $T^nG(1/T)$ de G.

## Le cas des entiers

Comme pour la multiplication, il est possible de poser le problème dans $\mathbb{Z}$ comme dans un anneau de polynômes. C'est ce dernier cas qui est le plus simple à étudier, puisqu'il n'y a pas à y gérer de retenue. On obtient cependant un résultat analogue sur les entiers.

**Théorème 4.2** Soit $M_{\mathbb{Z}}$ une fonction de multiplication pour $\mathbb{Z}$, et $f$ et $g$ deux entiers positifs avec $g \le f \le 2^n$. Soient $q$ et $r$ les quotient et reste de la division euclidienne de $f$ par $g$ :

$$f = qg + r \quad \text{avec} \quad 0 \le r < g.$$

On peut calculer $q$ et $r$ en $O(M_{\mathbb{Z}}(n))$ opérations binaires.

## Application aux calculs modulaires

Une application importante de la division euclidienne rapide est le calcul efficace dans des structures algébriques de type « quotient », par exemple dans n'importe quel corps fini.

**Corollaire 4.3** Si $A$ est un anneau et si $P$ est un polynôme unitaire de degré $n$ de $A[X]$, alors dans $A[X]/(P)$, l'addition et la multiplication par un élément de $A$ ont une complexité linéaire en $n$, la multiplication a une complexité en $O(M(n))$.

Le calcul de l'inverse (des éléments inversibles) dans $A[X]/(P)$ est un peu plus complexe. Un algorithme de complexité $O(M(n) \log n)$ est présenté au Chapitre 6.

*Démonstration.* Il s'agit de calculer modulo P. L'addition et la multiplication par un scalaire s'effectuent terme à terme, ainsi leur complexité est linéaire en le degré $n$ de P. Ensuite, pour multiplier deux éléments $A + (P)$ et $B + (P)$ de $A[X]/(P)$, on commence par multiplier $A$ et $B$ dans $A[X]$, puis on réduit le résultat modulo P. La complexité découle donc de ce qui précède. ■

La question du calcul dans $\mathbb{Z}/n\mathbb{Z}$ est légèrement plus délicate que son analogue dans $A[X]$, à cause des retenues. Malgré tout, les résultats s'étendent.

## 4.3 Suites récurrentes linéaires à coefficients constants

**Définition 4.1** Une suite $(a_n)_{n \ge 0}$ d'éléments de l'anneau $A$ est appelée suite récurrente linéaire à coefficients constants (srlcc) d'ordre $d$ si elle satisfait une récurrence
de la forme

$$a_{n+d} = p_{d-1}a_{n+d-1} + \cdots + p_0a_n, \quad n \ge 0,$$

où les $p_i$ sont des éléments de $\mathbb{A}$. Le polynôme $P = X^d - p_{d-1}X^{d-1} - \cdots - p_0$ de $\mathbb{A}[X]$ est appelé *polynôme caractéristique* de la suite $(a_n)_{n \ge 0}$.

Il faut noter que tout multiple d'un polynôme caractéristique est aussi caractéristique. Le pgcd de ces polynômes caractéristiques est appelé *polynôme minimal* de la suite.

### Exemples

**Exercice 4.1** Soit $M$ une matrice de taille $d \times d$ à coefficients dans $\mathbb{A}$, de polynôme caractéristique $\chi_M(X) = \det(XI_d - M)$, soient $u$ une matrice ligne et $v$ une matrice colonne. Montrer que la suite $(uM^n v)_{n \ge 0}$ est une suite récurrente linéaire à coefficients constants admettant $\chi_M$ comme polynôme caractéristique. ■

**Exercice 4.2** Soit $a_1, \dots, a_r$ des entiers strictement positifs. Soit $A_n$ et $B_n$ le nombre de $r$-uplets non ordonnés, resp. ordonnés, $(x_1, \dots, x_r) \in \mathbb{N}^r$, solutions de l'équation $a_1x_1 + \dots + a_rx_r = n$. Montrer que les suites $(A_n)$ et $(B_n)$ sont récurrentes linéaires à coefficients constants, admettant $(X^{a_1} - 1)\cdots(X^{a_r} - 1)$ et $X^{a_1} + \dots + X^{a_r} - 1$ comme polynômes caractéristiques. ■

**Exercice 4.3** Si $P \in \mathbb{A}[X]$ est de degré $d$, alors la suite $(P(n))_{n \ge 0}$ est une srlcc de polynôme caractéristique $(X - 1)^{d+1}$. ■

### Méthode naïve

L'approche naïve pour le calcul des N premiers termes d'une srlcc consiste tout simplement à dérouler la récurrence et requiert $O(dN)$ opérations dans $\mathbb{A}$.

### Exponentiation binaire

Pour calculer uniquement le N-ième terme de la suite, une alternative à cette méthode naïve est basée sur l'exponentiation binaire de la matrice compagnon associée à la suite. Par exemple, la récurrence (4.1) se récrit matriciellement

$$ \underbrace{\begin{pmatrix} F_N \\ F_{N+1} \end{pmatrix}}_{C} = \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} F_{N-1} \\ F_N \end{pmatrix} = C^N \begin{pmatrix} F_0 \\ F_1 \end{pmatrix}, \quad N \ge 1. $$

La puissance de la matrice constante $C$ se calcule alors récursivement par

$$ C^N = \begin{cases} (C^{N/2})^2, & \text{si N est pair,} \\ C \cdot (C^{\frac{N-1}{2}})^2, & \text{sinon.} \end{cases} $$

On en déduit que le calcul de $F_N$ peut être effectué sans calculer les précédents en $O(\log N)$ produits de matrices $2 \times 2$, soit $O(\log N)$ opérations dans $\mathbb{A}$.

Cette idée s'étend facilement à des récurrences d'ordre arbitraire.
**Exercice 4.4** Montrer que le N-ième terme de toute récurrence linéaire d'ordre $d$ à coefficients constants peut se calculer en $O(d^3 \log N)$ opérations dans $\mathbb{A}$.

L'utilisation d'une multiplication matricielle rapide abaisse cette complexité vis-à-vis de $d$ à $O(\text{MM}(d) \log N) = O(d^\theta \log N)$ opérations dans $\mathbb{A}$ (voir le Chapitre 8). La dépendance en $d$ reste cependant plus que quadratique. La section suivante présente une meilleure méthode, de complexité quasi-linéaire en $d$.

## Exponentiation modulaire

La complexité arithmétique de l'algorithme de la section précédente (Exercice 4.4) est quasi-optimale par rapport à l'indice N, mais n'est pas très bonne par rapport à l'ordre $d$ de la récurrence.

Une méthode plus efficace est obtenue en exploitant mieux la structure du problème. En effet, les matrices multipliées ne sont pas quelconques, mais sont des puissances d'une matrice compagnon. Cette structure peut être exploitée grâce à l'observation suivante.

**Lemme 4.4** Soit $(a_n)_{n \ge 0}$ une suite récurrente linéaire à coefficients constants de polynôme caractéristique P(X). Soit $\mathbb{B} = \mathbb{A}[X]/(P)$ et $x$ la classe de X dans $\mathbb{B}$. Alors la matrice compagnon C de P est à la fois

1. telle que $(a_{n+1}, \dots, a_{n+d}) = (a_n, \dots, a_{n+d-1}) \cdot C$ pour tout $n \ge 0$;

2. matrice de l'application $\mathbb{A}$-linéaire de multiplication par $x$ dans $\mathbb{B}$ dans la base canonique $\{1, x, \dots, x^{d-1}\}$.

*Démonstration.* Il suffit d'écrire la matrice pour le voir. Si $P = X^d - p_{d-1}X^{d-1} - \dots - p_0$, la matrice C vaut par définition

$$C = \begin{pmatrix} 0 & 0 & \dots & 0 & p_0 \\ 1 & 0 & \dots & 0 & p_1 \\ 0 & 1 & \dots & 0 & p_2 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & p_{d-1} \end{pmatrix}.$$

La première propriété traduit la définition de la récurrence :

$$a_{n+d} = p_{d-1}a_{n+d-1} + \dots + p_0a_n.$$

La seconde exprime que modulo $P$, $X^d = p_{d-1}X^{d-1} + \dots + p_0$.

**Corollaire 4.5** Avec les mêmes notations, pour tous entiers positifs $i$ et $k$,

$$a_{k+i} = (a_i, \dots, a_{d+i-1}) \cdot V_k,$$

où $V_k$ a pour coordonnées les coefficients du polynôme $X^k \mod P$.

*Démonstration.* Le vecteur $V_k$ est la première colonne de $C^k$, et d'après le lemme précédent, cette matrice est celle de la multiplication par $X^k$ dans $\mathbb{B}$. Les coefficients de la première colonne sont donc les coordonnées de $x^k \cdot 1$ dans la base $\{1, x, \dots, x^{d-1}\}$. ■
**Entrée** $(p_0, \dots, p_{d-1})$ les coefficients d'une récurrence linéaire;
$(a_0, \dots, a_{d-1})$ des conditions initiales; un entier N.

**Sortie** Le terme $a_N$ de la suite définie par

$$a_{n+d} = p_{d-1}a_{n+d-1} + \dots + p_0a_n.$$

1. Poser $P = X^d - p_{d-1}X^{d-1} - \dots - p_0$ et $x = X \mod P$.

2. Calculer récursivement $Q = x^N = q_0 + q_1x + \dots + q_{d-1}x^{d-1}$ par

$$x^k = \begin{cases} (x^{k/2})^2, & \text{si } k \text{ est pair,} \\ x \cdot (x^{\frac{k-1}{2}})^2, & \text{sinon.} \end{cases}$$

3. Renvoyer $q_0a_0 + \dots + q_{d-1}a_{d-1}$.

**Algorithme 4.2** – Calcul du N-ième terme d'une srlcc.

Le résultat de complexité est alors le suivant.

**Théorème 4.6** Soit $(a_n)$ une suite récurrente linéaire à coefficients constants, donnée par une récurrence d'ordre $d$ et des conditions initiales $a_0, \dots, a_{d-1}$. Soit $N \ge d$. Le calcul du N-ième terme $a_N$ peut se faire en seulement $O(M(d) \log N)$ opérations dans $\mathbb{A}$.

*Démonstration.* La méthode est présentée dans l'Algorithme 4.2. Chaque étape du calcul récursif a lieu dans $\mathbb{A}[X]/(P)$ et coûte $O(M(d))$ opérations dans $\mathbb{A}$ d'après le Corollaire 4.3, le nombre d'étapes est en $\log N$; la dernière multiplication ne prend que $O(d)$ opérations. La complexité est donc $O(M(d) \log N)$. ■

## 4.4 Développement de fractions rationnelles

Il est également possible de calculer efficacement les N premiers éléments d'une suite récurrente linéaire à coefficients constants, ou de manière équivalente, le développement de Taylor des fractions rationnelles.

### Séries rationnelles et suites récurrentes

Ces deux notions sont deux aspects de la même question.

**Lemme 4.7** Soit $A(X) = \sum_{n \ge 0} a_n X^n$ la série génératrice de la suite $(a_n)_{n \ge 0}$. Les assertions suivantes sont équivalentes :

i. La suite $(a_n)$ est une srlcc, ayant P pour polynôme caractéristique, de degré $d$.

ii. $A(X) = N_0/\bar{P}$ où $\bar{P} = P(1/X)X^d$, pour un certain $N_0 \in \mathbb{K}[X]$ avec $\deg(N_0) < d$. De plus, si P est le polynôme minimal de $(a_n)$, alors $d = \max\{1 + \deg(N_0), \deg(\bar{P})\}$ et $\text{pgcd}(N_0, \bar{P}) = 1$.

*Démonstration.* Supposons que P s'écrive $P = g_d X^d + \dots + g_0$. Pour l'équivalence, on utilise uniquement le fait que le coefficient de $X^{d+i}$ dans $\bar{P} \cdot A(X)$ est égal à $g_d a_{i+d} +$
$\cdots + g_0 a_i$ et que $\bar{P}(0) = g_d \neq 0$.

Soit maintenant P le polynôme minimal de $(a_n)$. On a $\deg(\bar{P}) \leq d$ avec égalité si et seulement si $g_0 \neq 0$, c'est-à-dire $X \nmid P$. Donc $d \geq \max\{1 + \deg(N_0), \deg(\bar{P})\}$. Supposons par l'absurde que cette inégalité est stricte. Alors $X \mid P$ et on a que $P/X$ est aussi polynôme caractéristique de $(a_n)$, ce qui contredit la minimalité de P. Donc $d = \max\{1 + \deg(N_0), \deg(\bar{P})\}$.

Soit enfin $u := \text{pgcd}(N_0, \bar{P})$. Alors $P_u := P/\bar{u}$ est un polynôme de degré $d - \deg(u)$ qui est caractéristique de $(a_n)$, car $\bar{P}/u = \bar{P}_u$ et $(\bar{P}/u) \cdot A(X) = N_0/u$ est un polynôme de degré $< d - \deg(u)$. Par la minimalité, cela implique que $\deg(u) = 0$, donc $N_0$ et $\bar{P}$ sont bien premiers entre eux. ■

Cette équivalence et le Théorème 4.6 ont pour conséquence immédiate un résultat de complexité.

**Corollaire 4.8** Soit $F(X)/G(X)$ dans $A(X)$ de degré au plus $d$, avec $G(0)$ inversible. Le N-ième coefficient du développement en série

$$F(X)/G(X) = \sum_{n \ge 0} a_n X^n$$

peut être calculé en $O(M(d) \log N)$ opérations.

## Développement

On peut calculer le développement de Taylor à l'ordre N d'une fraction rationnelle de degré $d \leq N$ en $O(M(N))$ opérations, en utilisant la méthode de Newton, décrite au Chapitre 3. Mais $O(M(N))$ est en général beaucoup plus gros que les $O(dN)$ opérations requises par la méthode naïve. Une complexité linéaire en N tout en ne croissant pas trop par rapport à $d$ est en fait possible.

**Théorème 4.9** Soit $F(X)/G(X)$ dans $A(X)$ avec numérateur et dénominateur de degrés au plus $d$, et $G(0)$ inversible. Le développement de Taylor de $F(X)/G(X)$ à précision $N \geq d$ peut se calculer en $O(NM(d)/d)$ opérations d'anneau $(+, -, \times)$ dans $A$.

Si l'anneau $A$ permet la FFT, cette estimation de complexité devient $O(N \log d)$ ce qui est quasi-optimal, simultanément vis-à-vis de N et de $d$.

Le point de départ est le Corollaire 4.5 qui montre comment calculer un coefficient de la suite en fonction de coefficients distants. L'idée est de calculer les coefficients par tranches de longueur $d$. Il ne suffit pas d'utiliser $d$ fois le corollaire, mais il faut encore une fois exploiter une structure supplémentaire, qui provient ici de ce que les coefficients cherchés sont consécutifs. Cette idée est résumée par le résultat suivant.

**Lemme 4.10** Soit $F(X)/G(X)$ dans $A(X)$ avec $G(0)$ inversible. Soit $d$ le degré de $G$ et soit $k \geq 0$. Alors les coefficients $a_k, a_{k+1}, \dots, a_{k+d-1}$ du développement de
**Entrée** Deux polynômes $F(X), G(X)$ avec $\deg G = d, G(0) \neq 0$ et $\deg F < d$, un entier $N > d$.

**Sortie** Les coefficients $a_0, \dots, a_N$ du développement $F(X)/G(X) = \sum_{n \ge 0} a_n X^n$.

1. Calculer les $2d-1$ premiers coefficients $a_0, \dots, a_{2d-2}$.

2. Poser $y_1 = X^{d-1} \mod \tilde{G}(X)$ et $y^* = X^d \mod \tilde{G}(X)$, où $\tilde{G}(X) = X^d G(1/X)$.

3. Pour $i = 2, \dots, \lceil N/d \rceil$ :

a. Calculer $y_i = y^* y_{i-1} \mod \tilde{G}(X)$.

b. Calculer $P = (a_{2d-2} + \dots + a_0 X^{2d-2}) y_i$.

c. Extraire $a_{id-1}, a_{id}, \dots, a_{(i+1)d-2}$ les coefficients des polynômes $X^{2d-2}, \dots, X^d, X^{d-1}$ dans P.

4. Renvoyer $a_0, \dots, a_N$.

**Algorithme 4.3** – Développement en série d'une fraction rationnelle.

$F(X)/G(X) = \sum_i a_i X^i$ sont les coefficients de $X^{2d-2}, \dots, X^d, X^{d-1}$ du produit

$$ (a_{2d-2} + \dots + a_0 X^{2d-2})(X^k \mod G(1/X)X^d). $$

**Démonstration.** Le polynôme $G(1/X)X^d$ est polynôme caractéristique de la suite des coefficients d'après le Lemme 4.7. Le coefficient de $X^{k+i}$ pour $i = 0, \dots, d-1$ dans le produit ci-dessus est égal au produit scalaire $(A_i | V_k)$ où $A_i = (a_i, \dots, a_{d+i-1})$, d'après le Corollaire 4.5. ■

**Preuve du Théorème 4.9.** La méthode est résumée dans l'Algorithme 4.3. La première étape se calcule par itération de Newton en $O(M(d))$ opérations dans $\mathbb{A}$, ce qui représente un coût négligeable sur l'ensemble. L'obtention de $y^*$ est une réciture à partir de $\tilde{G}$ et ne demande pas d'opération. Ensuite, chaque itération utilise un produit dans $\mathbb{A}[X]/(\tilde{G})$ et un produit d'un polynôme de degré $2d-2$ par un polynôme de degré $d-1$. Au total, sont donc effectuées $O(NM(d)/d)$ opérations dans $\mathbb{A}$. ■

La traduction de ce résultat au niveau des srlcc est immédiate.

**Corollaire 4.11** Soit $(a_n)$ une suite récurrente linéaire à coefficients constants, donnée par une récurrence d'ordre $d$, et des conditions initiales $a_0, \dots, a_{d-1}$. Soit $N \ge d$. Alors, les termes $a_0, \dots, a_N$ peuvent être calculés en $O(NM(d)/d)$ opérations dans $\mathbb{A}$.

## 4.5 Applications

### Évaluation d'un polynôme sur une progression arithmétique

**Corollaire 4.12** Un polynôme $P$ de degré $d$ peut être évalué aux $N+1 \gg d$ points $b, a+b, 2a+b, \dots, Na+b$ au prix total de $O(NM(d)/d)$ opérations arithmétiques.
*Démonstration.* L'opérateur $\Delta_a(P) = P(X+a) - P(X)$ fait décroître le degré des polynômes. Il s'ensuit que $\Delta_a^{d+1}P = 0$ et donc $P(an + b)$ est une suite récurrente linéaire d'ordre $d+1$ et de polynôme caractéristique $(X-1)^{d+1}$. ■

## Propriétés de clôture

La classe des srlcc admet de nombreuses propriétés de clôture : si $\mathbf{a} = (a_n)_n$ et $\mathbf{b} = (b_n)_n$ sont deux srlcc de polynômes caractéristiques P et Q, alors

1. la somme $\mathbf{a} + \mathbf{b} = (a_n + b_n)_n$ et le produit de Cauchy $\mathbf{a} \star_C \mathbf{b}$, dont le terme général est $\sum_{i=0}^n a_i b_{n-i}$, sont deux srlcc de polynôme caractéristique PQ;

2. le produit d'Hadamard $\mathbf{a} \odot \mathbf{b} = (a_n b_n)_n$ est une srlcc de polynôme caractéristique égal au produit composé $P \otimes Q$ défini au Chapitre 3;

3. la suite $(\sum_{i=0}^n \binom{n}{i} a_i b_{n-i})_n$ est une srlcc de polynôme caractéristique égal à la somme composée $P \oplus Q$ définie au Chapitre 3.

Ils se calculent donc tous en bonne complexité (quasi-linéaire en la taille de la sortie).

**Exercice 4.5** Prouver les assertions précédentes. ■

## Tests de primalité

Une application du calcul rapide d'un terme d'une récurrence est une famille de tests probabilistes de primalité, de complexité polynomiale. L'idée est de construire une suite récurrente $(a_n)$ d'entiers telle que la primalité de $n$ soit équivalente (ou presque équivalente) à $a_n \equiv 0 \mod n$.

Un cas particulier important en est le test de Fermat (pour lequel $a_n = a^{n-1} - 1$) implanté dans la plupart des systèmes de calcul formel. Bien qu'il soit probabiliste (si $n$ ne passe pas le test, $n$ est composé, mais si $n$ passe le test, alors il est premier seulement avec une grande probabilité), sa grande simplicité le rend souvent préférable à d'autres algorithmes sophistiqués.

**Exercice 4.6** Soit $(a_n)$ une srlcc d'ordre $d$. Montrer qu'il existe des constantes entières $c_0, c_1, \dots, c_d$ telles que $p$ divise $c_0 + c_1 a_{p-1} + \dots + c_d a_{p-d}$ dès lors que $p$ est un nombre premier. De plus, pour tout premier $p$, les constantes $c_i$ mod $p$ peuvent être trouvées en $O(M(d))$ opérations arithmétiques dans $\mathbb{A} = \mathbb{Z}/p\mathbb{Z}$. [Indication : pour $\mathbb{A}$ une matrice carrée d'entiers et $p$ un nombre premier, la trace de $A^p$ et celle de $\mathbb{A}$ sont congrues modulo $p$.] ■

Par exemple, si $(F_n)$ est la suite de Fibonacci, alors $p$ divise $F_{p-2} + 3F_{p-1} - 1$ dès lors que $p$ est premier et la réciproque est vraie avec une bonne probabilité. L'exercice précédent fournit un test de primalité similaire au test de Fermat. D'après le Théorème 4.6, son coût est de $O(M(d)\log p)$ opérations arithmétiques dans $\mathbb{Z}/p\mathbb{Z}$, soit $O(M(d)M_{\mathbb{Z}}(\log p)\log p)$ opérations binaires.

**Exercice 4.7** Soient $a$ et $N$ deux entiers premiers entre eux. Montrer que $N$ est premier si et seulement si $X^N + a = (X + a)^N \mod N$ dans $\mathbb{Z}[X]$. ■

**Exercice 4.8** Montrer que si $N$ est premier, $0 \le a < N$ et $P(X) \in \mathbb{Z}[X]$, alors $X^N + a =$
$(X + a)^N \bmod P(X)$ dans $\mathbb{Z}/N\mathbb{Z}[X]$; si de plus, $P(X)$ est de degré $r = O(\log^c N)$, pour un $c > 0$, alors cette égalité peut être testée « en temps polynomial », c'est-à-dire en un nombre d'opérations binaires polynomial en $\log N$. ■

## Notes

La suite de Fibonacci, introduite en 1202 par Leonardo Pisano (mieux connu sous le pseudonyme de Fibonacci) dans un problème récréatif décrivant la croissance d'une population de lapins jouit de riches propriétés algébriques, arithmétiques et combinatoires. Par exemple : (a) $F_n$ est le nombre de façons différentes de paver un rectangle $2 \times (n-1)$ au moyen de dominos $2 \times 1$; (b) $(F_n)_n$ est une suite de divisibilité, c'est-à-dire que $F_n$ divise $F_m$ dès lors que $n$ divise $m$; (c) si $n$ est impair, alors

$$F_n = 2^{n-1} \prod_{k=1}^{\frac{n-1}{2}} \left( \frac{1}{4} + \cos^2 \frac{k\pi}{n} \right).$$

**Exercice 4.9** Prouver les assertions (a)–(c). ■

Malgré sa simplicité, la suite de Fibonacci fait l'objet de nombreux problèmes ouverts. Par exemple, on ignore s'il existe une infinité de nombres de Fibonacci premiers. Les nombres de Fibonacci sont omniprésents en mathématiques et en informatique<fnref n="2" /> : ils interviennent aussi bien dans l'analyse de l'algorithme d'Euclide pour le calcul du plus grand commun diviseur de deux entiers, que dans la solution négative du dixième problème de Hilbert par Matiyasevich<fnref n="3" />.

Historiquement, le premier algorithme rapide pour la division des polynômes est dû à Moenck et Borodin [MB72]. Son point clé est que le quotient de la division euclidienne de deux polynômes ne dépend que de leurs coefficients de poids fort. Cette remarque est également à la base du calcul rapide de pgcd, étudié au Chapitre 6.

La paternité de l'algorithme de la Section 4.2 revient à Strassen [Str73]. Une alternative, de même complexité asymptotique, à l'algorithme esquissé en Section 4.2 pour les calculs modulaires a été proposée par Montgomery [Mon85].

Calculer les N premiers termes d'une srlcc de polynôme caractéristique fixé est une opération linéaire en les conditions initiales ; c'est le dual de l'opération de division d'un polynôme de degré N par un polynôme fixé. Les conséquences algorithmiques de ce fait seront décrites au Chapitre 12.

Le théorème de Skolem-Mahler affirme que pour toute srlcc $(a_n)$, l'ensemble de ses zéros (les indices $i$ pour lesquels $a_i = 0$) est la réunion d'un ensemble fini et d'un nombre fini de suites arithmétiques. Son étude est une question subtile. Par exemple, déterminer si l'ensemble des zéros est vide est un problème NP-difficile [BP02]. Le livre d'Everest, van der Poorten, Shparlinski et Ward [Eve+03] est une excellente référence pour en savoir plus sur les questions reliées aux srlcc. Sur ces questions, il existe aussi quelques articles synthétiques et bien écrits [CMP87; Poo89].

<footnote n="2">2. Le journal *The Fibonacci Quarterly* est entièrement dédié à l'étude de leurs propriétés.</footnote>

<footnote n="3">3. Ce problème proposait de trouver un algorithme pour décider si un système d'équations diophan- tiennes (polynômes à coefficients entiers) admet une solution en nombres entiers.</footnote>
L'Exercice 4.1 est le point clé de la méthode de Wiedemann pour la résolution de systèmes linéaires creux, traitée au Chapitre 9.

Le calcul rapide d'un terme de la suite de Fibonacci par exponentiation binaire de la matrice compagnon associée relève du folklore mathématique. Sa généralisation (Exercice 4.4) est décrite par Miller et Brown [MB66], mais était probablement connue bien avant.

Le Théorème 4.9 et le Corollaire 4.11 sont dus à Fiduccia [Fid85] et Shoup [Sho91a]. Leur récente généralisation au cas des matrices polynomiales est à la base du meilleur algorithme pour la résolution de systèmes linéaires à coefficients polynomiaux, exposé au Chapitre 11.

Il n'existe pas de test déterministe de primalité basés sur le test modulaire d'un terme d'une récurrence à coefficients constants. Par contre, on peut caractériser un nombre premier N à l'aide de suites qui vérifient des récurrences à coefficients polynomiaux (comme la factorielle, via le test de Wilson $(N-1)! = -1 \pmod N$). Malheureusement, cela ne fournit pas d'algorithmes efficace. Le premier algorithme déterministe qui prouve la primalité en temps polynomial est très récent [AKS04]; il part de la caractérisation de type Fermat donnée dans l'Exercice 4.8, et exhibe une constante c et un polynôme P tels que la primalité de N est impliquée par la vérification de l'identité de l'Exercice 4.8 pour seulement $r^{1/2} \log N$ valeurs de a.

Bibliographie

AKS04 AGRAWAL, Manindra, Neeraj KAYAL et Nitin SAXENA (2004). « PRIMES is in P ». In: *Annals of Mathematics. Second Series*, vol. 160, n°2, p. 781-793.

BP02 BLONDEL, Vincent D. et Natacha PORTIER (2002). « The presence of a zero in an integer linear recurrent sequence is NP-hard to decide ». In: *Linear Algebra and its Applications*, vol. 351/352. Fourth special issue on linear systems and control, p. 91-98.

CMP87 CERLIENCO, L., M. MIGNOTTE et F. PIRAS (1987). « Suites récurrentes linéaires. Propriétés algébriques et arithmétiques ». In: *L'Enseignement Mathématique*. II, vol. 33, p. 67-108.

Eve+03 EVEREST, Graham, Alf van der POORTEN, Igor SHPARLINSKI et Thomas WARD (2003). *Recurrence sequences*. Vol. 104. Mathematical Surveys and Monographs. American Mathematical Society.

Fid85 FIDUCCIA, C. M. (1985). « An efficient formula for linear recurrences ». In: *SIAM Journal on Computing*, vol. 14, n°1, p. 106-112.

MB66 MILLER, J. C. P. et D. J. Spencer BROWN (1966). « An algorithm for evaluation of remote terms in a linear recurrence sequence ». In: *Computer Journal*, vol. 9, p. 188-190.

MB72 MOENCK, R. T. et A. BORODIN (1972). « Fast modular transforms via division ». In: *Thirteenth Annual IEEE Symposium on Switching and Automata Theory*, p. 90-96.

Mon85 MONTGOMERY, Peter L. (1985). « Modular multiplication without trial division ». In: *Mathematics of Computation*, vol. 44, n°170, p. 519-521.
Poo89 POORTEN, A. J. van der (1989). « Some facts that should be better known, especially about rational functions ». In : *Number theory and applications*. Proceedings of a Conference held at Banff, AB, 1988. Dordrecht : Kluwer, p. 497–528.

Sho91a SHOUP, V. (1991). « A fast deterministic algorithm for factoring polynomials over finite fields of small characteristic ». In : *ISSAC'91 : International Symposium on Symbolic and Algebraic Computation*. ACM Press, p. 14–21.

Str73 STRASSEN, V. (1972/73). « Die Berechnungskomplexität von elementar-symmetrischen Funktionen und von Interpolationskoeffizienten ». In : *Numerische Mathematik*, vol. 20, p. 238–251.
<figure><img src="image_15.png" /></figure>
# 5. Calculs modulaires, évaluation et interpolation

## Résumé

Le concept d'algorithme modulaire est central en calcul formel. Il permet de pallier le phénomène d'explosion des expressions intermédiaires. Des cas particuliers importants sont l'évaluation multipoint et l'interpolation, pour lesquels il existe des algorithmes rapides qui peuvent être vus comme des généralisations de la FFT.

## 5.1 Introduction

Un ingrédient essentiel en calcul formel est l'utilisation de différents types de représentation pour les objets manipulés. Par exemple, un polynôme est classiquement codé par la liste de ses coefficients, mais on peut très bien le représenter par les valeurs qu'il prend en un nombre suffisant de points. D'ailleurs, nous avons vu au Chapitre 2 que c'est bien cette seconde représentation qui est la plus adaptée pour le calcul efficace, via la FFT, du produit de polynômes. Par ailleurs, d'autres opérations (comme la division polynomiale, traitée au Chapitre 4) se font mieux dans la première représentation. D'autres exemples de codages alternatifs déjà rencontrés dans cet ouvrage sont l'écriture des entiers en différentes bases de numération, ou encore la représentation des suites récurrentes linéaires à coefficients constants, soit par leurs éléments, soit par leurs séries génératrices, soit par une récurrence et des conditions initiales.

L'exemple de la FFT montre à quel point il est crucial d'examiner dans quelle représentation un problème donné est plus facile à traiter, et aussi, de trouver des algorithmes rapides pour la conversion d'une représentation à l'autre. Une incarna- tion de ce concept général est la notion d'algorithme modulaire, dont la méthodologie
« évaluation-interpolation » est un cas particulier très important. L'approche modulaire consiste à choisir des modules $m_i$, à faire des calculs modulo chaque $m_i$ et à reconstruire tout à la fin le résultat à l'aide d'une version effective du théorème des restes chinois. Pour garantir l'unicité du résultat et la correction de ce schéma, il suffit que les modules soient suffisamment nombreux et indépendants. Typiquement, s'il s'agit d'entiers, ils peuvent être choisis premiers entre eux et tels que leur produit dépasse le résultat final. En particulier, pour mettre en place une approche modulaire, on a besoin de bornes a priori sur la taille de l'objet à calculer.

Cette approche porte ses fruits surtout lorsque les coefficients dans les calculs intermédiaires sont beaucoup plus gros que ceux du résultat final. Il s'agit du phénomène d'explosion des expressions intermédiaires, qui se présente par exemple dans le calcul du déterminant d'une matrice entière ou polynomiale, ou encore au cours du calcul du pgcd de polynômes à coefficients entiers.

**Exemple 5.1** — Calcul du déterminant d'une matrice polynomiale. Soit à calculer le déterminant d'une matrice $n \times n$, dont les éléments sont des polynômes de degré au plus $d$ (le cas particulier $d = 1$ correspond au calcul d'un polynôme caractéristique). Le point de départ est la remarque que le pivot de Gauss produit sur la $i$-ième ligne des fractions rationnelles de degré $2^i d$. Ainsi, sa complexité ne semble pas polynomiale par rapport à $n$. Une approche évaluation-interpolation résout cette anomalie. Le déterminant étant un polynôme de degré au plus $nd$, il suffit de choisir un ensemble de $nd + 1$ points, d'y évaluer les éléments de la matrice, de calculer les $nd + 1$ déterminants de matrices scalaires et de finir par une interpolation fournissant le determinant cherché. Le même raisonnement s'applique aux matrices entières.

*Choix des modules.* Dans certaines situations, le programmeur a le choix des modules. Dans l'exemple précédent, on peut choisir comme modules des polynômes de la forme $X - \omega^i$, dans le cas favorable où l'anneau de base contient une racine principale de l'unité $\omega$ d'ordre suffisamment élevé ($\approx 2dn$). En théorie, ce choix est donc possible dès lors que l'anneau de base permet une multiplication polynomiale à base de FFT. En pratique, il existe des plages de degrés pour lesquelles, même si elle est faisable, la FFT n'est pas encore rentable; ce choix est alors déconseillé. Dans d'autres situations, le choix des modules est intrinsèquement imposé par le problème à traiter; c'est le cas pour le calcul de la factorielle exposé au Chapitre 15.

## 5.2 Présentation, résultats

Les problèmes d'évaluation et d'interpolation sont inverses l'un de l'autre. Dans leur version standard, ils s'énoncent ainsi. Soient $a_0, \dots, a_{n-1}$ des points dans $\mathbb{A}$, où $\mathbb{A}$ est un anneau commutatif.

**Évaluation multipoint** Étant donné un polynôme $P$ dans $\mathbb{A}[X]$, de degré strictement inférieur à $n$, calculer les valeurs :

$$P(a_0), \dots, P(a_{n-1}).$$

**Interpolation** Étant donnés $b_0, \dots, b_{n-1} \in \mathbb{A}$, trouver un polynôme $P \in \mathbb{A}[X]$ de degré
strictement inférieur à $n$ tel que

$$P(a_0) = b_0, \dots, P(a_{n-1}) = b_{n-1}.$$

Le problème d'interpolation admet toujours une solution unique sous l'hypothèse technique suivante :

$$a_i - a_j \text{ est inversible dans } \mathbb{A} \text{ lorsque } i \neq j. \qquad (\text{H})$$

En effet, si $V = (a_{i-1}^{j-1})_{i,j=1}^n$ est la matrice de Vandermonde associée aux points $a_i$, dont le déterminant vaut $\det(V) = \prod_{i<j} (a_j - a_i)$, alors le problème d'interpolation se traduit par la résolution en l'inconnue $p$ du système linéaire $Vp = b$, où $b$ est le vecteur des $b_i$. Or, ce système admet une solution unique, puisque l'hypothèse (H) entraîne l'inversibilité du déterminant $\det(V)$ et donc aussi celle de la matrice $V$.

**Exercice 5.1** Montrer que $\det(V) = \prod_{i<j} (a_j - a_i)$.

■

Cette discussion suggère un algorithme naïf pour l'interpolation, produisant l'unique solution $\mathbf{p} = \mathbf{V}^{-1}\mathbf{b}$ du système $V\mathbf{p} = \mathbf{b}$ à l'aide du pivot de Gauss en $O(n^3)$ opérations dans $\mathbb{A}$, ou encore en $O(\text{MM}(n)) = O(n^\theta)$ opérations ($2 \le \theta < 3$), en utilisant l'algorithmique matricielle rapide (Chapitre 8). De manière analogue, l'évaluation multipoint de $P$ en les $a_i$ se traduit par le produit matrice-vecteur $V\mathbf{p}$, où $\mathbf{p}$ est le vecteur des coefficients de $P$; elle peut donc être effectuée de manière naïve en $O(n^2)$ opérations arithmétiques. En s'y prenant vraiment naïvement, l'interpolation est donc plus coûteuse que l'évaluation multipoint. Nous verrons un peu plus loin qu'une méthode exploitant la formule d'interpolation de Lagrange permet de résoudre le problème d'interpolation en complexité quadratique en $n$.

L'objectif de ce chapitre est de montrer qu'il est possible d'atteindre une complexité quasi-optimale, tant pour l'évaluation multipoint que pour l'interpolation, en utilisant des algorithmes de type « diviser pour régner » qui ramènent ces questions à des multiplications de polynômes. Les principaux résultats sont les suivants :

**Théorème 5.1** — « Évaluation-interpolation » . On peut effectuer l'évaluation et l'interpolation sur $n$ points vérifiant l'hypothèse (H) en utilisant $O(M(n) \log n)$ opérations $(+, -, \times)$ de $\mathbb{A}$. Cette complexité peut être abaissée à $O(M(n))$ si les points forment une progression géométrique de raison inversible dans $\mathbb{A}$.

En termes pratiques, si l'anneau de base $\mathbb{A}$ est un corps fini, les algorithmes rapides deviennent avantageux pour des degrés de l'ordre de quelques dizaines : c'est sensiblement mieux que pour les algorithmes rapides de type « diviser pour régner » utilisés aux Chapitres 6 et 7 pour le calcul rapide de pgcd ou d'approximants de Padé-Hermite; cela reflète une relative simplicité des algorithmes.

L'évaluation multipoint et l'interpolation sont des instances particulières des problèmes *modules multiples* et *restes chinois* s'énonçant comme suit :

**Modules multiples** Étant donnés des polynômes $P, m_1, \dots, m_r$ dans $\mathbb{A}[X]$ qui satisfont $\deg(P) < n = \sum_i \deg(m_i)$, calculer les restes :

$$P \bmod m_1, \dots, P \bmod m_r.$$
**Restes chinois** Étant donnés des polynômes $b_1, \dots, b_r, m_1, \dots, m_r$ dans $\mathbb{A}[X]$, avec $m_i$ unitaires, trouver $P \in \mathbb{A}[X]$ de degré strictement inférieur à $n = \sum_i \deg(m_i)$ tel que

$$P \bmod m_1 = b_1, \dots, P \bmod m_r = b_r.$$

Les techniques de ce chapitre s'y généralisent. Pour les restes chinois, l'unicité du résultat est garantie par l'hypothèse que les résultats $\mathrm{Res}(m_i, m_j)$ sont inversibles dans $\mathbb{A}$, pour tous $i \neq j$. Cette condition technique généralise la condition d'inversibilité des $a_i - a_j$ (se rapporter au Chapitre 6 pour la définition du résultant). On peut alors prouver le résultat suivant, laissé en exercice :

**Théorème 5.2** On peut résoudre les deux problèmes ci-dessus avec un algorithme qui utilise $O(M(n) \log n)$ opérations $(+, -, \times)$ dans $\mathbb{A}$.

Naturellement, on peut poser les questions précédentes dans le cas où l'anneau de polynômes $\mathbb{A}[X]$ est remplacé par l'anneau $\mathbb{Z}$, les polynômes $m_i$ sont remplacés par des modules entiers $a_i$. Il est possible d'obtenir le même type de résultats de complexité, cette fois en comptant les opérations binaires.

**Théorème 5.3** Soient $a_1, \dots, a_r$ des entiers strictement positifs de produit $A = a_1 \cdots a_r$ ayant au plus $n$ chiffres. Si $b$ est un entier d'au plus $n$ chiffres, alors on peut calculer $b \bmod a_1, \dots, b \bmod a_r$ en $O(M(\mathbb{Z}(n) \log n))$ opérations binaires.

Inversement, sous l'hypothèse que les $a_i$ sont premiers entre eux, à partir d'entiers $b_1, \dots, b_r$ tels que $b_i < a_i$, on peut calculer l'unique entier $b < A$ tel que $b \bmod a_i = b_i$ pour tout $i$, en $O(M(\mathbb{Z}(n) \log n))$ opérations binaires.

La suite de ce chapitre est organisée comme suit : dans la Section 5.3 est décrite la méthode d'interpolation de Lagrange, de complexité quadratique. La Section 5.4 est consacrée aux algorithmes rapides sous-jacents à la preuve du Théorème 5.1. Comme les techniques pour le cas des entiers sont les mêmes que pour les polynômes, nous laisserons la preuve du Théorème 5.3 à titre d'exercice.

## 5.3 Interpolation de Lagrange

Un algorithme de complexité quadratique pour l'interpolation polynomiale repose sur une écriture explicite du polynôme interpolant. Soient $b_i$ les valeurs à interpoler et $a_i$ les points d'interpolation vérifiant l'hypothèse (H). On peut alors écrire $P$ sous la forme :

$$P(X) = \sum_{i=0}^{n-1} b_i \prod_{0 \le j \le n-1, j \ne i} \frac{X - a_j}{a_i - a_j}.$$

Cette égalité est usuellement appelée la formule d'interpolation de Lagrange. Pour la prouver, il suffit d'observer que pour tout $i$, le produit s'annule en $a_j$ ($j \neq i$), et vaut 1 en $a_i$. Afin de simplifier l'écriture de la formule de Lagrange, posons

$$A = \prod_j (X - a_j) \quad \text{et} \quad A_i = \prod_{j \neq i} (X - a_j) = \frac{A}{X - a_i},$$
**Entrée** $a_0, \dots, a_{n-1} \in \mathbb{A}$ vérifiant (H) et $b_0, \dots, b_{n-1}$ dans $\mathbb{A}$.

**Sortie** L'unique polynôme $P \in \mathbb{A}[X]$ de degré strictement inférieur à $n$ tel que $P(a_i) = b_i$ pour tout $i$.

1. Initialiser $A$ à 1, et $P$ à 0.

2. Pour $i = 0, \dots, n-1$ faire
   $A = A \cdot (X - a_i)$.

3. Pour $i = 0, \dots, n-1$ faire
   $A_i = A/(X - a_i)$
   $q_i = A_i(a_i)$
   $P = P + b_i A_i / q_i$.

4. Renvoyer $P$.

**Algorithme 5.1** – Algorithme de Lagrange pour l'interpolation polynomiale.

pour $i = 0, \dots, n-1$. On a donc les relations $A(a_i) = 0$ et pour $j \neq i$, $A_j(a_i) = 0$. On obtient alors l'écriture

$$P = \sum_{i=0}^{n-1} b_i \frac{A_i(X)}{A_i(a_i)}.$$

Une approche directe pour l'interpolation revient à calculer tout d'abord les polynômes $A_i$, puis leurs valeurs en les points $a_i$, et enfin à effectuer les combinaisons linéaires nécessaires. Le seul point non trivial est le calcul des $A_i$ : si on n'y fait pas attention, on risque de sortir des bornes en $O(n^2)$ (par exemple, si on les calcule tous indépendamment, et chacun de manière quadratique). Pour faire mieux, on partage le gros des calculs, en calculant d'abord le polynôme $A$.

**Proposition 5.4** L'Algorithme 5.1 utilise $O(n^2)$ opérations dans $\mathbb{A}$.

*Démonstration.* La multiplication d'un polynôme de degré $d$ par un polynôme de degré 1 prend un temps linéaire en $d$, disons $Cd$ opérations, C étant une constante. Calculer $A$ demande donc $C(1 + 2 + \cdots + (n-1)) = O(n^2)$ opérations. Ensuite, chaque passage dans la seconde boucle prend un temps linéaire en $n$. Il y a $n$ passages, d'où à nouveau un coût quadratique. ■

## 5.4 Algorithmes rapides

### Les idées

L'idée fondamentale de l'algorithme d'évaluation multipoint rapide est d'utiliser une technique de type « diviser pour régner ». Supposons pour simplifier que $n$ soit pair et séparons l'ensemble des points $a_0, \dots, a_{n-1}$ en deux paquets, $a_0, \dots, a_{n/2-1}$ et $a_{n/2}, \dots, a_{n-1}$. Il est naturel d'associer à $P$ deux polynômes $P_0$ et $P_1$ de degrés strictement inférieurs à $n/2$, tels que :

$$P_0(a_0) = P(a_0), \dots, P_0(a_{n/2-1}) = P(a_{n/2-1}) \quad \text{et} \quad P_1(a_{n/2}) = P(a_{n/2}), \dots, P_1(a_{n-1}) = P(a_{n-1}).$$
<figure><img src="image_16.png" /><figcaption>FIGURE 5.2 – Arbre des sous-produits.</figcaption></figure>

Entrée $a_0, \dots, a_{n-1}$ dans $\mathbb{A}$, avec $n$ de la forme $n = 2^\kappa$ ($\kappa \in \mathbb{N}$).

Sortie L'arbre $\mathcal{B}$ des sous-produits associé aux points $a_i$.

1. Si $n=1$, renvoyer $X-a_0$.

2. Calculer récursivement l'arbre $\mathcal{B}_0$ des sous-produits associé à $a_0, \dots, a_{n/2-1}$.

3. Calculer récursivement l'arbre $\mathcal{B}_1$ des sous-produits associé à $a_{n/2}, \dots, a_{n-1}$.

4. Calculer $A = B_0 B_1$, où $B_0$ et $B_1$ sont respectivement les racines de $\mathcal{B}_0$ et de $\mathcal{B}_1$.

5. Renvoyer l'arbre de racine $A$ et dont les fils sont $\mathcal{B}_0$ et $\mathcal{B}_1$.

Algorithme 5.3 – Algorithme rapide pour le calcul de l'arbre des sous-produits.

Pour effectuer cette construction, le choix suivant s'impose :

$$P_0 = P \mod (X - a_0) \cdots (X - a_{n/2-1}), \quad P_1 = P \mod (X - a_{n/2}) \cdots (X - a_{n-1}).$$

On a donc ramené le calcul en degré $n$ à deux calculs en degré $n/2$, plus deux divisions euclidiennes; on va ensuite itérer ce processus. En bout de course, on retrouve les valeurs $P(a_i)$, car celles-ci peuvent également s'écrire $P \mod (X - a_i)$.

Remarquons que les polynômes par lesquels on effectue les divisions euclidiennes ne sont pas « gratuits » : on doit les (pré)calculer. Pour cela, l'idée est de les empiler dans une structure d'arbre binaire, dont les nœuds internes sont étiquetés par des polynômes.

## L'arbre des sous-produits

Pour simplifier, dans tout ce qui suit, on suppose que le nombre de points est une puissance de 2, $n = 2^\kappa$. Quand ce n'est pas le cas, on adapte toutes les constructions (on fabrique toujours un arbre binaire, mais il lui manque alors des nœuds internes
<figure><img src="image_17.png" /><figcaption>FIGURE 5.4 – Parcours descendant de l’arbre des sous-produits.</figcaption></figure>

à certains niveaux). On peut donner la définition récursive suivante pour cet arbre (noté $\mathcal{B}$), représenté dans la Figure 5.2 :

- Si $\kappa = 0$, $\mathcal{B}$ se réduit à un nœud unique qui contient le polynôme $X - a_0$.

- Sinon, soient $\mathcal{B}_0$ et $\mathcal{B}_1$ les arbres associés respectivement aux points $a_0, \dots, a_{2^{\kappa-1}-1}$ et $a_{2^{\kappa-1}}, \dots, a_{2^{\kappa}-1}$. Soient $B_0$ et $B_1$ les polynômes présents aux racines de $\mathcal{B}_0$ et $\mathcal{B}_1$. Alors $\mathcal{B}$ est l’arbre dont la racine contient le produit $B_0B_1$ et dont les fils sont $\mathcal{B}_0$ et $\mathcal{B}_1$.

L’intérêt de cette structure d’arbre est double : d’une part, il contient tous les polynômes par lesquels on va diviser lors de l’algorithme d’évaluation, d’autre part, le coût de son calcul s’amortit, et devient essentiellement linéaire en le degré.

**Proposition 5.5** L’Algorithme 5.3 calcule tous les polynômes apparaissant dans l’arbre $\mathcal{B}$ en $O(M(n) \log n)$ opérations arithmétiques dans $\mathbb{A}$.

*Démonstration.* Soit $T(n)$ le coût du calcul de l’arbre pour $n$ points. La définition récursive implique l’inégalité suivante pour $T(n)$ :

$$T(n) \le 2T\left(\frac{n}{2}\right) + M\left(\frac{n}{2}\right).$$

Le résultat s’ensuit en invoquant le théorème « diviser pour régner ».

<figure><img src="image_18.png" /></figure>

## Évaluation multipoint rapide

L’algorithme d’évaluation multipoint devient simple à écrire de manière récursive, si on s’autorise un peu de souplesse dans l’écriture, en supposant précalculés tous les nœuds de l’arbre. Comme illustré dans la Figure 5.4, cela correspond à faire *descendre* les restes de divisions euclidiennes dans l’arbre des sous-produits. Or, la division avec reste par un polynôme unitaire de degré $n$ peut s’effectuer en $O(M(n))$ opérations de $\mathbb{A}$ (Chapitre 4). Nous obtenons le résultat de complexité suivant :
**Entrée** $a_0, \dots, a_{n-1}$ dans $\mathbb{A}$, avec $n$ de la forme $n = 2^\kappa$ ($\kappa \in \mathbb{N}$), l'arbre $\mathcal{B}$ des sous-produits associé aux points $a_i$, P dans $\mathbb{A}[X]$, avec deg P < $n$.

**Sortie** $P(a_0), \dots, P(a_{n-1})$.

1. Si $n=1$, renvoyer P.

2. $P_0 = P \bmod (X - a_0) \cdots (X - a_{n/2-1})$.

3. $P_1 = P \bmod (X - a_{n/2}) \cdots (X - a_{n-1})$.

4. Calculer (par un appel récursif) $L_0 = P_0(a_0), \dots, P_0(a_{n/2-1})$.

5. Calculer (par un appel récursif) $L_1 = P_1(a_{n/2}), \dots, P_1(a_{n-1})$.

6. Renvoyer $L_0, L_1$.

**Algorithme 5.5** – Algorithme rapide d'évaluation multipoint par division répétée.

**Proposition 5.6** La complexité de l'Algorithme 5.5 est de $O(M(n) \log n)$ opérations arithmétiques dans $\mathbb{A}$.

*Démonstration.* Soit $T(n)$ le coût du calcul pour $n$ points. L'algorithme ci-dessus implique la récurrence suivante pour $T(n)$ :

$$T(n) \le 2T\left(\frac{n}{2}\right) + O(M(n)),$$

et le résultat en découle grâce au théorème « diviser pour régner ».

**Exercice 5.2** Si $\omega \in \mathbb{A}$ est une racine principale $n$-ième de l'unité, montrer que l'algorithme FFT vu au Chapitre 2 est un cas particulier de l'Algorithme 5.5.

**Exercice 5.3** Estimer la complexité de la variante de l'Algorithme 5.5 dans laquelle l'arbre des sous-produits n'est pas précalculé, les polynômes modulo lesquels s'effectuent les divisions étant calculés à la volée.

## Sommes de fractions

Un problème intermédiaire à traiter avant de résoudre efficacement le problème d'interpolation est celui du calcul rapide d'une somme de fractions. Rappelons les définitions introduites dans la Section 5.3. À partir des points $a_i$, nous y avions défini les polynômes

$$A = \prod_j (X - a_j) \quad \text{et} \quad A_i = \prod_{j \neq i} (X - a_j) = \frac{A}{X - a_i}.$$

Il était apparu que pour effectuer l'interpolation, il suffisait de savoir effectuer la tâche suivante : étant données des valeurs $c_i$, calculer le numérateur et le dénominateur de la fraction rationnelle

$$\sum_{i=0}^{n-1} \frac{c_i}{X - a_i} \qquad (5.1)$$
**Entrée** $a_0, \dots, a_{n-1} \in \mathbb{A}$ et $c_0, \dots, c_{n-1}$ dans $\mathbb{A}$, avec $n$ de la forme $n = 2^\kappa$ ($\kappa \in \mathbb{N}$).

**Sortie** $N, P \in \mathbb{A}[X]$ tels que $P = (X - a_0) \cdots (X - a_{n-1})$ et $\frac{N}{P} = \sum_{i=0}^{n-1} \frac{c_i}{X - a_i}$.

1. Si $n=1$, alors renvoyer $c_0$ et $X - a_0$.

2. Appeler récursivement l'algorithmme avec $a_0, \dots, a_{n/2-1}$ et $c_0, \dots, c_{n/2-1}$. On note $N_1$ et $P_1$ les polynômes renvoyés.

3. Appeler récursivement l'algorithmme avec $a_{n/2}, \dots, a_{n-1}$ et $c_{n/2}, \dots, c_{n-1}$ pour obtenir en sortie $N_2$ et $P_2$.

4. Renvoyer $N_1 P_2 + N_2 P_1$ et $P_1 P_2$.

## Algorithme 5.6 – Somme de fractions.

C'est à cette question que nous allons répondre maintenant ; à nouveau, on utilise une idée de type « diviser pour régner » : par deux appels récursifs, on calcule

$$S_1 = \sum_{i=0}^{n/2-1} \frac{c_i}{X - a_i} \quad \text{et} \quad S_2 = \sum_{i=n/2}^{n} \frac{c_i}{X - a_i}$$

et on renvoie $S = S_1 + S_2$. Soit $T(n)$ le coût du calcul pour $n$ points. La méthode est résumée dans l'Algorithmme 5.6. Son coût satisfait à la récurrence :

$$T(n) \le 2T\left(\frac{n}{2}\right) + O(M(n)),$$

d'où l'on déduit le résultat suivant, en faisant encore une fois appel au théorème « diviser pour régner ».

**Proposition 5.7** L'Algorithmme 5.6 calcule le dénominateur et le numérateur de la fraction rationnelle (5.1) en $O(M(n) \log n)$ opérations dans $\mathbb{A}$.

## Interpolation

À ce stade, l'interpolation ne pose plus de réelle difficulté. Au vu des formules de la Section 5.3, le polynôme interpolant les valeurs $b_i$ aux points $a_i$ est donné par :

$$P(X) = \sum_{i=0}^{n-1} b_i \frac{A_i(X)}{A_i(a_i)} = A(X) \times \sum_{i=0}^{n-1} \frac{b_i/A_i(a_i)}{X - a_i}.$$

D'après les paragraphes précédents, on sait comment calculer rapidement la fraction rationnelle $\sum c_i/(X - a_i)$; il ne reste plus qu'à calculer les constantes $c_i = b_i/A_i(a_i)$. Le calcul itératif des constantes $c_i$ reviendrait à évaluer $n$ polynômes différents, chacun en un point, et aurait donc une complexité trop élevée, quadratique en $n$. La proposition suivante permet de remplacer ces $n$ évaluations par une évaluation multipoint d'un unique polynôme en $n$ points, ce qui mène à un gain de complexité.
**Entrée** $a_0, \dots, a_{n-1} \in \mathbb{A}$ vérifiant (H) et $b_0, \dots, b_{n-1}$ dans $\mathbb{A}$.

**Sortie** L'unique polynôme $P \in \mathbb{A}[X]$ de degré strictement inférieur à $n$ tel que $P(a_i) = b_i$ pour tout $i$.

1. Calculer l'arbre des sous-produits associé aux points $a_i$ (de racine A).

2. Calculer $d_i = A'(a_i)$ pour tout $0 \le i \le n-1$, avec l'Algorithmme 5.5.

3. $(c_0, \dots, c_{n-1}) = (b_0/d_0, \dots, b_{n-1}/d_{n-1})$.

4. Calculer la somme $R$ de fractions

$$c_0/(X - a_0) + \dots + c_{n-1}/(X - a_{n-1})$$

avec l'Algorithmme 5.6.

5. Renvoyer le numérateur de $R$.

**Algorithme 5.7 –** Algorithme d'interpolation rapide.

**Proposition 5.8** Pour tout $i = 0, \dots, n-1$, on a $A_i(a_i) = A'(a_i)$.

*Démonstration.* On dérive $A$, ce qui fournit l'égalité :

$$A' = \sum_{i=0}^{n-1} \prod_{j \neq i} (X - a_j) = \sum_{i=0}^{n-1} A_i.$$

On a vu précédemment que $A_i(a_j) = 0$ pour $j \neq i$, ce qui donne le résultat.

L'algorithmme d'interpolation et l'estimation de son coût s'en déduisent facilement.

**Proposition 5.9** La complexité de l'Algorithmme 5.7 est $O(M(n) \log n)$ opérations de $\mathbb{A}$.

*Démonstration.* C'est une conséquence des Propositions 5.5, 5.6 et 5.7.

## Évaluation et interpolation sur une suite géométrique

Dans cette section nous montrons que tout polynôme de degré $n$ peut être évalué et interpolé sur des points en progression géométrique $\{1, q, \dots, q^{n-1}\}$ en $O(M(n))$ opérations.

**Proposition 5.10** Soit $q$ un élément inversible de $\mathbb{A}$, et $a_i = q^i$ pour $0 \le i < n$.

Alors :

1. On peut évaluer tout polynôme $P \in \mathbb{A}[X]$ de degré strictement inférieur à $n$ sur les points $a_0, \dots, a_{n-1}$ en $O(M(n))$ opérations dans $\mathbb{A}$.

2. Si les points $(a_i)$ vérifient l'hypothèse (H) on peut effectuer l'interpolation en les points $a_0, \dots, a_{n-1}$ en $O(M(n))$ opérations dans $\mathbb{A}$.
Démonstration. Soit $P = p_0 + p_1 X + \dots + p_{n-1} X^{n-1}$. Pour $i = 0, \dots, 2n-2$, on introduit les nombres triangulaires $t_i = i(i-1)/2$ et la suite $\beta_i = q^{t_i}$; pour $i = 0, \dots, n-1$ on construit $\gamma_i = p_i/\beta_i$. Tous les éléments $q^{t_i}$, $\gamma_i$ et $\beta_i$ peuvent être calculés en $O(n)$ opérations, en utilisant l'égalité $q^{t_{i+1}} = q^i q^{t_i}$. L'observation clé est que $ij = t_{i+j} - t_i - t_j$ et donc $q^{ij} = \beta_i^{-1} \beta_j^{-1} \beta_{i+j}$. L'algorithme d'évaluation repose alors sur la formule

$$P(q^i) = \sum_{j=0}^{n-1} p_j q^{ij} = \beta_i^{-1} \cdot \sum_{j=0}^{n-1} \gamma_j \beta_{i+j},$$

qui montre que les valeurs $P(q^i)$ sont, à des facteurs constants près, données par les coefficients de $X^{n-1}, \dots, X^{2n-2}$ dans le produit de $\sum_{i=0}^{n-1} \gamma_i X^{n-i-1}$ et $\sum_{i=0}^{2n-2} \beta_i X^i$.

Pour l'interpolation, on adapte l'Algorithme 5.7, de façon à économiser un facteur log en tirant parti de la structure des points $a_i = q^i$.

On commence par remarquer que la racine A(X) de l'arbre des sous-produits peut être calculée en $O(M(n))$ opérations dans $\mathbb{A}$, par un algorithme qui ne nécessite pas la construction de tout l'arbre des sous-produits, donc différent de l'Algorithme 5.3. En effet, puisque les points $a_i$ forment une suite géométrique, il est possible de calculer récursivement $A = \prod_{i=0}^{n-1} (X - a_i)$ par un seul appel récursif calculant $B_0 = \prod_{i=0}^{n/2-1} (X - a_i)$, suivi du calcul de $B_1 = \prod_{i=n/2}^{n-1} (X - a_i)$ à partir de $B_0$ par une homothétie de coût $O(n)$, et enfin du calcul du produit $A = B_0 B_1$.

Par ailleurs, observons que dans l'Algorithme 5.7, le calcul de tout l'arbre des sous-produits a été utilisé uniquement pour l'évaluation multipoint de $A'$. Or, par la première partie de la Proposition 5.10, l'évaluation de $A'$ sur les points $a_i$ peut s'effectuer en $O(M(n))$ opérations.

Il nous reste à prouver que la somme des fractions (5.1), qui devient ici

$$\frac{N}{A} = \sum_{i=0}^{n-1} \frac{c_i}{X - q^i} \qquad (5.2)$$

peut également être déterminée pour le même coût.

On adopte la stratégie suivante : on calcule le développement en série à l'ordre $n$ de la fraction (5.2). Celui-ci suffit pour en déduire son numérateur $N(X)$. Pour ce faire, on utilise la formule $1/(a-X) = \sum_{j \ge 0} a^{-j-1} X^j$ valable dans $\mathbb{A}[[X]]$ pour tout $a \in \mathbb{A}$ inversible. On obtient :

$$\sum_{i=0}^{n-1} \frac{c_i}{X - q^i} \quad \text{mod} \quad X^n = - \sum_{i=0}^{n-1} \left( \sum_{j=0}^{n-1} c_i q^{-i(j+1)} X^j \right) = - \sum_{j=0}^{n-1} C(q^{-j-1}) X^j,$$

où $C(X)$ est le polynôme $\sum_{i=0}^{n-1} c_i X^i$. Or, les points $q^{-1}, q^{-2}, \dots, q^{-n}$ étant en progression géométrique, on sait évaluer $C$ sur ces points en $O(M(n))$.

**Exercice 5.4** Montrer que tout l'arbre des sous-produits associé aux points $a_i = q^i$, $i = 0, \dots, n-1$, peut être calculé en $M(n) + O(n \log n)$ opérations dans $\mathbb{A}$.

Une approche évaluation-interpolation reposant sur la Proposition 5.10 permet d'obtenir un bon algorithme pour la multiplication des matrices polynoniales.
**Corollaire 5.11** Soit $\mathbb{K}$ un corps possédant au moins 2d éléments. Il est possible de multiplier deux matrices polynomiales de $M_n(\mathbb{K}[X])$ de degré au plus $d$ en $\mathbb{K}(n, d) = O(d \mathbb{K}(n) + n^2 \mathbb{K}(d))$ opérations dans $\mathbb{K}$.

Exercices

**Exercice 5.5** — Calcul rapide de sommes de puissances. Soient $x_1, x_2, \dots, x_n$ des éléments de $\mathbb{K}$ et soit $P_k = \sum_{i=1}^n x_i^k$ pour $k \ge 1$. Montrer que, à partir de la donnée des $x_1, \dots, x_n$, il est possible de calculer :

1. tous les $P_k$ dont les indices $k \le n$ sont une puissance de 2 en $O(n \log n)$ opérations dans $\mathbb{K}$.

2. tous les $P_1, \dots, P_n$ en $O(M(n) \log n)$ opérations dans $\mathbb{K}$.

**Exercice 5.6** Estimer la complexité de l’algorithme direct pour évaluer un polynôme de degré au plus $n$ et ses premières $n$ dérivées en un point. Imaginer un algorithme de complexité quasi-linéaire résolvant ce problème.

**Exercice 5.7** Soient $a, b, c, d \in \mathbb{A}$. Montrer qu'on peut évaluer tout polynôme de degré au plus $n$ en $\{ba^{2i} + ca^i + d, \ 0 \le i < n\}$, en $O(M(n))$ opérations dans $\mathbb{A}$.

**Exercice 5.8** Soit $\mathbb{K}$ un corps, $m_1, m_2 \in \mathbb{K}[X]$ des polynômes premiers entre eux, de degrés au plus $n$, et $v_1, v_2 \in \mathbb{K}[X]$ de degrés strictement inférieurs à $n$. Donner un algorithme de complexité $O(M(n) \log n)$ qui calcule $f \in \mathbb{K}[X]$ de degré strictement inférieur à $2n$ tel que $f = v_1 \bmod m_1$ et $f = v_2 \bmod m_2$. [Indication : utiliser un pgcd étendu.]

En déduire un algorithme de type « diviser pour régner » pour l’interpolation polynomiale en degré $n$, de complexité $O(M(n) \log^2 n)$.

**Exercice 5.9** — Multiplication de polynômes à deux variables par substitution de Kronecker. Soient $f$ et $g$ dans $\mathbb{K}[X, Y]$ de degrés au plus $d_X$ en $X$ et au plus $d_Y$ en $Y$.

1. Montrer qu'il est possible de calculer tous les coefficients de $h = fg$ avec un algorithme de complexité $O(M(d_X d_Y))$ opérations dans $\mathbb{K}$. Utiliser la substitution $X \leftarrow Y^{2d_Y+1}$ pour se ramener à une multiplication de polynômes à une variable.

2. Améliorer ce résultat en donnant un schéma d’évaluation-interpolation qui permette le calcul des coefficients de $h$ en $O(d_X M(d_Y) + d_Y M(d_X))$ opérations de $\mathbb{K}$ si $\mathbb{K}$ est assez grand.

Notes

Pour évaluer un polynôme $P(X) = p_{n-1}X^{n-1} + \cdots + p_0$ en un seul point $a$, le schéma de Horner $P(a) = (\cdots(p_{n-1}a + p_{n-2})a + p_{n-3})a + \cdots + p_0)$ minimise le nombre d’opérations (additions ou multiplications). Cet algorithme effectue $n-1$ additions et
$n-1$ multiplications dans $\mathbb{A}$, et on ne peut pas faire mieux en général (ceci, si les coefficients du polynôme sont algébriquement indépendants — c'est-à-dire s'ils se comportent comme des indéterminées). Ce résultat d'optimalité a été conjecturé par Ostrowski [Ost54] et prouvé par Pan [Pan66]. Par contre, pour un polynôme donné, il peut être possible de faire mieux : l'évaluation de $P(X) = X^{n-1}$ se fait en temps logarithmique en $n$.

Une forme particulière du théorème des restes chinois a été enoncée il y a plus de 2000 ans par le mathématicien chinois Sun Tsu. Le traitement moderne est dû à Gauss [Gau66]. Ce théorème admet une formulation très générale, en termes d'idéaux d'anneaux non nécessairement commutatifs : Si $R$ est un anneau et $I_1, \dots, I_r$ des idéaux (à gauche) de $R$ mutuellement premiers (c'est-à-dire $I_i + I_j = R$ pour $i < j$), alors l'anneau quotient $R/\cap_i I$ est isomorphe à l'anneau produit $\prod R/I_i$ via l'isomorphisme $x \mod I \mapsto (x \mod I_1, \dots, x \mod I_r)$. La formule d'interpolation de Lagrange est due à Waring [War79]. Les avantages pratiques de l'arithmétique modulaire et les premières applications en informatique du théorème des restes chinois ont été mis en évidence dans les années 1950 par Svoboda et Valach [SV55], et indépendamment par Garner [Gar59] et Takahasi et Ishibashi [TI61].

Le premier algorithme sous-quadratique, de complexité arithmétique $O(n^{1,91})$, pour l'évaluation multipoint est dû à Borodin et Munro [BM71] et repose sur une approche *pas de bébés / pas de géants* exploitant la multiplication sous-cubique des matrices de Strassen [Str69]. Horowitz [Hor72] a étendu ce résultat à l'interpolation. La Proposition 5.5 calcule efficacement les fonctions symétriques élémentaires de $n$ éléments $a_0, \dots, a_{n-1}$ de $\mathbb{A}$; elle est due également à Horowitz [Hor72]. S'inspirant de la FFT, Fiduccia [Fid72b] eut l'idée d'un algorithme pour l'évaluation multipoint à base de division polynomiale récursive. Ne disposant pas de division de complexité sous-quadratique, son algorithme récursif reste quadratique. Borodin et Moenck corrigent ce point dans deux articles successifs [BM74; MB72], reposant sur des algorithmes quasi-optimaux pour la division [MB72; Str73]. Montgomery [Mon92] améliore la constante dans la complexité $O(M(n) \log n)$ de l'évaluation multipoint, sous l'hypothèse que la FFT est utilisée pour la multiplication polynomiale. L'algorithme d'évaluation sur une suite géométrique exposé dans ce chapitre provient des travaux de Rabiner, Schafer, Rader [RSR69], et Bluestein [Blu70]. À l'origine, cet algorithme a permis de montrer qu'une $DFT_\omega$, pour $\omega$ racine primitive $n$-ième de l'unité, peut être effectuée en $O(n \log n)$ opérations même si $n$ n'est pas une puissance de 2. L'algorithme d'interpolation est une adaptation de celui de Mersereau [Mer74]. Les meilleures constantes actuellement connues dans les $O(\cdot)$ du Théorème 5.1 sont obtenus à l'aide du principe de transposition de Tellegen [BLS03; BS05].

On connaît peu de choses sur la complexité intrinsèque de l'évaluation multipoint et de l'interpolation en degré $n$. Strassen [Str73] a prouvé une borne inférieure de type $n \log n$. Shoup et Smolensky [SS97] ont montré que ces bornes restent essentiellement vraies même dans un modèle où des précalculs en quantité illimitée sur les points d'évaluation sont permis. Cependant, l'optimalité des meilleurs algorithmes connus, de complexité $O(n \log^2 n)$, reste un problème ouvert. Même dans le cas particulier où les points d'évaluation sont en progression arithmétique, on ne dispose d'aucun algorithme d'évaluation multipoint de complexité $O(M(n))$, ni d'une preuve qu'un tel algorithme n'existe pas. De même, pour des points quelconques $a_0, \dots, a_{n-1}$ on ne
connaît pas d'algorithmes de complexité $O(M(n))$ pour calculer le polynôme $\prod_i(X - a_i)$. Notons toutefois que de tels algorithmes existent dans les cas particuliers où les $a_i$ forment une progression arithmétique ou géométrique [BS05]. L'Exercice 5.8 provient de [HH71]; historiquement, ce fut le premier algorithme rapide pour les restes chinois. Les Exercices 5.6 et 5.7 sont tirés de [ASU75].

Bibliographie

ASU75 AHO, A. V., K. STEIGLITZ et J. D. ULLMAN (1975). « Evaluating polynomials at fixed sets of points ». In : *SIAM Journal on Computing*, vol. 4, n°4, p. 533-539.

BLS03 BOSTAN, Alin, Grégoire LECERF et Éric SCHOST (2003). « Tellegen's principle into practice ». In : *ISSAC'03 : International Symposium on Symbolic and Algebraic Computation*. Éd. par J. R. SENDRA. ACM Press, p. 37-44.

Blu70 BLUESTEIN, L. I. (1970). « A linear filtering approach to the computation of the discrete Fourier transform ». In : *IEEE Trans. Electroacoustics*, vol. AUS-18, p. 451-455.

BM71 BORODIN, A. et I. MUNRO (1971). « Evaluating polynomials at many points ». In : *Information Processing Letters*, vol. 1, n°2, p. 66-68.

BM74 BORODIN, A. et R. T. MOENCK (1974). « Fast modular transforms ». In : *Journal of Computer and System Sciences*, vol. 8, n°3, p. 366-386.

BS05 BOSTAN, Alin et Éric SCHOST (2005). « Polynomial evaluation and interpolation on special sets of points ». In : *Journal of Complexity*, vol. 21, n°4, p. 420-446.

Fid72b FIDUCCIA, C. M. (1972). « Polynomial evaluation via the division algorithm : the fast Fourier transform revisited. » In : *STOC'72 : ACM Symposium on Theory of Computing*, p. 88-93.

Gar59 GARNER, Harvey L. (1959). « The residue number system ». In : *IRE-AIEE-ACM'59 Western joint computer conference*. San Francisco, California : ACM, p. 146-153.

Gau66 GAUSS, Carl Friedrich (1966). *Disquisitiones arithmeticae*. Translated into English by Arthur A. Clarke, S. J. Yale University Press.

HH71 HEINDEL, L. E. et E. HOROWITZ (1971). « On decreasing the computing time for modular arithmetic ». In : *12th symposium on switching and automata theory*. IEEE, p. 126-128.

Hor72 HOROWITZ, E. (1972). « A fast method for interpolation using preconditioning ». In : *Information Processing Letters*, vol. 1, n°4, p. 157-163.

MB72 MOENCK, R. T. et A. BORODIN (1972). « Fast modular transforms via division ». In : *Thirteenth Annual IEEE Symposium on Switching and Automata Theory*, p. 90-96.

Mer74 MERSEREAU, Russell M. (1974). « An algorithm for performing an inverse chirp z-transform ». In : *IEEE Transactions on Audio and Electroacoustics*, vol. ASSP-22, n°5, p. 387-388.

Mon92 MONTGOMERY, P. L. (1992). « An FFT extension of the elliptic curve method of factorization ». Thèse de doctorat. University of California, Los Angeles CA.
Ost54 OSTROWSKI, A. M. (1954). « On two problems in abstract algebra connected with Horner's rule ». In: *Studies in mathematics and mechanics presented to Richard von Mises*. NY: Academic Press Inc., p. 40–48.

Pan66 PAN, V. Ya. (1966). « Methods of computing values of polynomials ». In: *Russian Mathematical Surveys*, vol. 21, n°1, p. 105–136.

RSR69 RABINER, L. R., R. W. SCHAFER et C. M. RADER (1969). « The chirp z-transform algorithm and its application ». In: *Bell System Technical Journal*, vol. 48, p. 1249–1292.

SS97 SHOUP, Victor et Roman SMOLENSKY (1996/97). « Lower bounds for polynomial evaluation and interpolation problems ». In: *Computational Complexity*, vol. 6, n°4, p. 301–311.

Str69 STRASSEN, V. (1969). « Gaussian elimination is not optimal ». In: *Numerische Mathematik*, vol. 13, p. 354–356.

Str73 — (1972/73). « Die Berechnungskomplexität von elementarsymmetrischen Funktionen und von Interpolationskoeffizienten ». In: *Numerische Mathematik*, vol. 20, p. 238–251.

SV55 SVOBODA, A. et M. VALACH (1955). « Operátorové obvody (Operational circuits) ». In: *Stroje na Zpracování Informací (Information Processing Machines)*, vol. 3, p. 247–295.

TI61 TAKAHASI, H. et Y. ISHIBASHI (1961). « A new method for exact calculation by a digital computer ». In: *Information Processing in Japan*, p. 28–42.

War79 WARING, E. (1779). « Problems concerning interpolations ». In: *Philosophical Transactions of the Royal Society of London*, vol. 59, p. 59–67.
<figure><img src="image_19.png" /></figure>
# 6. Pgcd et résultant

## Résumé

Les calculs de pgcd sont cruciaux pour la simplification des fractions, qu'il s'agisse de fractions d'entiers ou de fractions de polynômes. Les algorithmes efficaces de factorisation de polynômes reposent par ailleurs de manière essentielle sur des calculs de pgcd. Comme pour la multiplication, les algorithmes efficaces de pgcd sont plus complexes dans le cas des entiers que dans le cas des polynômes, à cause des retenues, et nous n'en détaillons donc que la version polynomiale. Dans les définitions et résultats de nature moins algorithmique, nous utilisons le cadre algébrique des anneaux euclidiens qui permet de traiter simultanément toutes les applications que nous avons en vue. L'algorithme d'Euclide pour le calcul du pgcd est probablement le plus vieil algorithme non-trivial encore en utilisation. Il est présenté en Section 6.1. Dans le cas des polynômes de $K[X]$, sa complexité est quadratique en nombre d'opérations dans le corps $K$. Le résultant est également lié à l'algorithme d'Euclide, ses propriétés et son calcul sont présentés en Section 6.2. En outre, une technique dite « des sous-résultants » permet de réduire l'explosion de la taille des coefficients intermédiaires dont souffre l'algorithme d'Euclide pour le pgcd dans $Q[X]$. Le chapitre se termine en Section 6.3 sur un algorithme de complexité quasi-optimale pour le pgcd dans $K[X]$, exploitant la multiplication polynomiale rapide.

*Dans ce chapitre, $\mathbb{A}$ désignera toujours un anneau intègre.*

## 6.1 Algorithme d'Euclide

### Le pgcd

Dans tout anneau intègre $\mathbb{A}$, on appelle *plus grand diviseur commun* (pgcd) de $A \in \mathbb{A}$ et de $B \in \mathbb{A}$ tout $G \in \mathbb{A}$ qui divise A et B et tel que tout diviseur de A et de B divise
aussi G. Cependant, dans ce degré de généralité, la notion est délicate : il existe des anneaux intègres contenant des éléments sans pgcd.

**Exercice 6.1** Montrer que dans $\mathbb{A} = \mathbb{Z}[i\sqrt{3}]$, les éléments $\mathbb{A} = 4$ et $\mathbb{B} = 2 + 2i\sqrt{3}$ admettent $G_1 = 2$ et $G_2 = 1 + i\sqrt{3}$ comme diviseurs, mais n'admettent aucun pgcd. ■

On appelle traditionnellement *anneau à pgcd* un anneau intègre dans lequel tous couple d'éléments admet un pgcd. Algébriquement, le bon cadre pour parler de pgcd est celui des anneaux factoriels ; d'un point de vue algorithmique, une classe d'anneaux à pgcd bien adaptée aux calculs est celle des anneaux euclidiens. Contrairement à l'usage, nous notons dans ce chapitre pgcd(A, B) tout pgcd de A et de B sans faire de choix de normalisation parmi les pgcd de A et de B. L'algorithme d'Euclide présenté dans cette section permet le calcul du pgcd dans $\mathbb{Z}$ ou dans l'anneau $\mathbb{K}[X]$, où $\mathbb{K}$ est un corps. Il repose sur l'existence d'une *division euclidienne* dans ces anneaux. Un anneau intègre $\mathbb{A}$ est appelé *anneau euclidien* s'il existe une fonction de taille $d : \mathbb{A} \to \mathbb{N} \cup \{-\infty\}$ telle que pour tout $a \in \mathbb{A}$, $b \in \mathbb{A} \setminus \{0\}$, il existe $q$ et $r$ dans $\mathbb{A}$ avec

$$a = qb + r, \quad d(r) < d(b).$$

En particulier, $d(0) < d(b)$ dès que $b$ est non nul. Dans le cas des entiers, la valeur absolue fournit une telle fonction $d$, et le degré remplit ce rôle pour les polynômes. La notation $r := a \bmod b$ sert dans les deux cas à définir $r$ à partir de $a$ et de $b$. L'entier $r$ est appelé *reste* de la division euclidienne de $a$ par $b$.

<figure><img src="image_20.png" /></figure>

Un élément d'un anneau est dit irréductible si les produits qui lui sont égaux font tous intervenir un élément inversible. Tout anneau euclidien est *principal* (c'est-à-dire composé des multiples d'un élément donné de l'anneau), et en particulier il est *factoriel* (c'est-à-dire que tout élément non nul y a une factorisation unique en irréductibles, à l'ordre près de la factorisation et à des multiplications près des facteurs par des inversibles de l'anneau). Le pgcd se lit aussi sur les factorisations. Si deux éléments A et B d'un anneau factoriel A se factorisent sous la forme

$$A = a m_1^{k_1} \cdots m_s^{k_s}, \quad B = b m_1^{\ell_1} \cdots m_s^{\ell_s},$$

où $a$ et $b$ sont inversibles dans $\mathbb{A}$, les $m_i$ sont irréductibles et les $k_i$ et $\ell_i$ sont des entiers positifs ou nuls, montrer qu'alors un pgcd de A et B est

$$G = m_1^{\min(k_1, \ell_1)} \cdots m_s^{\min(k_s, \ell_s)}.$$

Autrement dit, on rassemble dans $G$ les facteurs irréductibles communs de $A$ et $B$, avec la plus grande multiplicité possible. La factorisation dans $\mathbb{Z}$ ou dans $\mathbb{K}[X]$ est plus coûteuse que le pgcd et cette propriété n'est donc pas utilisée pour le calcul du pgcd.

## Calcul du pgcd

Étant donnés A et B dans un anneau euclidien $\mathbb{A}$, l'Algorithme 6.1 d'Euclide calcule une suite de restes successifs dont la taille décroît, jusqu'à atteindre le pgcd. La terminaison provient de la décroissance stricte de la taille à chaque étape. La correction de cet algorithme se déduit de la relation

$$\mathrm{pgcd}(F, G) = \mathrm{pgcd}(H, G) \quad \text{pour} \quad H := F \bmod G.$$
**Entrée** A et B dans A.

**Sortie** Un pgcd de A et B.

1. $R_0 := A; R_1 := B; i := 1$.

2. Tant que $R_i$ est non nul, faire :

$$R_{i+1} := R_{i-1} \mod R_i$$

$$i := i + 1.$$

3. Renvoyer $R_{i-1}$.

### Algorithme 6.1 – Algorithme d'Euclide.

Par récurrence, il s'ensuit que $\text{pgcd}(F, G) = \text{pgcd}(R_i, R_{i+1})$ pour tout $i$. Si en outre $R_{i+1}$ est nul, alors $\text{pgcd}(R_i, R_{i+1}) = R_i$, ce qui prouve la correction.

**Exemple 6.1** Soient $A = X^4 - 13X^3 + 2X^2 - X - 1$ et $B = X^2 - X - 1$ dans $\mathbb{Q}[X]$. La suite des restes est :

$$ \begin{aligned} R_0 &= X^4 - 13X^3 + 2X^2 - X - 1, & R_1 &= X^2 - X - 1, \\ R_2 &= -22X - 10, & R_3 &= -41/121, & R_4 &= 0, \end{aligned} $$

de sorte que $-41/121$ est un pgcd de A et B et donc 1 aussi.

**Théorème 6.1** L'algorithme d'Euclide calcule un pgcd de deux polynômes A et B dans $\mathbb{K}[X]$ en $O(\deg A \deg B)$ opérations dans $\mathbb{K}$.

*Démonstration.* La correction a été prouvée dans le cas général. Pour l'étude de complexité, nous supposons d'abord que $\deg A \ge \deg B$. D'après le Chapitre 3 (page 82), le calcul naïf de $P \mod Q$ peut être effectué en $(2 \deg Q + 1)(\deg P - \deg Q + 1) + 1$ opérations de $\mathbb{K}$. Il s'ensuit que le coût de l'algorithme d'Euclide est borné par la somme des $(2 \deg(R_i) + 1)(\deg R_{i-1} - \deg R_i + 1) + 1$, pour $i \ge 1$. Tous les $\deg(R_i)$ pour $i \ge 1$ sont majorés par $\deg B$, de sorte que le coût est borné par $(2 \deg B + 1) \sum_{i \ge 1} (\deg R_{i-1} - \deg R_i + 1) + \deg B \le (2 \deg B + 2)(\deg R_0 + \deg B) = O(\deg A \deg B)$.

Si le degré de B est supérieur à celui de A, la première étape ne coûte pas d'opération arithmétique, et la borne ci-dessus s'applique au reste du calcul. ■

**Exercice 6.2** Énoncer et prouver l'analogue entier du résultat du Théorème 6.1. ■

La borne de complexité quadratique reflète bien le comportement de l'algorithme : dans une exécution typique les quotients ont degré 1 à chaque itération, les degrés des restes successifs diminuent de 1 à chaque itération et leur calcul est linéaire; le nombre de coefficients intermédiaires calculés est quadratique et ils sont calculés aussi efficacement que possible par un algorithme qui les calcule tous.
## Pgcd étendu et inversion modulaire

### Relation de Bézout

Une relation de la forme $G = UA + VB$ qui donne le pgcd $G$ de deux éléments $A$ et $B$ de $\mathbb{A}$ avec deux cofacteurs $U$ et $V$ dans $\mathbb{A}$ est appelée une *relation de Bézout*. L'algorithme d'Euclide étendu est une modification légère de l'algorithme d'Euclide qui calcule non seulement le pgcd, mais aussi une relation de Bézout particulière,

$$UA + VB = G, \text{ avec } d(UG) < d(B) \text{ et } d(VG) < d(A). \qquad (6.1)$$

Une fois le pgcd $G$ choisi, les cofacteurs $U$ et $V$ sont rendus uniques par la contrainte sur les tailles. Dans de nombreuses applications, c'est davantage de ces cofacteurs $U$ et $V$ que du pgcd lui-même dont on a besoin. On parle alors de calcul de pgcd étendu. Ce calcul intervient par exemple de manière importante dans la manipulation des nombres algébriques (Exemple 6.3 ci-dessous), dans le calcul dans des algèbres quotient $\mathbb{K}[X]/(P)$ (Exemples 6.2 et 6.5 ci-dessous, cf. aussi la Section 4.2 en page 84), et dans le développement rapide des séries algébriques (Chapitre 14). Il intervient également dans d'autres algorithmes qui ne seront pas abordés dans cet ouvrage, par exemple pour le calcul des « restes chinois rapides », pour la factorisation de polynômes dans $\mathbb{Z}[X]$ ou $\mathbb{Q}[X]$ par des techniques de type « Hensel » (Chapitre 21), ou pour l'intégration symbolique (à la Hermite).

Mais, avant tout, les calculs de pgcd étendu permettent d'effectuer des *inversions modulaires*. Lorsque $A$ et $B$ sont premiers entre eux, (c'est-à-dire si $G$ est un élément inversible de $\mathbb{A}$), alors l'élément $V$ est un inverse de $B$ modulo $A$.

**Exemple 6.2** La relation de Bézout pour $a + bX$ et $1 + X^2$ dans $\mathbb{R}[X]$ s'écrit :

$$(a - bX)(a + bX) + b^2(1 + X^2) = a^2 + b^2.$$

L'inverse de $B = a + bX$ modulo $A = 1 + X^2$ vaut donc

$$V = \frac{a - bX}{a^2 + b^2}.$$

Puisque le corps des complexes s'obtient par le quotient $\mathbb{R}[X]/(X^2+1)$, cette relation n'est autre que la formule d'inversion familière

$$\frac{1}{a + ib} = \frac{a - ib}{a^2 + b^2} \quad \text{où } i = \sqrt{-1}.$$

Si $A \in \mathbb{A}$ est irréductible, alors $A/(A)$ est un corps, et le calcul de la relation de Bézout permet d'y effectuer la division : si $B \neq 0 \bmod A$, alors $N/B = NV \bmod A$.

**Exemple 6.3** La « simplification » de la fraction rationnelle

$$r = \frac{\phi^4 - \phi + 1}{\phi^7 - 1} = \frac{25\phi - 34}{29}$$
où $\phi$ est le « nombre d'or », défini par $A(\phi) = 0$ pour $A(X) = X^2 - X - 1$, est obtenue par les trois calculs suivants :

- calcul du reste $U = 13X + 7$ de la division euclidienne de $X^7 - 1$ par A,

- détermination de la relation de Bézout (13X - 20)U - 169A = 29,

- calcul du reste $V = 25X - 34$ de la division euclidienne de $(13X - 20)(X^4 - X + 1)$ par A.

Plus généralement, dans le cas où $A \in \mathbb{K}[X]$ est un polynôme irréductible, ces calculs montrent que le corps $\mathbb{K}(\alpha)$ des fractions rationnelles en $\alpha$ racine de A est un espace vectoriel dont une base est $1, \alpha, \alpha^2, ..., \alpha^{\text{deg } A-1}$. Les algorithmes d'Euclide et d'Euclide étendu fournissent un moyen de calcul dans cette représentation. Il est ainsi possible de manipuler de manière exacte les racines d'un polynôme de degré arbitraire sans « résoudre ».

**Exemple 6.4** Le même calcul qu'à l'Exemple 6.3 fonctionne sur des entiers :

$$r = \frac{25}{33} \equiv 5 \quad \text{mod} \quad 7$$

se déduit de $33 \mod 7 = 5, 3 \times 5 - 2 \times 7 = 1, 3 \times 25 \mod 7 = 5$.

**Exemple 6.5** Lorsque l'élément $A \in \mathbb{A}$ n'est pas irréductible, $\mathbb{A}/(A)$ n'est pas un corps. Cependant, les éléments inversibles peuvent y être inversés par le même calcul que ci-dessus. Lorsqu'un élément B non nul n'est pas inversible, la relation de Bézout se produit avec un G différent de 1. Il est alors possible de tirer parti de cette information (un facteur de A) en scindant le calcul d'une part sur G et d'autre part sur B/G.

**Exercice 6.3** Soit $p$ un nombre premier et soit $a$ un entier non multiple de $p$. Donner un algorithme pour le calcul de l'inverse de $a$ modulo $p$ utilisant le petit théorème de Fermat ($a^{p-1} = 1 \mod p$) et l'exponentiation binaire. Analyser la complexité binaire de cet algorithme et comparer avec l'approche par Euclide étendu. ■

## Calcul des cofacteurs

L'idée-clé est de suivre pendant l'algorithme d'Euclide la décomposition de chacun des restes $R_i$ sur A et B. Autrement dit, pour tout $i$, l'algorithme calcule des éléments $U_i$ et $V_i$ tels que

$$U_i A + V_i B = R_i,$$

dont les tailles sont bien contrôlées. Pour $i = 0$, il suffit de poser $U_0 = 1, V_0 = 0$, ce qui correspond à l'égalité $1 \cdot A + 0 \cdot B = A = R_0$. Pour $i = 1$, l'égalité $0 \cdot A + 1 \cdot B = B = R_1$ est obtenue avec $U_1 = 0, V_1 = 1$. Ensuite, la division euclidienne $R_{i-1} = Q_i R_i + R_{i+1}$ donne la relation

$$R_{i+1} = R_{i-1} - Q_i R_i = (U_{i-1} - Q_i U_i)A + (V_{i-1} - Q_i V_i)B,$$

qui pousse à définir $U_{i+1}$ par $U_{i-1} - Q_i U_i$ et $V_{i+1}$ par $V_{i-1} - Q_i V_i$. À partir de cette relation, une preuve par récurrence montre que les conditions de tailles de la relation
**Entrée** A et B dans A.

**Sortie** Un pgcd G de A et B, et les cofacteurs U et V de l'identité de Bézout UA + VB = G.

1. $R_0 := A; U_0 := 1; V_0 := 0; R_1 := B; U_1 := 0; V_1 := 1; i := 1$.

2. Tant que $R_i$ est non nul, faire :

a. Calculer le quotient $Q_i$ et le reste $R_{i+1}$ de la division euclidienne de $R_{i-1}$ par $R_i$.

b. $U_{i+1} := U_{i-1} - Q_i U_i$ ; $V_{i+1} := V_{i-1} - Q_i V_i$.

c. $i := i + 1$.

3. Renvoyer $R_{i-1}, U_{i-1}, V_{i-1}$.

Algorithme 6.2 – Algorithme d'Euclide étendu.

de Bézout sont satisfaites. La méthode est résumée dans l'Algorithme 6.2. À nouveau, dans le cas des polynômes ou des entiers, la complexité est quadratique.

## 6.2 Résultant

### Matrice de Sylvester

L'algorithme d'Euclide pour deux polynômes A, B à une variable est étroitement relié à la matrice de Sylvester. Si

$$A = a_m X^m + \cdots + a_0, \quad (a_m \neq 0), \quad \text{et} \quad B = b_n X^n + \cdots + b_0, \quad (b_n \neq 0),$$

appartiennent à $\mathbb{A}[X]$, cette matrice, carrée de taille $m+n$, s'écrit

$$\text{Syl}(A, B) = \begin{pmatrix} a_m & \cdots & a_1 & a_0 \\ & \ddots & & \vdots \\ b_n & \cdots & b_1 & b_0 \\ & \ddots & & \vdots \\ & & b_n & b_{n-1} & \cdots & b_0 \end{pmatrix}$$

avec les $n$ premières lignes contenant les coefficients de A et les $m$ suivantes contenant les coefficients de B. La transposée de cette matrice représente l'application linéaire

$$\varphi_{A,B} : (U, V) \mapsto AU + BV,$$

en se limitant aux paires $(U, V)$ telles que $\deg U < \deg B$ et $\deg V < \deg A$, et en utilisant la base $(X^{n-1}, 0), (X^{n-2}, 0), \dots, (1, 0), (0, X^{m-1}), (0, X^{m-2}), \dots, (0, 1)$ au départ et la base $(X^{m+n-1}, X^{m+n-2}, \dots, 1)$ à l'arrivée. Les combinaisons linéaires des lignes de la matrice de Sylvester donnent donc les coefficients des polynômes qui peuvent être obtenus comme image.

Lorsque $\mathbb{A}$ est un corps, au vu de la relation de Bézout (6.1), le pgcd de A et B peut être obtenu ainsi et il est caractérisé comme l'élément unitaire de plus petit
degré qui peut l'être. Ceci permet de le lire sur une forme échelonnée en ligne de la matrice Syl(A, B) (cf. Définition 8.3, page 170). Qui plus est, le rang de Syl(A, B) vaut la dimension de l'image de l'application $\varphi_{A,B}$, qui est égale à $m+n-\text{deg pgcd}(A, B)$. Par le théorème du rang, il s'ensuit que le degré du pgcd(A, B) est égal à la dimension du noyau de Syl(A, B), c'est-à-dire au nombre de lignes nulles de sa forme échelonnée.

**Proposition 6.2** Soient A, B deux polynômes de $K[X]$ où $K$ est un corps. La forme échelonnée en ligne de la matrice de Sylvester Syl(A, B) contient sur sa dernière ligne non nulle les coefficients d'un pgcd de A et B. En outre, le nombre de lignes nulles de cette matrice échelonnée est le degré du pgcd de A et de B.

**Exemple 6.6** Pour les polynômes $A = X^4 - X^3 - 7X^2 + 2X + 3$ et $B = X^3 - 4X^2 + 2X + 3$, la matrice de Sylvester et sa forme échelonnée en ligne sont :

$$ \begin{pmatrix} 1 & -1 & -7 & 2 & 3 & 0 & 0 \\ 0 & 1 & -1 & -7 & 2 & 3 & 0 \\ 0 & 0 & 1 & -1 & -7 & 2 & 3 \\ 1 & -4 & 2 & 3 & 0 & 0 & 0 \\ 0 & 1 & -4 & 2 & 3 & 0 & 0 \\ 0 & 0 & 1 & -4 & 2 & 3 & 0 \\ 0 & 0 & 0 & 1 & -4 & 2 & 3 \end{pmatrix} \quad \text{et} \quad \begin{pmatrix} 1 & -1 & -7 & 2 & 3 & 0 & 0 \\ 0 & 1 & -1 & -7 & 2 & 3 & 0 \\ 0 & 0 & 1 & -1 & -7 & 2 & 3 \\ 0 & 0 & 0 & -14 & 45 & -3 & -18 \\ 0 & 0 & 0 & 0 & -\frac{5}{7} & \frac{12}{7} & \frac{9}{7} \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{10} & -\frac{3}{10} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} $$

La dernière ligne non nulle donne le pgcd $X - 3$.

**Définition 6.1** Le résultant de A et B est le déterminant de la matrice de Sylvester Syl(A, B). Il est noté $\text{Res}(A, B)$, ou $\text{Res}_X(A, B)$ si l'on veut insister sur l'élimination de la variable X.

L'utilisation des résultats pour l'élimination demande de travailler dans ce chapitre avec des matrices à coefficients dans un anneau commutatif qui n'est pas nécessairement un corps. Les définitions et propriétés du déterminant persistent dans ce cas (c'est une forme multilinéaire alternée, le déterminant de la matrice identité vaut 1, on peut le développer par rapport à une ligne ou une colonne, il ne change pas par addition de combinaisons linéaires de lignes ou de colonnes, et s'exprime comme somme sur les permutations).

Le résultant peut être vu comme une condition de cohérence pour le système formé par les deux polynômes A et B.

**Corollaire 6.3** Soient A et B deux polynômes de $K[X]$. Alors A et B sont premiers entre eux si et seulement si $\text{Res}(A, B) \neq 0$.

*Démonstration.* Le déterminant d'une matrice carrée se lit sur la forme échelonnée en ligne en faisant le produit des éléments diagonaux. Il est non nul si et seulement si tous les éléments diagonaux le sont et dans ce cas un pgcd non nul est sur la dernière ligne. À l'inverse, s'il existe un pgcd G non nul, alors les polynômes $X^kG$ montrent
<figure><img src="image_21.png" /><figcaption>FIGURE 6.3 – Le résultant calcule des projections.</figcaption></figure>

l'existence de lignes avec un coefficient de tête non nul dans chaque colonne de la forme échelonnée. ■

**Exemple 6.7** Le discriminant d'un polynôme $A = a_m X^m + \dots + a_1 X + a_0$ de degré $m$ à coefficients dans un corps $K$ est le polynôme disc$(A)$ défini par

$$ \mathrm{Res}(A, A') = (-1)^{m(m-1)/2} a_m \mathrm{disc}(A). $$

Il s'annule lorsque ces deux polynômes ont une racine commune dans une clôture algébrique de $K$, c'est-à-dire, en caractéristique nulle, lorsque $A$ a une racine multiple.

**Exercice 6.4** Calculer le discriminant de $aX^2 + bX + c$ en prenant le déterminant de la matrice de Sylvester. ■

## Applications du résultant

### Calculs de projections

Algébriquement, la « résolution » d'un système polynomial se ramène souvent à une question d'élimination. Lorsqu'elle est possible, l'élimination successive des variables amène le système d'entrée sous une forme triangulaire. De proche en proche, la résolution d'un tel système se réduit alors à la manipulation de polynômes à une variable, pour lesquels compter, isoler, ..., les solutions est bien plus facile. Géométriquement, l'élimination correspond à une projection. Dans cette section nous montrons comment le résultant de deux polynômes permet de traiter ce genre de problèmes, pour les systèmes de deux équations en deux inconnues. Le schéma général est représenté en Figure 6.3. Cette opération est à la base d'une technique plus générale pour un nombre arbitraire d'équations et d'inconnues, la résolution géométrique, qui sera présentée au Chapitre 28.
<figure><img src="image_22.png" /><figcaption>FIGURE 6.4 – Courbes de l'Exemple 6.8.</figcaption></figure>

**Exemple 6.8** Les deux courbes de la Figure 6.4 ont pour équations

$$A = (X^2 + Y^2)^3 - 4X^2Y^2 = 0, \quad B = X^2(1+Y) - (1-Y)^3 = 0.$$

Ces deux polynômes sont irréductibles et leur pgcd, qui vaut 1, ne renseigne pas sur leurs racines communes. Les résultats par rapport à X et Y donnent en revanche des polynômes qui s'annulent sur les coordonnées de ces points d'intersection :

$$\begin{aligned} \operatorname{Res}_X(A, B) &= (4Y^7 + 60Y^6 - 152Y^5 + 164Y^4 - 95Y^3 + 35Y^2 - 9Y + 1)^2, \\ \operatorname{Res}_Y(A, B) &= 16X^{14} + 6032X^{12} - 1624X^{10} + 4192X^8 - 815X^6 - 301X^4 - 9X^2 + 1. \end{aligned}$$

Il n'y a que 4 points d'intersection visibles sur la Figure 6.4, les 10 autres ont au moins une coordonnée qui n'est pas réelle. Le caractère bicarré du résultant en Y provient de la symétrie de la figure par rapport à l'axe des Y. C'est pour cette même raison que le premier résultant est un carré.

Le résultat général est le suivant.

**Proposition 6.4** Soient $A = a_m Y^m + \dots$ et $B = b_n Y^n + \dots$ où les coefficients $a_i$ et $b_j$ sont dans $K[X]$ où $K$ est un corps algébriquement clos. Alors les racines du polynôme $\operatorname{Res}_Y(A, B) \in K[X]$ sont d'une part les abscisses des solutions du système $A = B = 0$, d'autre part les racines des termes de tête $a_m$ et $b_n$.

*Démonstration.* La preuve repose sur des résultats de la section suivante, qui sont indépendants de cette proposition.

Tout d'abord, les racines des termes de tête $a_m$ et $b_n$ sont racines du résultant d'après la formule de Poisson (Théorème 6.5 ci-dessous).
<figure><img src="image_23.png" /><figcaption>FIGURE 6.5 – Courbes de l’Exemple 6.9.</figcaption></figure>

Ensuite, d'après la Proposition 6.7 ci-dessous, il existe U et V dans $K[X, Y]$ tels que

$$ \mathrm{Res}_Y(A, B) = UA + VB. $$

En évaluant cette identité en $(x, y)$ tel que $A(x, y) = B(x, y) = 0$, on voit qu'un tel $x$ est racine de $\mathrm{Res}_Y(A, B)$. À l'inverse, si $x$ annule le résultant, alors d'après le Corollaire 6.3, les polynômes $A(x, Y)$ et $B(x, Y)$ ont un pgcd de degré au moins 1. Comme le corps est algébriquement clos, il existe alors $y \in K$ racine de ce pgcd tel que $A(x, y) = B(x, y) = 0$. ■

Graphiquement, les racines « dégénérées » du second cas se traduisent par la présence d'asymptotes verticales.

**Exemple 6.9** Avec $A = X^2Y + X + 1$, $B = XY - 1$, on obtient $\mathrm{Res}_Y(A, B) = -X(2X + 1)$ (asymptote en $X = 0$, « vraie » solution en $X = -\frac{1}{2}$). Les courbes correspondantes sont représentées en Figure 6.5 : celle pour A en bleu, celle pour B en rouge.

Une fois calculé, le résultant donnant les coordonnées des abscisses des intersections des deux courbes, il est souhaitable d'obtenir les ordonnées correspondantes. Dans le cas simple où l'avant-dernier reste de l'algorithme d'Euclide est de degré 1, il fournit bien une telle paramétrisation.

**Exemple 6.10** Sur les deux polynômes A et B de l'Exemple 6.8, vus comme poly-
nômes dans $Q(X)[Y]$, la suite des restes est

$$ 
\begin{aligned}
& (X^4 - 13X^2 + 15)Y^2 - 4(X^2 + 2)(X^2 + 3)Y + (X^6 - 2X^4 - 8X^2 + 10), \\
& (4X^8 - 344X^6 - 1243X^4 - 301X^2 - 21)Y + (84X^8 - 372X^6 + 169X^4 + 143X^2 + 15), \\
& 1.
\end{aligned}
 $$

Un calcul de pgcd du coefficient de Y dans l'avant-dernier reste avec $\text{Res}_Y(A, B)$ montre que ces deux polynômes sont premiers entre eux. Pour chaque racine X de $\text{Res}_Y(A, B)$ il y a donc un unique point d'intersection aux deux courbes, d'ordonnée donnée par

$$ Y = - \frac{84X^8 - 372X^6 + 169X^4 + 143X^2 + 15}{4X^8 - 344X^6 - 1243X^4 - 301X^2 - 21}. $$

En utilisant un calcul de pgcd étendu, ceci peut être récrit, comme expliqué plus haut, en un polynôme de degré au plus 13 en X.

Le même calcul sur A et B vus comme polynômes dans $Q(Y)[X]$ donne B comme dernier reste avant le pgcd. Ceci correspond aux deux points ayant la même projection sur l'axe des Y.

### Implicitation

Une autre application géométrique du résultant est l'implicitation de courbes du plan. Étant donnée une courbe paramétrée

$$ X = A(T), \quad Y = B(T), \quad A, B \in \mathbb{K}(T), $$

il s'agit de calculer un polynôme non trivial en X et Y qui s'annule sur la courbe. Il suffit pour cela de prendre le résultant en T des numérateurs de $X - A(T)$ et de $Y - B(T)$.

**Exemple 6.11** La courbe « en fleur » de la Figure 6.4, plus savamment appelée « quadrifolium », peut aussi être donnée sous la forme

$$ X = \frac{4T(1 - T^2)^2}{(1 + T^2)^3}, \quad Y = \frac{8T^2(1 - T^2)}{(1 + T^2)^3}. $$

Il suffit d'effectuer le calcul du résultant

$$ \text{Res}_T((1 + T^2)^3 X - 4T(1 - T^2)^2, (1 + T^2)^3 Y - 8T^2(1 - T^2)) $$

pour retrouver (à un facteur constant près) le polynôme A de l'Exemple 6.8.

### Propriétés

L'essentiel des propriétés du résultant est fourni par le théorème suivant.
**Théorème 6.5** — Formule de Poisson. Si les polynômes $A$ et $B$ de $\mathbb{A}[X]$ s'écrivent

$$A = a(X - \alpha_1) \cdots (X - \alpha_m), \quad B = b(X - \beta_1) \cdots (X - \beta_n),$$

alors, le résultant de $A$ et $B$ vaut

$$\begin{aligned} \operatorname{Res}(A, B) &= a^n b^m \prod_{i,j} (\alpha_i - \beta_j) = (-1)^{mn} b^m \prod_{1 \le j \le n} A(\beta_j) \\ &= a^n \prod_{1 \le i \le m} B(\alpha_i) = (-1)^{mn} \operatorname{Res}(B, A). \end{aligned}$$

**Démonstration.** Il suffit de prouver la première égalité. Les trois suivantes en sont des conséquences immédiates.

Le facteur $a^n b^m$ est une conséquence de la multilinéarité du déterminant. On peut donc se restreindre au cas où $a = b = 1$. Il est commode de considérer le cas générique où les $\alpha_i$ et $\beta_j$ sont des indéterminées et où l'anneau $\mathbb{A}$ est $\mathbb{Z}[\alpha_1, \dots, \alpha_m, \beta_1, \dots, \beta_n]$. Si le résultat, qui est une identité entre polynômes, est vrai dans cet anneau, il l'est aussi par spécialisation pour des valeurs arbitraires des $\alpha_i$ et $\beta_j$. Le Corollaire 6.3 avec $\mathbb{K} = \mathbb{Q}(\alpha_1, \dots, \alpha_m, \beta_1, \dots, \beta_n)$ montre que le polynôme obtenu en remplaçant la variable $\beta_j$ par la variable $\alpha_i$ dans le résultant est le polynôme identiquement nul. Il s'ensuit que le produit des $\alpha_i - \beta_j$ divise le résultant. Par ailleurs, le degré en $\alpha_i$ de chacune des $n$ premières lignes de la matrice de Sylvester est 1, et ce degré est nul pour les $m$ autres lignes. Ceci donne une borne $n$ pour le degré en chaque $\alpha_i$ du résultant et de la même manière une borne $m$ pour son degré en chaque $\beta_j$. Il s'ensuit que le résultant est égal au produit des $\alpha_i - \beta_j$ à un facteur constant près. Ce facteur est indépendant des $\alpha_i$ et $\beta_j$. En choisissant $\alpha_1 = \cdots = \alpha_m = 0$, c'est-à-dire $A = X^m$, la matrice de Sylvester est triangulaire, et son déterminant vaut $B(0)^m$, ce qui donne le facteur 1 et conclut la preuve. ■

**Exemple 6.12** Si $A = a \prod_{i=1}^m (X - x_i)$, la formule de Poisson entraîne que le discri-
minant $\operatorname{disc}(A)$ de $A$ est égal à :

$$(-1)^{\frac{m(m-1)}{2}} a^{m-2} \prod_{i=1}^{m} A'(x_i) = (-1)^{\frac{m(m-1)}{2}} a^{m-2} \prod_{i=1}^{m} a \prod_{j \neq i} (x_i - x_j) = a^{2m-2} \prod_{i<j} (x_i - x_j)^2.$$

**Corollaire 6.6** — Multiplicativité du résultant. Pour tous polynômes $A, B, C$
de $\mathbb{A}[X]$, $\operatorname{Res}(AB, C) = \operatorname{Res}(A, C)\operatorname{Res}(B, C)$.

**Démonstration.** L'anneau $\mathbb{A}$ étant intègre, il possède un corps de fractions qui a lui-
même une clôture algébrique dans laquelle les polynômes $A, B$ et $C$ peuvent s'écrire
sous la forme utilisée dans le théorème précédent. L'identité sur les résultats (qui
appartiennent à $\mathbb{A}$ par définition) est alors vérifiée en considérant les produits deux à
deux de racines. ■
Une dernière propriété utile du résultant est la suivante.

**Proposition 6.7** Il existe U et V dans $\mathcal{A}[X]$ avec $\deg U < \deg B$ et $\deg V < \deg A$, tels que le résultant de $A$ et de $B$ s'écrive $\mathrm{Res}(A, B) = UA + VB$.

*Démonstration.* En ajoutant à la dernière colonne de la matrice de Sylvester le produit des colonnes précédentes par des puissances adaptées de $X$, on fait apparaître dans la dernière colonne les polynômes $X^{\deg B-1}A, \dots, XA, A$ et $X^{\deg A-1}B, \dots, XB, B$, sans avoir changé le déterminant. Le développement du déterminant par rapport à la dernière colonne permet alors de conclure. ■

L'algorithme d'Euclide étendu montre que le résultant $\mathrm{Res}(A, B)$ (qui est un multiple par un élément du corps des fractions $\mathbb{K}$ de l'un des restes euclidiens « standards », peut s'écrire comme une combinaison polynomiale de $A$ et $B$, à coefficients dans $\mathbb{K}[X]$. La Proposition 6.7 montre que cette combinaison se fait « sans division », c'est-à-dire avec des coefficients dans $\mathcal{A}[X]$.

### Spécialisation du résultant

Spécialisation et résultant ne commutent pas toujours, comme le montre l'exemple suivant.

**Exemple 6.13** Le résultant de $aX^2 + bX + c$ et $dX + e$ vaut $ae^2 + cd^2 - bde$, qui vaut $d(cd - be)$ en $a = 0$, alors que le résultant de $bX + c$ et $dX + e$ vaut $be - cd$.

Le résultat suivant fournit une condition simple pour que le résultant de la spécialisation de deux polynômes soit égal à la spécialisation de leur résultant, propriété qui sera utile pour les calculs sur les systèmes polynomiaux de la Partie V.

**Proposition 6.8** Soient $\mathcal{A}$ et $\mathbb{B}$ deux anneaux commutatifs avec des unités respectives $1_{\mathcal{A}}$ et $1_{\mathbb{B}}$. Soit $\varphi$ un morphisme de $\mathcal{A}$ dans $\mathbb{B}$ qui envoie $1_{\mathcal{A}}$ sur $1_{\mathbb{B}}$, et que l'on étend naturellement à $\mathcal{A}[X]$ en l'appliquant coefficient par coefficient. Soient $\mathcal{A}$ et $\mathbb{B}$ deux polynômes de $\mathcal{A}[X]$. Si $\mathcal{A}$ est unitaire alors $\varphi(\mathrm{Res}(\mathcal{A}, \mathbb{B})) = \mathrm{Res}(\varphi(\mathcal{A}), \varphi(\mathbb{B}))$.

*Démonstration.* Le lemme ci-dessous donne une expression du résultant comme le déterminant de l'endomorphisme de multiplication par $B$ dans $\mathcal{A}[X]/(A(X))$, de matrice notée $M(\mathcal{A}, \mathbb{B})$ dans la base $1, X, \dots, X^{\deg(\mathcal{A})-1}$. Sa $i$-ième colonne est le reste $R_i$ de la division euclidienne de $X^i B$ par $\mathcal{A}$. Comme $\mathcal{A}$ est unitaire, les coefficients de $R_i$ sont des polynômes en les coefficients de $\mathcal{A}$ et $\mathbb{B}$. Par ailleurs, le reste de la division euclidienne de $\varphi(X^i B)$ par $\varphi(\mathcal{A})$ coïncide avec $\varphi(R_i)$. On en déduit que $\varphi(M(\mathcal{A}, \mathbb{B})) = M(\varphi(\mathcal{A}), \varphi(\mathbb{B}))$. La conclusion vient du fait que $\varphi(\det(M(\mathcal{A}, \mathbb{B}))) = \det(\varphi(M(\mathcal{A}, \mathbb{B})))$. ■

**Lemme 6.9** Soit $\mathcal{A}$ un anneau commutatif avec une unité notée $1_{\mathcal{A}}$, et soient $\mathcal{A}$ et $\mathbb{B}$ deux polynômes de $\mathcal{A}[X]$. Si $\mathcal{A}$ est unitaire, alors le résultant $\mathrm{Res}(\mathcal{A}, \mathbb{B})$ est égal au déterminant de l'endomorphisme de multiplication par $\mathbb{B}$ dans $\mathcal{A}[X]/(\mathcal{A}(X))$.
*Démonstration.* Notons $A = 1_A X^m + a_{m-1} X^{m-1} + \cdots + a_0$, et $B = b_n X^n + \cdots + b_0$, avec $b_n \neq 0$. La matrice de Sylvester $Syl(A, B)$ s'écrit :

$$ \mathrm{Syl}(A, B) = \begin{pmatrix} 1_A & \cdots & a_1 & a_0 & & \\ & \ddots & & \ddots & & \\ & & 1_A & a_{m-1} & \cdots & a_0 \\ b_n & \cdots & b_1 & b_0 & & \\ & \ddots & & \ddots & & \\ & & b_n & b_{n-1} & \cdots & b_0 \end{pmatrix} $$

En notant $R_i = r_{i,m-1}X^{m-1} + \cdots + r_{i,0}$ le reste de la division euclidienne de $X^i B$ par $A$, et en effectuant les opérations sur les lignes de $Syl(A, B)$ correspondant à ces divisions, on constate que le déterminant de cette matrice est égal au déterminant de la matrice

$$ \begin{pmatrix} r_{m-1,m-1} & r_{m-1,m-2} & \cdots & r_{m-1,0} \\ r_{m-2,m-1} & r_{m-2,m-2} & \cdots & r_{m-2,0} \\ \vdots & \vdots & \ddots & \vdots \\ r_{0,m-1} & r_{0,m-2} & \cdots & r_{0,0} \end{pmatrix}, $$

dont la transposée n'est autre que la matrice de multiplication par $B$ dans $A[X]/(A(X))$, pour la base $X^{m-1}, \dots, X, 1$.

## Calcul en coût quadratique

En utilisant la Définition 6.1 et les résultats du Chapitre 8, le résultant de deux polynômes $A, B$ de $K[X]$ peut se calculer naïvement en $O(MM(n)) = O(n^\theta)$ opérations dans $K$, où $n = \max(\deg(A), \deg(B))$ et où $2 \le \theta \le 3$ est un exposant réalisable pour la multiplication matricielle. Un algorithme plus efficace, de complexité seulement quadratique, s'obtient comme conséquence du Théorème 6.5 et est laissé en exercice.

**Exercice 6.5 — Algorithme de type Euclide pour le résultant.** Soient $A$ et $B$ deux polynômes de $K[X]$ de degrés $m = \deg(A)$ et $n = \deg(B)$.

1. Soit $A = QB + R$ la division euclidienne dans $K[X]$ de $A$ par $B$. Si $r = \deg(R)$, et si $b_n$ est le coefficient dominant de $B$, montrer que

$$ \mathrm{Res}(A, B) = (-1)^{mn} \cdot b_n^{m-r} \cdot \mathrm{Res}(B, R). $$

2. Écrire un algorithme de type Euclide (c'est-à-dire à base de divisions répétées), qui calcule le résultant $\mathrm{Res}(A, B)$ en $O(mn)$ opérations dans $K$.

## Calcul avec des nombres algébriques

L'idée que les polynômes sont de bonnes structures de données pour représenter leur racines amène à chercher des algorithmes pour effectuer les opérations de base sur ces racines, comme la somme ou le produit. Le résultant répond à cette attente.
**Proposition 6.10** Soient $A = \prod_i (X - \alpha_i)$, $B = \prod_j (X - \beta_j)$ et $C$ des polynômes unitaires de $K[X]$, avec $C$ et $A$ premiers entre eux. Alors

$$ \mathrm{Res}_{X}(A(X), B(T-X)) = \prod_{i,j} (T - (\alpha_i + \beta_j)), $$

$$ \mathrm{Res}_{X}(A(X), B(T+X)) = \prod_{i,j} (T - (\beta_j - \alpha_i)), $$

$$ \mathrm{Res}_{X}(A(X), X^{\deg B}B(T/X)) = \prod_{i,j} (T - \alpha_i \beta_j), $$

$$ \mathrm{Res}_{X}(A(X), C(X)T - B(X)) = \mathrm{Res}(A, C) \prod_{i} \left( T - \frac{B(\alpha_i)}{C(\alpha_i)} \right). $$

*Démonstration.* C'est une application directe du Théorème 6.5.

Cette proposition fournit une preuve effective du résultat suivant.

**Corollaire 6.11** L'ensemble des éléments algébriques sur $K$ est un corps.

**Exemple 6.14** On sait que $\sqrt{2}$ est racine de $X^2 - 2$, tout comme $\sqrt{3}$ est racine de $X^2 - 3$. Un polynôme ayant pour racine $\sqrt{2} + \sqrt{3}$ est donné par

$$ \mathrm{Res}_{X}(X^2 - 2, (T - X)^2 - 3) = T^4 - 10T^2 + 1. $$

Ces opérations ne distinguent pas les racines de $A$ et $B$, les quatre racines du résultant sont donc $\pm\sqrt{2} \pm\sqrt{3}$.

**Exemple 6.15** Ramanujan a prouvé l'identité suivante

$$ \sqrt[3]{\cos \frac{2\pi}{7}} + \sqrt[3]{\cos \frac{4\pi}{7}} + \sqrt[3]{\cos \frac{8\pi}{7}} = \sqrt[3]{\frac{5-3\sqrt{7}}{2}}. $$

Pour prouver automatiquement cette identité, une méthode consiste à construire un polynôme annulant son membre gauche, puis à isoler la racine correspondant numériquement à sa valeur, et à vérifier qu'il s'agit bien du membre droit.

En posant $x = \exp(2i\pi/7)$, les trois racines cubiques sont annulées par les polynômes

$$ P_1 = X^3 - \frac{x+x^{-1}}{2}, \quad P_2 = X^3 - \frac{x^2+x^{-2}}{2}, \quad P_3 = X^3 - \frac{x^4+x^{-4}}{2}. $$

Un premier polynôme $P_4(Y) = \mathrm{Res}_X(P_1(Y-X), P_2(X))$ de degré 9 est d'abord obtenu, puis $P_5 = \mathrm{Res}_X(P_4(Y-X), P_3(X))$ de degré 27, et enfin le résultant de $x^7 - 1$ et du
numérateur de P₅ vu maintenant comme polynôme en x fournit

$$32768Y^6(Y - 3)(Y^2 + 3)^2(Y^2 + 3Y + 9)(Y^2 + 3Y + 3)^3 \\ \times (Y^2 - 3Y + 3)^3(4Y^9 - 30Y^6 + 75Y^3 + 32)^6 \\ \times (64Y^{18} + 768Y^{15} + 7728Y^{12} + 17152Y^9 + 16008Y^6 + 3864Y^3 + 343)^6.$$

Une évaluation numérique montre que c'est le facteur

$$4Y^9 - 30Y^6 + 75Y^3 + 32$$

qui annule le membre gauche (il faut faire un peu attention : Ramanujan utilise la convention que la racine cubique d'un réel négatif est un réel négatif). On peut alors résoudre ce polynôme par les formules de Cardan pour découvrir le membre droit : la seule racine réelle est

$$-\sqrt[3]{\frac{3\sqrt{7}-5}{2}}.$$

On a donc non seulement prouvé l'identité, mais découvert automatiquement le membre droit.

**Exercice 6.6** Montrer à l'aide de résultats les identités suivantes :

$$\frac{\sin \frac{2\pi}{7}}{\sin^2 \frac{3\pi}{7}} - \frac{\sin \frac{\pi}{7}}{\sin^2 \frac{2\pi}{7}} + \frac{\sin \frac{3\pi}{7}}{\sin^2 \frac{\pi}{7}} = 2\sqrt{7} \quad \text{et} \quad \frac{\sin^2 \frac{3\pi}{7}}{\sin \frac{2\pi}{7}} - \frac{\sin^2 \frac{2\pi}{7}}{\sin \frac{\pi}{7}} + \frac{\sin^2 \frac{\pi}{7}}{\sin \frac{3\pi}{7}} = 0.$$

<figure><img src="image_24.png" /></figure>

## Calcul du résultant en deux variables

On dispose à l'heure actuelle d'algorithmes rapides (de complexité quasi-optimale) pour le calcul du résultant avec une seule variable. Le meilleur algorithme connu actuellement pour le calcul du résultant en deux variables est une méthode d'évaluation-interpolation qui se ramène à une seule variable (voir l'Exercice 6.10). Cet algorithme n'est pas quasi-optimal. En revanche, pour les trois premières opérations de la Proposition 6.10, des algorithmes quasi-optimaux à base de multiplication rapide de séries existent, ils ont été présentés au Chapitre 3 (Section 3.3, page 68).

## Sous-résultants

Le nombre d'opérations dans le corps des coefficients n'est pas une mesure suffisante de la complexité des calculs lorsque les opérations elles-mêmes ont une complexité variable. C'est le cas pour le calcul de pgcd dans $Q[X]$ et dans $K(Y)[X]$. Dans ces corps de coefficients, on constate empiriquement les phénomènes suivants :

- l'algorithme d'Euclide amène à faire des divisions, et introduit des dénomineurs au cours du calcul;

- la taille des coefficients croît rapidement;

- ces coefficients peuvent souvent se « simplifier ».
**Exemple 6.16** L'exécution de l'algorithme d'Euclide sur les polynômes

$$A = 115X^5 + 7X^4 + 117X^3 + 30X^2 + 87X + 44, \quad B = 91X^4 + 155X^3 + 3X^2 + 143X + 115$$

produit les restes successifs suivants :

$$\frac{3601622}{8281}X^3 - \frac{1196501}{8281}X^2 + \frac{151912}{637}X + \frac{2340984}{8281}, \\ \frac{189886027626841}{12971681030884}X^2 - \frac{57448278681703}{3242920257721}X - \frac{17501090665331}{3242920257721}, \\ \frac{3748556212578804983806085060}{4354148470945709877351001}X + \frac{141833360915123969328014892}{334934497765054605950077}$$

En multipliant ces restes par des constantes bien choisies, on obtient les restes :

$$3601622X^3 - 1196501X^2 + 1974856X + 2340984,$$

$$22930325761X^2 - 27749440252X - 8453612204,$$

$$288979986761465X + 142143002707719.$$

Il est souhaitable d'éviter le calcul sur les rationnels (ou les fractions rationnelles) : leur addition requiert des calculs de multiplications et éventuellement de pgcd. Une première idée est d'éviter totalement les divisions en utilisant des pseudo-restes.

**Définition 6.2** Si A et B sont des polynômes de $A[X]$ et si $b$ est le coefficient de tête de B, le pseudo-reste $\bar{R}$ de A et B est défini par

$$b^{\deg A - \deg B + 1} A = \bar{Q}B + \bar{R},$$

où $\bar{Q}, \bar{R} \in A[X]$ avec $\deg \bar{R} < \deg \bar{Q}$.

Remplacer les calculs de restes de l'algorithme d'Euclide par des calculs de pseudo-restes évite d'introduire des dénominateurs. Cette idée seule n'est pas suffisante : les coefficients de ces pseudo-restes croissent trop.

**Exemple 6.17** Sur le même exemple, la modification de l'algorithme d'Euclide produit la suite de pseudo-restes

$$\begin{aligned} & 3601622X^3 - 1196501X^2 + 1974856X + 2340984, \\ & 189886027626841X^2 - 229793114726812X - 70004362661324, \\ & 257057096083899261191107914552182660X \\ & \qquad + 126440823512296156588839626542149756. \end{aligned}$$

Une possibilité pour éviter cette croissance est de diviser ces polynômes par le pgcd de leurs coefficients à chaque étape. C'est ce que nous avons fait dans l'Exemple 6.16 ci-dessus. On parle alors de suite de pseudo-restes primitifs. Cependant, ces calculs de
**Entrée** A et B dans $A[X]$ avec $\deg B < \deg A$.

**Sortie** Le dernier sous-résultant non nul.

1. Initialiser $f = g = s = 1$.

2. Tant que $\deg B > 0$ faire :

a. $d = \deg A - \deg B$.

b. Calculer le pseudo reste R de A par B.

c. Si les degrés de A et B sont impairs, faire $s = -s$.

d. Faire $A := B$ et mettre dans B le quotient de la division de R par $fg^d$.

e. Copier le coefficient dominant de A dans $f$.

f. Mettre à jour $g$ en lui affectant le quotient de la division de $f^d$ par $g^{d-1}$.

3. $d = \deg A$;

4. Renvoyer le quotient de la division de $sB^d$ par $g^{d-1}$.

**Algorithme 6.6 – Algorithme des sous-résultants.**

pgcd de coefficients sont trop coûteux.

Le calcul des sous-résultants, donné dans l'Algorithme 6.6, est une modification de l'algorithmme d'Euclide qui évite les divisions, tout en prévoyant des facteurs communs qui apparaissent dans les coefficients de la suite des pseudo-restes de sorte à limiter leur croissance. Ces pseudo-restes sont appelés des sous-résultants. Dans cet algorithme, on se contente de renvoyer le dernier sous-résultant non nul; on conçoit aisément comment modifier l'algorithmme pour obtenir *n'importe* quel sous-résultant (spécifié par des conditions de degré, par exemple).

**Exemple 6.18** Toujours sur les mêmes polynômes, la suite des sous-résultants redonne les polynômes que nous avions déjà calculés par simplification :

$$ 
\begin{aligned}
& 3601622X^3 - 1196501X^2 + 1974856X + 2340984, \\
& 22930325761X^2 - 27749440252X - 8453612204, \\
& 288979986761465X + 142143002707719.
\end{aligned}
$$

La suite des sous-résultants de cet exemple est donc primitive : les facteurs prédits par l'algorithmme du sous-résultant ont suffi à éliminer tout facteur commun entre les coefficients des pseudo-restes.

Comme le résultant, les sous-résultants sont liés à la matrice de Sylvester et à l'algorithmme d'Euclide. Une formule de Cramer sur une sous-matrice de la matrice de Sylvester donne le résultat suivant.

**Proposition 6.12** Toutes les divisions effectuées au cours de l'algorithmme des sous-résultants sont exactes.

La preuve de ce résultat est technique et omise ici.
<figure><img src="image_25.png" /><figcaption>FIGURE 6.7 – Paramétrisation des solutions d'un système algébrique par les racines d'un résultant.</figcaption></figure>

En corollaire, on remarque en outre qu'il est assez facile de déduire de ce résultat des estimations sur la « taille » des coefficients qui apparaissent au cours de l'algorithme. Un cas particulier intéressant est celui où A est l'anneau de polynômes $K[Y]$. Alors, si A et B ont des degrés totaux *m* et *n*, tous les sous-résultants ont des degrés bornés par *mn*. Ces bornes permettent un calcul du résultant par évaluation-interpolation dès que l'on dispose d'un algorithme efficace à une variable.

### Calcul de la paramétrisation

Les sous-résultants permettent de finir la « résolution » (commencée en page 118) des systèmes de deux polynômes à deux variables. Pour simplifier, nous faisons l'hypothèse que pour tout *x* racine du résultant $R = \text{Res}_Y(A, B)$, il n'y a qu'un seul *y* tel que $A(x, y) = B(x, y) = 0$. Dans ce cas, on peut montrer que le sous-résultant $S_1$ (de degré 1 en Y) est non nul ; écrivons-le sous la forme $S_1 = P_0(X)Y + Q_0(X)$. Comme pour le résultant, il existe des polynômes U et V dans $K[X, Y]$ tels qu'on ait l'égalité

$$AU + BV = P_0(X)Y + Q_0(X).$$

On en déduit donc que toutes les solutions $(x, y)$ du système $A = B = 0$ satisfont l'équation $P_0(x)y = Q_0$. Autrement dit, en écartant les solutions dégénérées où $P_0(x) = 0$, on obtient l'ordonnée des points solutions en évaluant la fraction rationnelle $-Q_0/P_0$ sur les racines du résultant R. Autrement dit, cette procédure permet de décrire les solutions du système sous la forme

$$\begin{cases} y = Q(x) \\ R(x) = 0. \end{cases}$$

Géométriquement, ce polynôme Q permet de retrouver $y$ à partir de $x$, il permet donc d'effectuer l'opération symbolisée sur la Figure 6.7, la paramétrisation des solutions du système par les racines de R.
**Exemple 6.19** Sur les mêmes polynômes qu’à l’Exemple 6.8, le sous-résultant redonne la paramétrisation calculée à partir de l’algorithme d’Euclide, mais en étant plus économe dans les opérations sur les coefficients.

## 6.3 Algorithme d’Euclide rapide

Le point central de cette section est la présentation d'un algorithme rapide pour le pgcd, qui est un algorithme récursif à base de multiplications et de divisions euclidiennes rapides. Historiquement, c'est d'abord un algorithme rapide pour le pgcd d'entiers (le « demi-pgcd ») qui est apparu, il a ensuite été étendu aux polynômes. Nous décrivons maintenant l'algorithme dit du « demi-pgcd ». Pour simplifier, nous nous plaçons dans le cas du pgcd de polynômes, mais cette approche s'étend au cas des entiers.

Soient donc A et B dans $K[X]$, avec $n = \text{deg} A$ et $\text{deg} B < n$ (ce qui n'est pas une forte restriction). Posons comme précédemment $R_0 = A$, $R_1 = B$. Soient ensuite $R_i$ les restes successifs de l'algorithme d'Euclide et soit en particulier $R_N = \text{pgcd}(A, B)$ le dernier non nul d'entre eux. D'après les égalités en page 115, il existe une matrice $2 \times 2$, notée $M_{A,B}$, de cofacteurs telle que :

$$ M_{A,B} \begin{pmatrix} R_0 \\ R_1 \end{pmatrix} = \begin{pmatrix} U_N & V_N \\ U_{N+1} & V_{N+1} \end{pmatrix} \begin{pmatrix} R_0 \\ R_1 \end{pmatrix} = \begin{pmatrix} R_N \\ 0 \end{pmatrix}. $$

L'algorithme de pgcd rapide calcule d'abord cette matrice; à partir de là, on en déduit aisément le pgcd, pour $O(M(n))$ opérations supplémentaires.

Notons qu'une étape de l'algorithme d'Euclide classique peut elle aussi s'écrire sous une forme matricielle. En effet, si $Q_i$ est le quotient de la division de $R_{i-1}$ par $R_i$, on peut écrire :

$$ \begin{pmatrix} 0 & 1 \\ 1 & -Q_i \end{pmatrix} \begin{pmatrix} R_{i-1} \\ R_i \end{pmatrix} = \begin{pmatrix} R_i \\ R_{i+1} \end{pmatrix}. $$

Ainsi, la matrice $M_{A,B}$ est un produit de matrices élémentaires de ce type. Elle est inversible.

### Demi-pgcd : définition

Comme étape intermédiaire pour obtenir une division euclidienne rapide, on utilise l'algorithme dit du « demi-pgcd » (half-gcd en anglais), qui permet de faire des « pas de géant » dans la liste des restes successifs.

Le demi-pgcd est défini comme suit. Les degrés des restes successifs $R_i$ décroissent strictement, donc il existe un unique indice $j$ tel que :

$$ - \text{deg} R_j \geq \frac{n}{2}; $$

$$ - \text{deg} R_{j+1} < \frac{n}{2}. $$

On a vu en page 115 qu'il existe $U_j, V_j, U_{j+1}, V_{j+1}$ tels que :

$$ U_j R_0 + V_j R_1 = R_j \quad \text{et} \quad U_{j+1} R_0 + V_{j+1} R_1 = R_{j+1}, $$
**Entrée** R₀ = A et R₁ = B dans $\mathbb{K}[X]$ avec deg B < deg A = n.

**Sortie** La matrice M<sub>A,B</sub>.

1. Soit M<sub>dpgcd</sub> = dpgcd(A, B).

2. Calculer R<sub>j</sub> et R<sub>j+1</sub>.

3. Si R<sub>j+1</sub> = 0, renvoyer M<sub>dpgcd</sub>.

4. Soit M la matrice de la division euclidienne de R<sub>j</sub> par R<sub>j+1</sub>, et R<sub>j+2</sub> le reste correspondant.

5. Si R<sub>j+2</sub> = 0, renvoyer MM<sub>dpgcd</sub>.

6. Calculer M<sub>R<sub>j+1</sub>,R<sub>j+2</sub></sub> par un appel récursif.

7. Renvoyer M<sub>R<sub>j+1</sub>,R<sub>j+2</sub></sub>MM<sub>dpgcd</sub>.

Algorithme 6.8 – Algorithme d'Euclide rapide via le demi-pgcd rapide.

ces polynômes étant en outre uniques si on rajoute les conditions de degré voulues
sur la relation de Bézout. Ceci peut se récrire sous la forme matricielle suivante :

$$
\begin{pmatrix} U_j & V_j \\ U_{j+1} & V_{j+1} \end{pmatrix} \begin{pmatrix} R_0 \\ R_1 \end{pmatrix} = \begin{pmatrix} R_j \\ R_{j+1} \end{pmatrix}.
$$

Pour fixer les idées, en première approximation, on peut estimer que les polynômes
U<sub>j</sub>, V<sub>j</sub>, U<sub>j+1</sub>, V<sub>j+1</sub> ont des degrés de l'ordre de n/2 (attention, ce n'est qu'une estimation,
pas une égalité).

Notons M<sub>dpgcd</sub> la matrice ci-dessus donnant les restes R<sub>j</sub> et R<sub>j+1</sub>. L'algorithme du demi-pgcd (noté dpgcd ci-dessous) a pour vocation de calculer cette matrice. Avant d'en étudier le fonctionnement, voyons comment il permet d'obtenir le calcul du pgcd étendu. L'idée est qu'un appel à dpgcd en degré n permet d'obtenir les restes de degré approximativement n/2 ; alors, un appel en degré n/2 permet d'obtenir les restes de degré approximativement n/4, et ainsi de suite jusqu'à trouver le pgcd. Cette approche est détaillée dans l'Algorithme 6.8. À l'étape (4), R<sub>j</sub> et R<sub>j+1</sub> ont par définition des degrés qui encadrent n/2, mais on n'a pas de borne supérieure sur le degré de R<sub>j</sub>, qui peut être proche de n. Aussi, pour obtenir deux polynômes de degré au plus n/2 pour l'appel récursif, il est nécessaire d'effectuer une étape de division euclidienne.

Soit H(n) la complexité de l'algorithmme de demi-pgcd sur des entrées de degré au plus n. Pour estimer la complexité de l'algorithmme de pgcd rapide, on fait l'hypothèse que H(n + n') ≥ H(n) + H(n') et que M(n) est négligeable devant H(n). Ces hypothèses sont vérifiées pour l'algorithmme proposé en Section 6.3.

**Proposition 6.13** *L'Algorithme 6.8 calcule M<sub>A,B</sub> en O(H(n)) opérations.*

*Démonstration.* La validité de l’algorithme est immédiate, puisque toutes les matrices calculées ne sont au fond que des matrices qui font la transition de deux restes successifs à deux autres restes successifs. Multiplier ces matrices permet de composer ces opérations de transition.

Estimons la complexité. Le coût du dpgcd est H(n). Ensuite, toutes les opérations des étapes (2) à (5) ont une complexité en O(M(n)), en utilisant pour l'étape (4) une
**Entrée** A et B dans $K[X]$, avec $\deg B < \deg A = n$.

**Sortie** La matrice $M_{\text{dpgcd}}$ de A et B.

1. Soit $m = \lceil \frac{n}{2} \rceil$. Si $\deg B < m$, renvoyer la matrice identité.

2. Soient $f$ et $g$ les quotients respectifs de A et B par $X^m$ (qui ont des degrés approximativement $n/2$).

3. Soient $M = \text{dpgcd}(f, g)$, et $A'$, $B'$ définis par

$$M\begin{pmatrix} A \\ B \end{pmatrix} = \begin{pmatrix} A' \\ B' \end{pmatrix}.$$

4. Si $\deg B' < m$, renvoyer M.

5. Soient $C' = A'$ mod $B'$, et $M'$ la matrice de la division euclidienne correspondante.

Informellement, $B'$ a un degré au pire de l'ordre de $\frac{3n}{4}$. Par contre, $A'$ peut avoir un degré plus grand. L'étape de division euclidienne permet de se ramener à des polynômes de degrés de l'ordre de $\frac{3n}{4}$.

6. Soient $\ell = 2m - \deg B'$, $b$ et $c$ les quotients respectifs de $B'$ et $C'$ par $X^\ell$.

$\ell$ est de l'ordre de $\frac{n}{4}$, de sorte que $b$ et $c$ ont des degrés de l'ordre de $\frac{n}{2}$.

7. Soit $M'' = \text{dpgcd}(b, c)$. Renvoyer la matrice $M''M'M$.

**Algorithme 6.9 – Algorithme du demi-pgcd (dpgcd) rapide.**

division euclidienne rapide (cf. Théorème 4.1 en page 83). On effectue ensuite un appel récursif, puis de nouveau des opérations dont le coût est en $O(M(n))$. Ainsi, la complexité $B(n)$ de l'algorithme de pgcd rapide satisfait à la récurrence :

$$B(n) \leq B(n/2) + H(n) + CM(n),$$

C étant une constante. Le théorème « diviser pour régner » permet de conclure. ■

Autrement dit, le coût du pgcd est, à une constante près, le même que celui du demi-pgcd. Il devient donc légitime de se consacrer uniquement à ce dernier.

### Demi-pgcd : algorithme

Pour mettre en place une stratégie « diviser pour régner » pour le demi-pgcd, on utilise un moyen de « couper les polynômes » en deux. L'idée est de ne garder que les coefficients de poids forts des polynômes d'entrée et d'utiliser le fait que le quotient de la division euclidienne de deux polynômes ne dépend que de leurs coefficients de poids forts, d'une façon tout à fait quantifiée par la suite.

Remarquons qu'on accède aux coefficients de poids fort d'un polynôme en prenant son quotient par un polynôme de la forme $X^k$: par exemple, si A est

$$X^{10} + 2X^9 + 5X^8 + 3X^7 + 11X^6 + \dots,$$
le quotient de A par X⁶ est

$$X^4 + 2X^3 + 5X^2 + 3X + 11.$$

Voyons comment cette idée permet d'obtenir notre algorithme. Étant donnés deux polynômes A et B de degrés de l'ordre de $n$, on leur associe $f$ et $g$ de degré approximativement $\frac{n}{2}$, en nommant $f$ et $g$ les quotients respectifs de A et B par $X^{\frac{n}{2}}$: on a jeté les coefficients de petits degrés.

On calcule la matrice M du demi-pgcd de $f$ et $g$ (grosso modo, les éléments de cette matrice ont degré $\frac{n}{4}$). La remarque ci-dessus permet de montrer qu'on obtient ainsi une matrice de transition vers les restes de A et B eux-mêmes.

On applique donc cette matrice aux polynômes initiaux; on obtient des restes B' et C' dont les degrés sont de l'ordre de $\frac{3n}{4}$. On effectue alors une seconde fois le même type d'opération, avec un appel récursif en degré $\frac{n}{2}$, pour gagner à nouveau $\frac{n}{4}$, et finalement arriver en degré $\frac{n}{2}$.

Les choix exacts des degrés de troncature sont subtils, et on admettra que les valeurs données dans l'Algorithme 6.9 sont correctes.

**Proposition 6.14** L'Algorithme 6.9 calcule la matrice du demi-pgcd de A et B en $O(M(n) \log n)$ opérations dans $\mathbb{K}$.

*Démonstration.* Comme indiqué plus haut, on admet que les choix des degrés de troncature permettent d'assurer la validité de l'algorithme. Il est plus facile d'en effectuer l'analyse de complexité. Le coût $H(n)$ peut être majoré par deux fois le coût en degré au plus $\frac{n}{2}$ (les deux appels ont lieu lors des étapes (3) et (7)), plus un certain nombre de multiplications de polynômes et une division euclidienne. En remarquant que tous les produits se font en degré au plus $n$, on en déduit la récurrence

$$H(n) \leq 2H(n/2) + CM(n),$$

où C est une constante. La conclusion découle comme d'habitude du lemme « diviser pour régner ».

**Proposition 6.15** Le calcul du résultant de deux polynômes A et B à coefficients dans $\mathbb{K}$ peut se faire en $O(M(n) \log n)$ opérations dans $\mathbb{K}$.

*Démonstration.* Il s'agit d'adapter l'Algorithme 6.9 pour accélérer l'algorithme de coût quadatrique vu dans l'Exercice 6.5. Nous laissons ces détails techniques en exercice.

### Complément : calcul d'un reste choisi

On peut tirer un raffinement fort utile de l'algorithme de demi-pgcd : le calcul d'un reste particulier (sélectionné par une condition de degré), ainsi que des cofacteurs associés.
**Théorème 6.16** Soient $R_0 = A$ et $R_1 = B$, avec $\deg A = n$ et $\deg B < \deg A$, et soient $R_i$ les restes successifs de l'algorithme d'Euclide appliqué à A et B.

Soit $\ell < n$, et $R_j$ le premier reste de degré inférieur à $\ell$. On peut calculer $R_j$ et les cofacteurs associés $U_j$ et $V_j$ en $O(M(n) \log n)$ opérations de $K$.

*Démonstration succinte.* Si $\ell$ est plus petit que $\frac{n}{2}$, on fait un calcul de dpgcd pour se ramener en degré $\frac{n}{2}$, et on effectue un appel récursif. Si $\ell$ est plus grand que $\frac{n}{2}$, on fait un calcul de dpgcd sur les polynômes divisés par $X^{n-2\ell}$. ■

## Le cas des entiers

Les mêmes idées algorithmiques s'étendent au calcul sur les entiers, mais il est beaucoup plus difficile d'écrire un algorithme correct dans ce cas. On se contentera d'énoncer les résultats de complexité suivants :

**Théorème 6.17** Soient $R_0 = A \in \mathbb{N}$, $R_1 = B \in \mathbb{N}$, avec $B \le A$ et $R_i$ les restes de l'algorithme d'Euclide appliqué à A et B. On peut calculer :

- le pgcd de A et B,

- le premier reste inférieur à un entier $\ell < A$,
ainsi que les cofacteurs associés en $O(M_Z(\log A) \log \log A)$ opérations binaires.

## Exercices

**Exercice 6.7** — Factorisation sans carré. Soit $K$ un corps de caractéristique nulle. Si $F \in K[X]$ se factorise en irréductibles comme $\prod_i f_i^{\alpha_i}$, le polynôme $\bar{F} = \prod_i f_i$ est appelé la partie sans carré de F. Si $n = \deg(F)$, montrer qu'on peut calculer les coefficients de la partie sans carré de F à partir des coefficients de F en $O(M(n) \log n)$ opérations dans $K$. ■

**Exercice 6.8** Soient $f, g \in K[X]$ des polynômes unitaires.

1. Soit $N \in \mathbb{N} \setminus \{0\}$. Montrer que l'unique polynôme unitaire de $K[X]$ dont les racines sont les puissances N-ièmes des racines de $f$ peut être obtenu par un calcul de résultant.

2. Si $f$ est le polynôme minimal d'un nombre algébrique $\alpha$, montrer qu'on peut déterminer un polynôme annulateur de $g(\alpha)$ à l'aide d'un résultant.

3. Calculer le polynôme minimal sur $\mathbb{Q}$ de $\alpha = \sqrt[3]{4} + \sqrt[3]{2} + 1$. ■

**Exercice 6.9** Soit $f \in K[X]$ un polynôme unitaire de degré $d \ge 1$. Pour $N \ge 1$, on note $G_N(f)$ l'unique polynôme unitaire de degré $d$ dont les racines sont les puissances N-ièmes des racines de $f$.

1. Exprimer $G_N(f)$ à l'aide d'un résultant de polynômes à deux variables.

2. Justifier l'appartenance à $K[X]$ du polynôme $G_N(f)$.

3. Utiliser la question (1) pour donner un algorithme pour le calcul de $G_N(f)$;
estimer sa complexité.

4. Montrer que $G_2(f)$ peut se calculer en $O(M(d))$ opérations dans $K$.

5. Si $N$ est une puissance entière de 2, montrer comment on peut calculer $G_N(f)$ en $O(M(d) \log N)$ opérations dans $K$.

**Exercice 6.10 — Résultants en deux variables.** Soient $F, G \in K[X, Y]$ ayant des degrés en $X$ et en $Y$ bornés par un entier $D \ge 1$. Montrer que le polynôme $R(Y) := \text{Res}_X(F, G)$ est de degré au plus $2D^2$. Donner un algorithme de type évaluation-interpolation pour le calcul des coefficients de $R$ et estimer son coût en opérations dans $K$.

**Exercice 6.11 — Résultant de Rothstein-Trager.** Soit $F = A/B$ une fraction rationnelle dans $Q(X)$, où $A$ et $B$ sont deux polynômes de $Q[X]$, premiers entre eux et tels que $\text{deg}(A) < \text{deg}(B)$. On suppose que tous les pôles $p_1, \dots, p_\ell$ de $F$ sont simples. La décomposition en éléments simples de $F$ s'écrit $F = \sum_{i=1}^\ell r_i / (X - p_i)$, où $r_i$ est le résidu de $F$ au pôle $p_i$.

1. Montrer que pour tout $i$, le résidu $r_i$ vaut $A(p_i)/B'(p_i)$.

2. Expliciter, à l'aide d'un résultant, un polynôme $R$ qui annule tous les résidus $r_1, \dots, r_\ell$.

3. Donner un algorithme permettant le calcul des coefficients de $R$ à partir de ceux de $A$ et de $B$, et estimer sa complexité arithmétique en fonction du degré de $B$.

**Notes**

Pour les déterminants sur des anneaux, nous renvoyons au livre de Lang [Lan02]. L'introduction du résultant et aussi des sous-résultants remonte à Bézout et Sylvester [Béz64; Syl39; Syl40]. Ces notions ont été introduites en calcul formel dans les années 1965 par Collins [Col66; Col67; Col71], Brown et Traub [Bro71; BT71]. L'article de synthèse de von zur Gathen et Lücking retrace leur histoire [GL03]. De nombreuses propriétés des résultats et de l'algorithm e d'Euclide sont établies de manière élémentaire à l'aide des fonctions symétriques par Lascoux [Las03]. L'article d'El Kahoui est aussi une référence très abordable pour l'algorithme des sous-résultants [El 03]. On peut aussi consulter les livres suivants pour une approche un peu plus complète : [GCL92; GG03; Knu97; Lan02].

L'idée de base l'algorithme de demi-pgcd (le quotient de la division euclidienne de deux entiers ne dépend que de leurs chiffres de poids fort) est due à Lehmer [Leh38]. Knuth [Knu71] en a fait le premier algorithme sous-quadratique, qui calcule le gcd de deux entiers de taille binaire $n$ en complexité binaire $O(M_Z(n) \log^4 n)$. Schönhage [Sch71] a montré comment abaisser cette complexité à $O(M_Z(n) \log n)$. Son algorithme a été adapté au cas polynomial par Moenck [Moe73], mais seulement pour le cas où la suite des restes est normale (de degré décroissant régulièrement de 1). Brent, Gustavson et Yun [BGY80; GY79] ont donné un algorithme de type LKS (Lehmer-Knuth-Schönhage) pour le cas polynomial général, de complexité arithmétique $O(M(n) \log n)$ en degré $n$. Cet algorithme utilise $O(n \log n)$ opérations non scalaires. Strassen [Str81; Str83] a montré que $\Omega(n \log n)$ opérations non scalaires sont
aussi nécessaires, du moins pour une entrée générique, donc l’algorithme LKS est
optimal de ce point de vue.

La proposition 6.15 est due à Schwartz [Sch80b]. Schönhage a donné un algorithme
(probabiliste) pour le calcul rapide du pgcd de deux polynômes de $\mathbb{Z}[X]$. Pour deux
polynômes de degré borné par $n$ et coefficients d’au plus $\ell$ chiffres, l’algorithme a
une complexité binaire $\tilde{O}(n(n+\ell))$, il est donc quasi-optimal (en la taille binaire de
l’entrée) lorsque $\ell > n$ [Sch88b].

L'Exercice 6.7 est inspiré par un résultat de Yun [Yun76], qui montre qu'il est
possible de déterminer toute la factorisation sans carré en la même complexité. En fait,
les problèmes du calcul du pgcd et de la factorisation sans carré sont équivalents du
point de vue de la complexité [Yun77]. Ces résultats sont abordés dans le Chapitre 18.

Il est très difficile de trouver une référence complète, lisible et sans erreur sur les
algorithmes de pgcd rapide. Le Chapitre 11 de la troisième édition de l'ouvrage de
von zur Gathen et Gerhard [GG13] fournit une bonne partie des résultats techniques
implicitement utilisés dans la Section 6.3. Les notes du livre de Yap [Yap00] nous ont
beaucoup inspirés, mais elles contiennent des erreurs pour le pgcd entier. On pourra
également consulter les articles plus récents suivants : [Möl08; PW02; SZ04].

Ces algorithmes de pgcd rapide sont aussi assez délicats à implanter de manière ef-
ficace. Par exemple, pour le pgcd de polynômes, une bonne implantation doit prendre
en compte les algorithmes de multiplication utilisés. Si on utilise la FFT, il est possible
d'économiser des transformées de Fourier directes et inverses; il faut alors régler le
nombre de points d'évaluation en fonction des degrés du résultat final, et pas des
résultats intermédiaires, etc. Si l'on travaille sur un corps fini « raisonnable » (les
coefficients faisant de quelques bits jusqu'à quelques dizaines voire centaines de bits),
on peut estimer que le seuil au-delà duquel utiliser l'algorithme rapide se situe autour
des degrés 200 ou 300. Il faut noter que très peu de systèmes disposent de telles
implantations (c'est le cas pour MAGMA, MATHEMAGIX, et la bibliothèque NTL). En ce
qui concerne le pgcd des entiers, la situation est sensiblement plus délicate, toujours
à cause du problème des retenues (une bonne partie des algorithmes présentés dans
les ouvrages de référence sont incorrects, à cause de ce problème). Pour donner une
idée, dans les systèmes MAGMA, MATHEMATICA, ou des bibliothèques telles que GMP
(code toujours en développement et utilisé entre autres par les entiers de MAPLE,
MATHEMATICA et SAGE), le seuil se situe autour de nombres de 30 000 bits (10 000
chiffres décimaux). Les notions de ce chapitre — pgcd, pgcd étendu, demi-pgcd, résul-
tant, sous-résultants — peuvent aussi être définies dans le contexte non commutatif
des récurrences et équations différentielles. Certaines de ces généralisations seront
présentées dans le cas linéaire au Chapitre 30 et on trouvera plus de notes en page 558.

Bibliographie

Béz64 BÉZOUT, É. (1764). « Recherches sur le degré des équations résultantes de l'évanouissement des inconnues ». In : *Histoire de l'académie royale des sciences*, p. 288-338.

BGY80 BRENT, Richard P., Fred G. GUSTAVSON et David Y. YUN (1980). « Fast solution of Toeplitz systems of equations and computation of Padé approximants ». In : *Journal of Algorithms*, vol. 1, n°3, p. 259-295.
Bro71 BROWN, W. S. (1971). « On Euclid's algorithm and the computation of polynomial greatest common divisors ». In : *SYMSAC'71 : ACM Symposium on Symbolic and Algebraic Manipulation*. Los Angeles, California, United States : ACM, p. 195-211.

BT71 BROWN, W. S. et J. F. TRAUB (1971). « On Euclid's algorithm and the theory of subresultants ». In : *Journal of the Association for Computing Machinery*, vol. 18, p. 505-514.

Col66 COLLINS, G. E. (1966). « Polynomial remainder sequences and determinants ». In : *The American Mathematical Monthly*, vol. 73, p. 708-712.

Col67 COLLINS, George E. (1967). « Subresultants and reduced polynomial remainder sequences ». In : *Journal of the Association for Computing Machinery*, vol. 14, p. 128-142.

Col71 — (1971). « The calculation of multivariate polynomial resultants ». In : *Journal of the Association for Computing Machinery*, vol. 18, p. 515-532.

El 03 EL KAHOUI, M'hammed (2003). « An elementary approach to subresultants theory ». In : *Journal of Symbolic Computation*, vol. 35, n°3, p. 281-292.

GCL92 GEDDES, Keith O., Stephen R. CZAPOR et George LABAHN (1992). *Algorithms for computer algebra*. Kluwer Academic Publishers.

GG03 GATHEN, Joachim von zur et Jürgen GERHARD (2003). *Modern computer algebra*. 2e éd. Cambridge University Press.

GG13 — (2013). *Modern computer algebra*. 3e éd. Cambridge University Press.

GL03 GATHEN, Joachim von zur et Thomas LÜCKING (2003). « Subresultants revisited ». In : *Theoretical Computer Science*, vol. 297, n°1-3, p. 199-239.

GY79 GUSTAVSON, Fred G. et David Y. Y. YUN (1979). « Fast algorithms for rational Hermite approximation and solution of Toeplitz systems ». In : *IEEE Transactions on Circuits and Systems*, vol. 26, n°9, p. 750-755.

Knu71 KNUTH, Donald E. (1971). « The analysis of algorithms ». In : *Actes du Congrès International des Mathématiciens 1970*. Vol. 3. Nice : Gauthier-Villars, p. 269-274.

Knu97 — (1997). *The art of computer programming*. 3e éd. Vol. 2 : Numerical Algorithms. Computer Science and Information Processing. Addison-Wesley Publishing Co.

Lan02 LANG, Serge (2002). *Algebra*. 3e éd. Vol. 211. Graduate Texts in Mathematics. Springer-Verlag.

Las03 LASCOUX, Alain (2003). *Symmetric functions and combinatorial operators on polynomials*. Vol. 99. CBMS Regional Conference Series in Mathematics. Published for the Conference Board of the Mathematical Sciences, Washington, DC.

Leh38 LEHMER, D. H. (1938). « Euclid's Algorithm for large numbers ». In : *The American Mathematical Monthly*, vol. 45, n°4, p. 227-233.

Moe73 MOENCK, R. T. (1973). « Fast computation of GCDs ». In : *STOC'73 : ACM Symposium on Theory of Computing*. Austin : ACM, p. 142-151.
Möl08 MÖLLER, Niels (2008). « On Schönhage's algorithm and subquadratic integer GCD computation ». In: Mathematics of Computation, vol. 77, n°261, p. 589-607.

PW02 PAN, V. Y. et X. WANG (2002). « Acceleration of Euclidean algorithm and extensions ». In: ISSAC'02: International Symposium on Symbolic and Algebraic Computation. Éd. par Teo MORA. Lille: ACM Press, p. 207-213.

Sch71 SCHÖNHAGE, A. (1971). « Schnelle Berechnung von Kettenbruchentwicklungen ». In: Acta Informatica, vol. 1, p. 139-144.

Sch80b SCHWARTZ, J. T. (1980). « Fast probabilistic algorithms for verification of polynomial identities ». In: Journal of the Association for Computing Machinery, vol. 27, n°4, p. 701-717.

Sch88b SCHÖNHAGE, A. (1988). « Probabilistic computation of integer polynomial GCDs ». In: Journal of Algorithms, vol. 9, n°3, p. 365-371.

Str81 STRASSEN, V. (1981). « The computational complexity of continued fractions ». In: SYMSAC'81. Snowbird, Utah, United States: ACM, p. 51-67.

Str83 — (1983). « The computational complexity of continued fractions ». In: SIAM Journal on Computing, vol. 12, n°1, p. 1-27.

Syl39 SYLVESTER, James Joseph (1839). « On rational derivation from equations of coexistence, that is to say, a new and extended theory of elimination, 1839-40 ». In: The Collected mathematical papers of James Joseph Sylvester. Baker, H. F., p. 40-53.

Syl40 — (1840). « A method of determining by mere inspection the derivatives from two equations of any degree ». In: Philosophical Magazine Series 3, vol. 16, n°101, p. 132-135.

SZ04 STEHLÉ, D. et P. ZIMMERMANN (2004). « A binary recursive GCD algorithm ». In: ANTS-VI. Vol. 3076. Lecture Notes in Computer Science. Springer-Verlag, p. 411-425.

Yap00 YAP, Chee (2000). Fundamental problems in algorithmic algebra. Oxford University Press.

Yun76 YUN, David Y. Y. (1976). « On square-free decomposition algorithms ». In: SYMSAC'76: ACM Symposium on Symbolic and Algebraic Computation. Yorktown Heights, New York, United States: ACM, p. 26-35.

Yun77 — (1977). « On the equivalence of polynomial GCD and square-free factorization problems ». In: MACSYMA Users' Conference. NASA, Langley Res. Center, Washington D.C., United States, p. 65-70.
# 7. Approximants de Padé et de Padé-Hermite

## Résumé

Un premier problème abordé dans ce chapitre est le calcul d'approximants de Padé. Plus généralement, on s'intéresse à la *reconstruction rationnelle*, dont un autre cas particulier important est l'interpolation des fractions rationnelles. En Section 7.1, nous ramenons la reconstruction rationnelle à l'algorithme d'Euclide étendu, et l'appliquons à la reconnaissance d'une suite récurrente linéaire à partir de ses premiers termes (algorithme de Berlekamp-Massey.) Les applications en sont nombreuses : pour la résolution de systèmes différentiels à coefficients constants (Chapitre 13), pour le calcul de polynômes minimaux de matrices creuses (Chapitre 9), pour la résolution de systèmes linéaires à coefficients polynomiaux (Chapitre 11).

Un second problème, traité en Section 7.2, est le calcul d'approximants de Padé-Hermite. Il peut s'effectuer grâce à un algorithme qui généralise le calcul des approximants de Padé. L'approximation de Padé-Hermite permet de reconstruire des polynômes annulant une série algébrique, et des opérateurs différentiels annulant une série différentiellement finie. Il permet également de reconstruire une récurrence à coefficients polynomiaux d'ordre arbitraire à partir des premiers termes de la suite.

## 7.1 Reconstruction rationnelle

**Définition 7.1** Soit $\mathbb{K}$ un corps, $A \in \mathbb{K}[X]$ un polynôme de degré $n > 0$ et $B \in \mathbb{K}[X]$ de degré $< n$. Pour un $k \in \{1, \dots, n\}$ fixé, la *reconstruction rationnelle de B modulo A*
est la recherche d'un couple de polynômes $(R, V) \in \mathbb{K}[X]^2$ vérifiant :

$$ \text{pgcd}(V, A) = 1, \quad \text{deg}(R) < k, \quad \text{deg}(V) \le n-k \quad \text{et} \quad \frac{R}{V} \equiv B \quad \text{mod } A. \quad (\text{RR}) $$

On parle d'approximation de Padé lorsque $A = X^n$, et d'interpolation de Cauchy lorsque $A = \prod_i (X - u_i)$ avec $u_1, \dots, u_n \in \mathbb{K}$ deux à deux distincts.

<figure><img src="image_26.png" /></figure>

R Si $k=n$, alors clairement $(R, V) = (B, 1)$ est une solution du problème (RR).

Si $k < n$, il est possible que (RR) n'admette aucune solution. C'est le cas par exemple en prenant $n=3, k=2$ et $A = X^3, B = X^2 + 1 \in \mathbb{K}[X]$. En effet, si $V(X) = aX + b$ avec $b \neq 0$, alors $R = (aX + b)(X^2 + 1) = bX^2 + aX + b \pmod{X^3}$, ce qui est incompatible avec $\text{deg}(R) \le 1$.

Par ailleurs, si (RR) admet une solution $(R, V)$, alors la fraction rationnelle $R/V$ est forcément unique. En effet, si $(R_1, V_1) \in \mathbb{K}[X]^2$ est une autre solution de (RR), alors $R_1/V_1 \equiv R/V \pmod{A}$, donc $A$ divise $R_1V - V_1R$. Or, le polynôme $R_1V - V_1R$ ayant un degré strictement inférieur à celui de $A$, il doit être identiquement nul. Donc les fractions $R/V$ et $R_1/V_1$ coïncident.

### Un problème plus simple

Si $R, V \in \mathbb{K}[X]$ sont tels que $(R, V)$ est solution du problème (RR), alors $(R, V)$ vérifie aussi le problème plus simple

$$ \text{deg}(R) < k, \quad \text{deg}(V) \le n-k \quad \text{et} \quad R \equiv VB \pmod{A}, \quad (\text{RRS}) $$

où, à la différence de (RR), on a mis de côté la contrainte sur le pgcd de $A$ et de $V$.

Remarquons tout de suite que le problème (RRS) admet toujours une solution non triviale $(R, V) \neq (0, 0)$ : en effet, il se traduit en termes d'algèbre linéaire en un système linéaire homogène ayant $k + (n - k + 1) = n + 1$ inconnues (les coefficients de $R$ et de $V$) et $n$ équations. Par ailleurs, la solution $R/V$ du problème (RRS) est encore unique (même idée de preuve que pour (RR) : si $A$ divise à la fois $R - VB$ et $R_1 - V_1B$, alors il divise également la combinaison linéaire $RV_1 - R_1V = (R - VB)V_1 - (R_1 - V_1B)V$).

### Lien avec le pgcd étendu

Nous allons prouver que le problème (RRS) peut être résolu en utilisant l'algorithme d'Euclide étendu décrit en page 116. Nous en déduirons ensuite une procédure de décision et calcul pour (RR).

L'intuition du lien entre le problème (RRS) et le calcul du pgcd étendu s'appuie sur la remarque simple suivante : la congruence $R \equiv VB \pmod{A}$ équivaut à l'existence d'un polynôme $U$ tel que $R = UA + VB$; or, cette dernière égalité est une relation de type Bézout.

Par exemple, dans le cas particulier $k=1$, si $A$ et $B$ sont premiers entre eux, alors on prend $R=1$ et la relation de Bézout $R = UA + VB$ avec $\text{deg}(V) < \text{deg}(A) = n$ fournit la réponse. Sinon, on prend $R=0$, et comme le ppcm de $A$ et $B$ a un degré $< m+n$, le polynôme $V = \text{ppcm}(A, B)/B$ est de degré au plus $n-1$ et vérifie $VB = 0 \pmod{A}$.
## Calcul de la reconstruction rationnelle

L'algorithme d'Euclide étendu, décrit page 116, calcule une suite de restes successifs dont le degré décroît, ainsi qu'une suite de cofacteurs $(U_i, V_i)$. Les éléments de la i-ième itération vérifient l'identité de Bézout $U_iA + V_iB = R_i$. L'écriture matricielle de l'itération

$$ \begin{pmatrix} U_i & V_i \\ U_{i+1} & V_{i+1} \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ 1 & -Q_i \end{pmatrix} \times \begin{pmatrix} U_{i-1} & V_{i-1} \\ U_i & V_i \end{pmatrix} $$

permet de déduire aisément que $U_iV_{i+1} - V_iU_{i+1} = (-1)^i$ quel que soit $i$, et en particulier que $U_i$ et $V_i$ sont premiers entre eux. Cela entraîne également que $\text{pgcd}(R_i, V_i) = \text{pgcd}(A, V_i)$.

La clé du résultat principal de cette section (Théorème 7.2 ci-dessous) est le lemme suivant, qui montre que les cofacteurs $V_i$ obtenus par l'algorithme d'Euclide étendu ont des degrés modérés et bien maîtrisés.

**Lemme 7.1** Soient $n_0 = \text{deg}(A) > n_1 = \text{deg}(B) > n_2 > \dots > n_\ell$ les degrés des restes $R_i$ dans l'algorithme d'Euclide étendu exécuté avec A et B. Alors :

$$ \text{deg}(V_i) = n - n_{i-1} \quad \text{pour } 1 \le i \le \ell + 1. \qquad (7.1) $$

*Démonstration.* On procède par récurrence forte sur $i$. L'initialisation est claire, car $\text{deg}(V_1) = 0$. Supposons l'égalité (7.1) vérifiée pour $1 \le i \le k$.

Comme $\text{deg}(V_{k-1}) = n - n_{k-2} < n - n_{k-1} = \text{deg}(V_k)$, le degré de $V_{k+1} = V_{k-1} - Q_k V_k$ est égal à $\text{deg}(Q_k V_k) = \text{deg}(Q_k) + \text{deg}(V_k) = (n_{k-1} - n_k) + (n - n_{k-1}) = n - n_k$. Ceci donne le résultat pour $i = k + 1$ et conclut la preuve. ■

R

Dans une situation « générique », les quotients successifs $Q_i$ ont tous degré 1, et donc $n_{i+1} = n_i - 1$ (la suite des restes est appelée *normale*), si bien que le lemme précédent implique alors l'égalité $\text{deg}(V_i) = i - 1$.

Le résultat suivant montre que problème de reconstruction rationnelle (RR) peut être résolu à l'aide de l'algorithme d'Euclide étendu.

**Théorème 7.2** Soit $A \in \mathbb{K}[X]$ de degré $n > 0$ et soit $B \in \mathbb{K}[X]$ de degré $< n$. Soit $k \in \{1, 2, \dots, n\}$ et soient $R_j, U_j, V_j$ le reste et les cofacteurs obtenus à la j-ième itération de l'algorithme d'Euclide étendu exécuté sur A et B. Si $j$ est choisi minimal tel que $\text{deg}(R_j) < k$, alors on a :

1. Il existe une solution $(R, V) \neq (0, 0)$ de (RRS), à savoir $(R, V) = (R_j, V_j)$. Si, de plus, $\text{pgcd}(R_j, V_j) = 1$, alors $(R, V)$ est aussi solution de (RR).

2. Si (RR) admet une solution et si $R/V \in \mathbb{K}(X)$ en est une forme irréductible, alors il existe une constante $\alpha \in \mathbb{K} \setminus \{0\}$ telle que $R = \alpha R_j$ et $V = \alpha V_j$.

Le problème (RR) admet donc une solution si, et seulement si, $\text{pgcd}(A, V_j) = 1$.

*Démonstration.* Par construction, $R_j$ vaut bien $U_jA + V_jB \equiv V_jB$ modulo A. Par ailleurs, le Lemme 7.1 montre que le degré de $V_j$ vaut $n - \text{deg}(R_{j-1})$ et est donc borné par $n - k$ (par la minimalité de $j$). Donc $(R, V) = (R_j, V_j)$ vérifie bien (RRS).
**Entrée** A, B ∈ K[X] avec n = deg(A) > deg(B), et k ∈ {1, ..., n}.

**Sortie** Une solution (R, V) de (RR), ou 0 si une telle solution n'existe pas.

1. $R_0 = A$; $U_0 = 1$; $V_0 = 0$; $R_1 = B$; $U_1 = 0$; $V_1 = 1$; $i = 1$.

2. Tant que $\deg(R_i) \ge k$, faire :

a. Calculer la quotient $Q_i$ et le reste $R_{i+1}$ de la division de $R_{i-1}$ par $R_i$.

b. $U_{i+1} = U_{i-1} - Q_i U_i$; $V_{i+1} = V_{i-1} - Q_i V_i$.

c. $i = i + 1$.

3. Si $\text{pgcd}(A, V_i) = 1$, renvoyer $(R_i, V_i)$; sinon renvoyer 0.

**Algorithme 7.1 – Algorithme de reconstruction rationnelle.**

De plus, le $\text{pgcd}(A, V_j)$ est égal au $\text{pgcd}(R_j, V_j)$. Par conséquent, si ce dernier vaut 1, alors $V_j$ est inversible modulo A et donc $(R, V) = (R_j, V_j)$ vérifie aussi (RR).

L'unicité (à une constante près) est plus délicate à prouver. Supposons que $(R, V)$ est une solution de (RR). Alors R s'écrit UA + VB pour un certain $U \in K[X]$. Nous allons prouver que $(R, V) = (\alpha R_j, \alpha V_j)$ pour un certain $\alpha \in K[X] \setminus \{0\}$.

D'abord, les polynômes $U_j$ et $V_j$ vérifient $U_jV = UV_j$ : supposons le contraire et considérons le système linéaire

$$ \begin{pmatrix} U_j & V_j \\ U & V \end{pmatrix} \times \begin{pmatrix} A \\ B \end{pmatrix} = \begin{pmatrix} R_j \\ R \end{pmatrix}. $$

Le système étant supposé de Cramer, A s'écrit comme quotient de deux déterminants

$$ A = \begin{vmatrix} R_j & V_j \\ R & V \end{vmatrix} \begin{vmatrix} U_j & V_j \\ U & V \end{vmatrix}. $$

L'hypothèse $\deg(R_j) < k \le \deg(R_{j-1})$ et le Lemme 7.1 entraînent que le degré du membre droit de l'égalité précédente est majoré par

$$ \begin{align*} \deg(R_j V - R V_j) &\le \max\{k-1 + \deg(V), \deg(R) + n - \deg(R_{j-1})\} \\ &\le \max\{n-1, (k-1) + n - \deg(R_{j-1})\} = n-1, \end{align*} $$

ce qui contredit $\deg(A) = n$, concluant la preuve de $U_j V = UV_j$. Cette égalité implique que $V_j$ divise $U_j V$, et puisque $U_j$ et $V_j$ sont premiers entre eux, il s'ensuit que $V_j$ divise V. En écrivant $V = \alpha V_j$, avec $\alpha \in K[X]$, on obtient $UV_j = VU_j = \alpha U_j V_j$, et donc $U = \alpha U_j$. Enfin, $R = UA + VB = \alpha(U_j A + V_j B) = \alpha R_j$. Les polynômes R et V étant premiers entre eux, il s'ensuit que $\alpha$ est une constante de $K \setminus \{0\}$, et que donc $\text{pgcd}(R_j, V_j) = 1$. ■

L'Algorithmme 7.1 résume cette approche ; la ressemblance avec l'algorithme d'Euclide étendu est claire.
**Entrée** Une borne $N \in \mathbb{N}$ sur le degré du polynôme minimal de la suite $(a_n)_{n \ge 0}$ et les $2N$ premiers termes $a_0, \dots, a_{2N-1} \in \mathbb{K}$.

**Sortie** Le polynôme générateur minimal de $(a_n)_{n \ge 0}$.

1. $A = a_0 + a_1 X + \dots + a_{2N-1} X^{2N-1}$.

2. Calculer $R, V \in \mathbb{K}[X]$ la solution de (7.3) telle que $V(0) = 1$.

3. $d = \max\{1 + \deg(R), \deg(V)\}$. Renvoyer $\tilde{V} = V(1/X)X^d$.

**Algorithme 7.2 – Algorithme de Berlekamp–Massey.**

### Approximants de Padé

Si $n > 0, k \in \{1, 2, \dots, n\}$ et si $B \in \mathbb{K}[X]$ est de degré $< n$, alors un approximant de Padé pour $B$ de type $(k-1, n-k)$ est une fraction rationnelle $R/V \in \mathbb{K}[X]$ telle que

$$X \nmid V, \quad \deg(R) < k, \quad \deg(V) \le n-k \quad \text{et} \quad \frac{R}{V} \equiv B \quad \text{mod} \ X^n. \qquad (7.2)$$

Le Théorème 7.2 (pour $A = X^n$) entraîne la conséquence suivante.

**Corollaire 7.3** Soit $B \in \mathbb{K}[X]$ de degré $< n$. Soient $R_j, U_j, V_j$ le reste et les cofacteurs obtenus à la $j$-ième itération de l'algorithme d'Euclide étendu exécuté avec $X^n$ et B. Si $j$ est choisi minimal tel que $\deg(R_j) < k$, alors on a :

1. L'équation (7.2) admet une solution si et seulement si $\text{pgcd}(R_j, V_j) = 1$.

2. Si $\text{pgcd}(R_j, V_j) = 1$, alors $R_j/V_j$ est l'unique approximant de Padé pour B de type $(k-1, n-k)$.

L'algorithme qui s'en déduit est un cas particulier de l'Algorithme 7.1, en prenant $A = X^n$ et en remplaçant le test $\text{pgcd}(A, V_i) = 1$ de l'étape (3) par $V_i(0) \ne 0$.

### Algorithme de Berlekamp–Massey

L'algorithme décrit dans cette section permet de deviner des récurrences à coefficients constants d'ordre arbitraire à partir de termes consécutifs de la suite.

On se donne dans un corps $\mathbb{K}$ les $2N$ premiers éléments d'une suite récurrente linéaire $(a_n)_{n \ge 0}$ pour laquelle on sait qu'il existe un polynôme caractéristique de degré $d \le N$, c'est-à-dire un polynôme

$$f(X) = f_d X^d + \dots + f_0, \quad \text{tel que} \quad f_d a_{n+d} + \dots + f_0 a_n = 0, \quad \text{pour tout } n \ge 0.$$

Le problème est de calculer le polynôme minimal (c'est-à-dire le polynôme caractéristique de degré minimal) de $(a_n)_{n \ge 0}$. Le Lemme 4.7 (page 87) montre que ce calcul équivaut à la résolution du problème de Padé

$$\frac{R}{V} \equiv A \quad \text{mod} \ X^{2N}, \quad X \nmid V, \deg(R) < N, \deg(V) \le N \quad \text{et} \quad \text{pgcd}(R, V) = 1. \qquad (7.3)$$

En effet, il implique que $(R, V) = (N_0, \bar{P})$ est solution de (7.3). La méthode qui en découle est résumée dans l'Algorithme 7.2.
### Interpolation rationnelle de Cauchy

L'interpolation des fractions rationnelles, appelée aussi interpolation de Cauchy, est une généralisation naturelle de l'interpolation polynomiale (de Lagrange).

Soient $k \in \{1, \dots, n\}$, $v_0, \dots, v_{n-1} \in \mathbb{K}$ et soient $u_0, \dots, u_{n-1}$ des points distincts de $\mathbb{K}$. On cherche une fraction rationnelle $R/V \in \mathbb{K}(X)$ avec $R, V \in \mathbb{K}[X]$ satisfaisant aux contraintes

$$V(u_i) \neq 0, \quad \frac{R(u_i)}{V(u_i)} = v_i \quad \text{pour } 0 \le i < n, \quad \deg(R) < k, \deg(V) \le n-k. \quad (7.4)$$

Le Théorème 7.2 implique la procédure de décision et de calcul suivante : on calcule le polynôme interpolant $B$ et le polynôme $A = (X - u_0) \cdots (X - u_{n-1})$. On exécute l'algorithme d'Euclide étendu sur $A$ et $B$ : on note $R_j, U_j, V_j$ le reste et les cofacteurs obtenus à la $j$-ième itération. On choisit ensuite $j$ minimal tel que $\deg(R_j) < k$. Si $\text{pgcd}(A, V_j) = 1$, alors $R_j/V_j$ est l'unique forme canonique de la fraction rationnelle cherchée. Sinon, le problème (7.4) n'a pas de solution.

Une application importante de l'interpolation rationnelle est la reconnaissance de récurrences linéaires d'ordre 1 à coefficients polynomiaux.

**Corollaire 7.4 —** *Devinette de récurrences hypergéométriques.* Soit $(a_n)$ une suite d'éléments non nuls de $\mathbb{K}$ vérifiant la récurrence $p_1(n)a_{n+1} + p_0(n)a_n = 0$, dont les coefficients (inconnus) $p_0(X)$ et $p_1(X)$ sont des polynômes de $\mathbb{K}[X]$, premiers entre eux (récurrence hypergéométrique). Si l'on connaît une borne $d$ sur le degré de $p_0$ et de $p_1$ telle que $p_0(j) \neq 0$ pour $0 \le j \le 2d+1$, alors on peut déterminer $p_0$ et $p_1$ par un calcul d'interpolation rationnelle, à partir des $2d+1$ premiers termes de la suite $(a_n)$.

### Algorithmes rapides

En combinant les résultats du Théorème 7.2 ci-dessus, du Théorème 6.16 en page 134, et du Chapitre 5 pour le calcul rapide des polynômes $A$ et $B$ définis en Section 7.1, on déduit le résultat de complexité suivant.

**Théorème 7.5** Soient $A \in \mathbb{K}[X]$ de degré $n > 0$, $B \in \mathbb{K}[X]$ de degré $< n$ et $k \in \{1, 2, \dots, n\}$. Soient $v_0, \dots, v_{n-1} \in \mathbb{K}$, et soient $u_0, \dots, u_{n-1}$ des points distincts de $\mathbb{K}$. Il est possible de calculer en $O(M(n) \log n)$ opérations dans $\mathbb{K}$ :

1. une solution du problème reconstruction rationnelle (RR);

2. un approximant de Padé pour $B$ de type $(k-1, n-k)$;

3. une fraction rationnelle $F = R/V \in \mathbb{K}(X)$ telle que $\deg(V) \le n-k$, $\deg(R) < k$, et $F(u_i) = v_i$ pour $0 \le i < n$.

### Le cas entier

L'an analogue pour les entiers de la reconstruction rationnelle est également calculable en complexité quasi-optimale. Si $n$ et $B$ sont deux entiers, $0 \le B < n$, et étant donné $k \le n$, il s'agit de calculer des entiers $R$ et $V$ tels que

$$\text{pgcd}(n, V) = 1 \quad \text{et} \quad B = \frac{R}{V} \mod n, \quad |R| < k, \quad 0 \le V \le \frac{n}{k}. \quad (7.5)$$
**Théorème 7.6** À partir de $k, n$ et B, on peut calculer R et V satisfaisant à (7.5) en $O(M_Z(\log n) \log \log n)$ opérations binaires.

### Le cas des récurrences

Une conséquence utile du Théorème 7.5 à la reconnaissance rapide des suites est contenue dans le corollaire suivant. Sa généralisation aux cas des suites polynomialement récurrentes arbitraires sera traitée en Section 7.2.

**Corollaire 7.7** Soit $(a_n)$ une suite d'éléments de $\mathbb{K}$ vérifiant

- soit une récurrence linéaire à coefficients constants d'ordre au plus N;

- soit une récurrence linéaire d'ordre 1 à coefficients polynomiaux de degré au plus N.

À partir des premiers $2N + 1$ termes de la suite, on peut retrouver les coefficients de la récurrence en $O(M(N) \log N)$ opérations dans $\mathbb{K}$.

## 7.2 Approximants de Padé–Hermite

**Définition 7.2** Soit $\mathbf{F} = {}^t(f_1, \dots, f_n)$ un vecteur de $n \ge 1$ séries formelles de $\mathbb{K}[[X]]$, et soit $\mathbf{d} = (d_1, \dots, d_n) \in \mathbb{N}^n$. Un vecteur non nul $\mathbf{P} = (P_1, \dots, P_n)$ de polynômes de $\mathbb{K}[X]$ est appelé *approximant de Padé-Hermite de type* $\mathbf{d}$ de $\mathbf{F}$ si :

1. La valuation $\text{val}(\mathbf{P} \cdot \mathbf{F})$ de la série $\mathbf{P} \cdot \mathbf{F} = \sum_{i=1}^n P_i f_i$ est au moins égale à $\sigma = \sum (d_i + 1) - 1$;

2. $\text{deg}(P_i) \le d_i$ pour tout $1 \le i \le n$.

L'entier $\sigma$ est alors appelé l'ordre de l'approximant.

**Exemple 7.1** Dans la terminologie de la section précédente, si $R/V \in \mathbb{K}(X)$ est un approximant de Padé de type $(k, n-k)$ de $B \in \mathbb{K}[[X]]$, alors $(R, V)$ est un approximant de Padé-Hermite pour $(-1, B)$, de type $(k, n-k)$.

**Exercice 7.1** Soient A, B deux polynômes de $\mathbb{K}[X]$, avec $n = \text{deg}(A) > \text{deg}(B)$. Montrer que $(R, V)$ est solution du problème de reconstruction rationnelle (RRS) défini en page 140 si et seulement s'il existe un polynôme $U \in \mathbb{K}[X]$ tel que $(R, V, U)$ soit un approximant de Padé-Hermite de $(-1, B, A)$ de type $(k-1, n-k, n-1)$. ■

Le premier résultat de cette section montre l'existence des approximants de Padé-Hermite et fournit également un premier algorithme pour leur calcul.

**Théorème 7.8** Tout vecteur de séries formelles $\mathbf{F} = (f_1, \dots, f_n) \in \mathbb{K}[[X]]^n$ admet un approximant de Padé-Hermite de type $\mathbf{d} = (d_1, \dots, d_n) \in \mathbb{N}^n$ donné.

*Démonstration.* On procède par coefficients indéterminés : en écrivant $P_i = \sum_{j=0}^{d_i} P_{i,j} X^j$, on obtient un système linéaire homogène à $\sigma = \sum_i (d_i + 1) - 1$ équations en les $\sigma + 1$
inconnues $p_{i,j}$. Puisqu'il a moins d'équations que d'inconnues, ce système admet forcément une solution non triviale.

L'algorithme qui découle de la preuve du Théorème 7.8 repose sur la recherche d'un élément non trivial dans le noyau d'une matrice à coefficients dans $K$, de taille $\sigma \times (\sigma + 1)$. En utilisant les résultats du Chapitre 8, on aboutit donc à un algorithme de complexité $O(MM(\sigma)) = O(\sigma^\theta)$ pour $2 \le \theta \le 3$.

### Algorithme de Derksen : idées et résultats préliminaires

Le but de la suite de ce chapitre est de présenter un algorithme plus efficace, dû à Derksen, de complexité seulement quadratique en $\sigma$. En anticipant un peu, cet algorithme revient à effectuer une sorte de pivot de Gauss sur une matrice à coefficients dans $K[X]$, mais de taille bien plus petite que $\sigma$ (seulement linéaire en $\max(\mathbf{d})$). Dans le cas particulier $n = 2$, l'algorithme de Derksen a essentiellement la même complexité que le calcul d'approximants de Padé via l'algorithme d'Euclide étendu vu en Section 7.1.

Pour simplifier la présentation, on se restreint dans la suite au cas où le type de l'approximant cherché est de la forme $\mathbf{d} = (d, ..., d) \in \mathbb{N}^n$, pour un certain $d \in \mathbb{N}$. L'idée de l'algorithme de Derksen est de construire non pas un seul approximant de Padé-Hermite, mais toute une famille de tels approximants, et cela de manière incrémentale. Plus exactement, pour $s = 0, 1, ..., $, il construit une base d'une forme spéciale, appelée « base minimale », du $K[X]$-module

$$V_s = \{ \mathbf{P} \in K[X]^n \mid \mathrm{val}(\mathbf{P} \cdot \mathbf{F}) \ge s \}.$$

Comme nous le verrons plus loin, une telle base de $V_\sigma$ pour $\sigma = nd + n - 1$ contiendra alors nécessairement un approximant de Padé-Hermite de type $\mathbf{d} = (d, ..., d)$ de $\mathbf{F}$.

Observons que, grâce aux inclusions $X^s K[X]^n \subseteq V_s \subseteq K[X]^n$, le module $V_s$ est libre de rang $n$. Précisons ce que l'on entend par « base minimale ». Pour ce faire, introduisons d'abord une notion de *degré* et de *type* d'un vecteur de polynômes.

**Définition 7.3 — Degré et type d'un vecteur de polynômes.** Pour un vecteur de polynômes $\mathbf{P} = (P_1, ..., P_n) \in K[X]^n$, son *degré* $\mathrm{deg}(\mathbf{P})$ et son *type* $\mathrm{type}(\mathbf{P})$ sont donnés par :

$$\mathrm{deg}(\mathbf{P}) = \max\{\mathrm{deg}(P_1), ..., \mathrm{deg}(P_n)\} \text{ et } \mathrm{type}(\mathbf{P}) = \max\{i \mid \mathrm{deg}(\mathbf{P}) = \mathrm{deg}(P_i)\}.$$

**Définition 7.4 — Base minimale.** Soit $V \subseteq K[X]^n$ un sous-module libre de rang $n$. Une suite $Q_1, ..., Q_n$ est appelée « base minimale » de $V$ si, pour tout $i$, le vecteur $Q_i$ est non nul, de type $i$, et de degré minimal parmi les éléments de $V \setminus \{0\}$ de type $i$.

Le résultat suivant précise le lien entre une base minimale et l'approximation de Padé-Hermite.
**Lemme 7.9** Soit $d \ge 1$. Supposons que $Q_1, \dots, Q_n$ est une base minimale de $V_{nd+n-1}$. Soit $\ell$ tel que $\deg(Q_\ell)$ est minimal. Alors $Q_\ell$ est un approximant de Padé–Hermite de type $(d, \dots, d)$ pour $F$.

*Démonstration.* Le Théorème 7.8 montre l'existence d'un vecteur non nul $\mathbf{P}$ de $V_{nd+n-1}$ tel que $\deg(\mathbf{P}) \le d$. Soit $i$ le type de $\mathbf{P}$. La suite d'inégalités

$$ \deg(Q_\ell) \le \deg(Q_i) \le \deg(\mathbf{P}) \le d, $$

prouve que $Q_\ell$ est un approximant de Padé–Hermite de type $(d, \dots, d)$ de $F$. $\blacksquare$

Le résultat suivant montre qu'une base minimale est nécessairement une base du $\mathbb{K}[X]$-module $V$ au sens usuel.

**Théorème 7.10** Soit $V \subseteq \mathbb{K}[X]^n$ un sous-module libre de rang $n$. Toute base minimale de $V$ est une base du $\mathbb{K}[X]$-module $V$.

*Démonstration.* Montrons d'abord qu'il s'agit d'un système de générateurs. Soit $W = \mathbb{K}[X]Q_1 + \dots + \mathbb{K}[X]Q_n \subseteq V$. On suppose par l'absurde que $V \ne W$. Soit $\mathbf{P} \in V \setminus W$ un élément minimal dans $V \setminus W$ pour l'ordre $\mathbf{P} < \mathbf{Q}$ défini par

$$ \begin{gathered} \deg(\mathbf{P}) < \deg(\mathbf{Q}), \quad \text{ou} \\ \deg(\mathbf{P}) = \deg(\mathbf{Q}) \quad \text{et} \quad \text{type}(\mathbf{P}) < \text{type}(\mathbf{Q}). \end{gathered} \tag{7.6} $$

Autrement dit, $\mathbf{P}$ est de type minimal parmi les éléments de degré minimal de $V \setminus W$.

Soit $i$ le type de $\mathbf{P}$. Les relations type($\mathbf{P}$) = type($Q_i$) et $\deg(\mathbf{P}) \ge \deg(Q_i)$ entraînent l'existence d'un monôme $q \in \mathbb{K}[X]$ de degré $\deg(q) = \deg(\mathbf{P}) - \deg(Q_i)$ tel que type($\mathbf{P} - qQ_i$) < type($\mathbf{P}$). Puisque $\deg(\mathbf{P} - qQ_i) \le \deg(\mathbf{P})$, on obtient que

$$ \mathbf{P} - q\mathbf{Q}_i < \mathbf{P}. $$

Par la minimalité de $\mathbf{P}$, il s'ensuit que $\mathbf{P} - q\mathbf{Q}_i$ appartient à $W$, donc $\mathbf{P} \in W$, ce qui contredit le choix de $\mathbf{P}$.

Pour conclure la preuve, montrons que les $Q_i$ forment une famille libre. Si $\sum_i a_i Q_i = 0$ est une combinaison polynomiale nulle des $Q_i$, on a que pour tout $i$, le vecteur $a_i Q_i$ est de type $i$. L'assertion découle du lemme suivant. $\blacksquare$

**Lemme 7.11** Si $\mathbf{P}$ et $\mathbf{Q}$ ont des types distincts, alors le type de $\mathbf{P} + \mathbf{Q}$ est soit le type de $\mathbf{P}$ soit celui de $\mathbf{Q}$.

*Démonstration.* Supposons $j = \text{type}(\mathbf{P}) > i = \text{type}(\mathbf{Q})$. Si $\deg(\mathbf{P}) \ge \deg(\mathbf{Q})$, alors on a $\text{type}(\mathbf{P} + \mathbf{Q}) = j$ et si $\deg(\mathbf{P}) < \deg(\mathbf{Q})$, puis $\text{type}(\mathbf{P} + \mathbf{Q}) = i$. $\blacksquare$
## Algorithme de Derksen : fonctionnement

L'idée de l'algorithme est de construire de proche en proche une base minimale de $V_s$, partant de la base minimale des $Q_k = (0, 0, ..., 0, 1, 0, ..., 0)$ (avec 1 en position $k$) de $V_0$. D'après le Lemme 7.9, l'élément de degré minimal dans une base minimale de $V_{nd+n-1}$ fournit un approximant de Padé-Hermite de type $(d, ..., d)$ de $F$.

Le résultat suivant montre comment construire une base minimale de $V_{s+1}$ à partir d'une base minimale de $V_s$, et il constituera le cœur de l'itération.

**Théorème 7.12** Soit $Q_1, ..., Q_n$ une base minimale de

$$V_s = \{ \mathbf{P} \in \mathbb{K}[X]^n \mid \mathrm{val}(\mathbf{P} \cdot \mathbf{F}) \geq s \}.$$

1. Si $\mathrm{val}(Q_i \cdot F) \geq s+1$ quel que soit $i$, alors $V_{s+1}$ et $V_s$ coïncident, et $\{Q_1, ..., Q_n\}$ est une base minimale de $V_{s+1}$.

2. Supposons que $1 \leq i \leq n$ est tel que les deux conditions suivantes soient réunies :

— $\mathrm{val}(Q_i \cdot F) = s$;

— si $\mathrm{val}(Q_\ell \cdot F) = s$ pour un $\ell \neq i$, alors $Q_i < Q_\ell$, où $\ell < s$ est l'ordre (7.6).

Alors :

a. pour $\ell \neq i$, il existe un scalaire $\lambda_\ell \in \mathbb{K}$ tel que $\tilde{Q}_\ell = Q_\ell - \lambda_\ell Q_i$ vérifie $\mathrm{val}(\tilde{Q}_\ell \cdot F) > s$;

b. en posant $\tilde{Q}_i = XQ_i$, la suite $\tilde{Q}_1, ..., \tilde{Q}_n$ forme une base minimale de $V_{s+1}$.

**Démonstration.** L'inclusion $V_{s+1} \subseteq V_s$ est évidente, et inversement, le Théorème 7.10 montre que tout $\mathbf{P} \in V_s$ s'écrit comme combinaison linéaire $\sum_i a_i Q_i$. Ainsi $\mathbf{P} \cdot \mathbf{F} = \sum_i a_i (Q_i \cdot \mathbf{F})$ est de valuation au moins $s+1$, donc $\mathbf{P} \in V_{s+1}$. Ceci prouve le premier point.

Prouvons maintenant l'assertion (2.a). Si $\mathrm{val}(Q_\ell \cdot F) > s$, on pose $\lambda_\ell = 0$; si $\mathrm{val}(Q_\ell \cdot F) = s$, alors $Q_\ell \cdot F = c_\ell X^s + \dots$ et $Q_i \cdot F = c_i X^s + \dots$, avec $c_i \neq 0$, et alors $\lambda_\ell = c_\ell/c_i$ convient.

Pour l'assertion (2.b), commençons par montrer que

$$\tilde{Q}_1, \dots, \tilde{Q}_{i-1}, Q_i, \tilde{Q}_{i+1}, \dots, \tilde{Q}_n$$

reste une base minimale de $V_s$. Il suffit pour cela de montrer que pour $\ell \neq i$, le vecteur $\tilde{Q}_\ell$ a même type et même degré que $Q_\ell$. Si $\mathrm{val}(Q_\ell \cdot F) > s$, cela est évident, car $\lambda_\ell = 0$ et donc $\tilde{Q}_\ell = Q_\ell$. Sinon, le choix de $i$ assure que $Q_\ell > Q_i$, et donc $Q_\ell$ et $Q_\ell - \lambda_\ell Q_i$ ont même degré et même type.

Montrons maintenant que $(\tilde{Q}_1, ..., \tilde{Q}_n)$ est une base minimale de $V_{s+1}$. Comme la multiplication par un polynôme ne change pas le type, celui de $\tilde{Q}_i = XQ_i$ est bien $i$. Il suffit donc de montrer que si $\mathbf{P} \in V_{s+1}$ est de type $\ell \in \{1, 2, ..., n\}$, alors $\mathrm{deg}(\mathbf{P}) \geq \mathrm{deg}(\tilde{Q}_\ell)$. Si $\ell \neq i$, ceci est une évidence : comme $\mathbf{P}$ appartient à $V_{s+1} \subseteq V_s$, et comme la suite $\tilde{Q}_1, ..., \tilde{Q}_{i-1}, Q_i, \tilde{Q}_{i+1}, ..., \tilde{Q}_n$ forme une base minimale de $V_s$, le degré de $\mathbf{P}$ est nécessairement au moins égal à $\mathrm{deg}(\tilde{Q}_\ell) = \mathrm{deg}(Q_\ell)$.

Dans la suite de la preuve, on peut donc supposer que $\mathbf{P} \in V_{s+1}$ est de type $i$, le but étant de montrer que $\mathrm{deg}(\mathbf{P}) \geq \mathrm{deg}(XQ_i)$. Le Théorème 7.10 montre que $\mathbf{P}$ s'écrit
**Entrée** $\mathbf{F} = (f_1, \dots, f_n) \in \mathbb{K}[[X]]^n$ et $d \ge 1$.

**Sortie** Un approximant de Padé-Hermite de $\mathbf{F}$, de type $(d, \dots, d)$.

1. Pour $k$ de $1$ à $n$ définir

$\mathbf{Q}_k = (0, 0, \dots, 0, 1, 0, \dots, 0)$, avec 1 en position $k$.

2. Pour $j$ de $0$ à $nd + n - 2$ faire :

a. $i = 0$

b. Pour $k$ de $1$ à $n$ faire :

i. $c_k = \text{coeff}(\mathbf{Q}_k \cdot \mathbf{F}, j)$

ii. Si $c_k \neq 0$ et $(\mathbf{Q}_k < \mathbf{Q}_i \text{ ou } i = 0)$, alors $i = k$.

c. Si $i \neq 0$ alors faire

i. $\mathbf{Q}_i = c_i^{-1} \mathbf{Q}_i$

ii. Pour $k$ de $1$ à $n$ faire

si $k \neq i$, alors $\mathbf{Q}_k = \mathbf{Q}_k - c_k \mathbf{Q}_i$.

iii. $\mathbf{Q}_i = X\mathbf{Q}_i$

3. $p = 1$

4. Pour $k$ de $2$ à $n$ faire :

si $\text{deg}(\mathbf{Q}_k) < \text{deg}(\mathbf{Q}_p)$, alors $p = k$.

5. Renvoyer $\mathbf{Q}_p$.

**Algorithme 7.3** – Algorithme de Derksen pour les approximants de Padé-Hermite.

$\mathbf{P} = \sum_{j \neq i} a_j \tilde{\mathbf{Q}}_j + a_i \mathbf{Q}_i$. Puisque $\text{type}(\mathbf{P}) = i$, le Lemme 7.11 entraîne $a_i \neq 0$. De plus, le degré de $\mathbf{P}$ est égal à celui de $a_i \mathbf{Q}_i$, d'après le Lemme 7.13 ci-dessous.

Comme $\text{val}(\mathbf{P} \cdot \mathbf{F}) > s$ et comme pour $k \neq i$, $\text{val}(\mathbf{Q}_k \cdot \mathbf{F}) > s$, on a nécessairement que $\text{val}(a_i) > 0$. En particulier $\text{deg}(a_i) > 0$ et donc $\text{deg}(\mathbf{P}) \geq 1 + \text{deg}(\mathbf{Q}_i)$. $\blacksquare$

**Lemme 7.13** Si $\text{type}(\mathbf{P}) = \text{type}(\mathbf{Q}) = i$ et $\text{type}(\mathbf{P} + \mathbf{Q}) < i$, alors $\text{deg}(\mathbf{P}) = \text{deg}(\mathbf{Q})$.

*Démonstration.* Soient $\mathbf{P} = (P_1, \dots, P_n)$ et $\mathbf{Q} = (Q_1, \dots, Q_n)$. L'hypothèse entraîne les égalités $\text{deg}(P_i) = \text{deg}(\mathbf{P})$ et $\text{deg}(Q_i) = \text{deg}(\mathbf{Q})$. Cela implique que $P_i$ et $Q_i$ ont le même degré; sinon, le type de $\mathbf{P} + \mathbf{Q}$ serait égal à $i$. $\blacksquare$

L'Algorithme 7.3 découle de la conjonction du Lemme 7.9 et du Théorème 7.12. Il est de complexité seulement quadratique en $\sigma$. En résumé, nous avons le résultat suivant.

**Théorème 7.14** Soit $\mathbf{F} = (f_1, \dots, f_n) \in \mathbb{K}[[X]]^n$. Il est possible de calculer un approximant de Padé-Hermite de type $(d, \dots, d)$ de $\mathbf{F}$ en $O(n\sigma^2) = O(n^3d^2)$ opérations dans $\mathbb{K}$.

**Exercice 7.2** Écrire les détails de la preuve du Théorème 7.14. $\blacksquare$
R

Dans le cas « générique » (dit normal), la sortie $\mathbf{Q} = \mathbf{Q}_1, \dots, \mathbf{Q}_n$ de l'algorithme de Derksen est de degré

$$ \left( \begin{array}{ccccc} d+1 & d & \cdots & d & d \\ d+1 & d+1 & \cdots & d & d \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ d+1 & d+1 & \cdots & d+1 & d \\ d & d & \cdots & d & d \end{array} \right). $$

C'est la dernière ligne $\mathbf{Q}_n$ qui fournit l'approximant désiré.

## Approximants de Padé-Hermite de type arbitraire

Pour calculer un approximant de Padé-Hermite de type $\mathbf{d} = (d_1, \dots, d_n)$, il suffit de remplacer, dans l'Algorithme 7.3, deg par $\deg_{\mathbf{d}}$ et type par $\text{type}_{\mathbf{d}}$, où :

$$ \deg_{\mathbf{d}}(\mathbf{P}) = \max_{1 \le i \le n} \{\deg(P_i) - d_i\}, \quad \text{type}_{\mathbf{d}}(\mathbf{P}) = \max\{i \mid \deg(P_i) - d_i = \deg_{\mathbf{d}}(\mathbf{P})\}. $$

## Applications

Si l'on sait qu'une série est rationnelle de degré au plus $d$, l'approximation de Padé de type $(d, d)$ suffit pour reconstruire la fraction rationnelle. Une question naturelle est de généraliser cette observation dans le cadre de l'approximation de Padé-Hermite. Les réponses sont partielles, mais entièrement satisfaisantes dans la pratique.

### Recherche de relations à coefficients polynomiaux entre séries formelles

Supposons qu'il existe une combinaison linéaire $\sum_{i=1}^n P_i f_i = 0$, à coefficients des polynômes $P_i(X) \in \mathbb{K}[X]$ de degrés au plus $d$. Par ailleurs, supposons calculé un approximant de Padé-Hermite $\mathbf{Q} = (\mathbf{Q}_1, \dots, \mathbf{Q}_n)$ de $\mathbf{F} = (f_1, \dots, f_n)$ de type $\mathbf{d} = (d, \dots, d)$, via l'algorithme de Derksen. La question est : quel lien y a-t-il entre $\mathbf{P} = (P_1, \dots, P_n)$ et $\mathbf{Q}$?

D'abord, dans le cas générique, la réponse est très simple : $\mathbf{P}$ et $\mathbf{Q}$ sont identiques, à un coefficient scalaire près. En effet, $\mathbf{Q} = \mathbf{Q}_n$ et $\mathbf{Q}_1, \dots, \mathbf{Q}_n$ est une base de $V_{nd+n-1}$ dont les degrés sont décrits page 150. En particulier, $\mathbf{P}$ doit être une combinaison linéaire des $\mathbf{Q}_i$. Des considérations sur les degrés, utilisant la forme bien particulière des degrés des $\mathbf{Q}_i$, mènent à la conclusion désirée.

En effet, si $(c_1(X), \dots, c_n(X)) \cdot {}^t(\mathbf{Q}_1, \dots, \mathbf{Q}_n) = \mathbf{P}$, alors $c_1 \mathbf{Q}_{11} + \dots + c_n \mathbf{Q}_{n1} = \mathbf{P}_1$, et comme $\deg(\mathbf{Q}_{11}) = d+1$ et $\deg(\mathbf{Q}_{j1}) = d$ pour $j > 1$ et $\deg(\mathbf{P}_1) \le d$, on obtient que $c_1 = 0$. De même, $c_2 = \dots = c_{n-1} = 0$ et $c_n$ doit être une constante $c \in \mathbb{K}$ telle que $\mathbf{P} = c\mathbf{Q}$.

Dans le cas général, un argument de noethérianité permet de prouver que pour $D \gg 0$, $V_{nD+n-1}$ contient la relation $\mathbf{P}$, qui sera trouvée par l'algorithme de Derksen. Seulement, on ne dispose pas de borne $a$ priori en fonction de $d$, sur le D minimal avec cette propriété. En effet, si on note

$$ W_j = \{\mathbf{Q} \in \mathbb{K}[X]^n \mid \text{val}(\mathbf{Q} \cdot \mathbf{F}) \ge j \text{ et } \deg(\mathbf{Q}) \le \deg(\mathbf{P})\}, $$

et $W_\infty = \bigcap_{j \ge 0} W_j$, alors $W_\infty$ contient toutes les relations de $\mathbf{F}$ en degré $d$, et en particulier $\mathbf{P}$. Puisque $W_0 \supseteq W_1 \supseteq W_2 \supseteq \dots$ est une suite décroissante d'espaces vectoriels de dimension finie, elle est stationnaire, donc il existe un N tel que $W_\infty = W_N = W_{N+1} = \dots$.
Le cas normal correspond à la situation où $\dim(W_{k+1}) = \dim(W_k) - 1$ pour chaque $k$ (noter la ressemblance avec la normalité de la suite des restes dans l’algorithme d’Euclide).

### Reconstruction d'équations algébriques et différentielles

Deux cas particuliers importants, pour lesquels on peut être encore plus précis, sont les approximants algébriques et différentiels.

Soit $A \in \mathbb{K}[[X]]$ une série formelle algébrique. Le problème d'approximation algébrique consiste à retrouver un polynôme $P(X, Y) \in \mathbb{K}[X, Y]$ tel que $P(X, A(X)) = 0$ à partir des premiers termes de $A$. Si un tel $P$, de degré $d$ en $X$ et $n$ en $Y$, existe, alors les coefficients des puissances de $Y$ formeront un approximant de Padé-Hermite de type $(d, \dots, d)$ du vecteur de séries $(1, A, \dots, A^n)$. La difficulté vient de ce que un calcul d'approximation de Padé-Hermite ne trouve, a priori, qu'un polynôme $Q \in \mathbb{K}[X, Y]$ tel que $Q(X, A(X)) = 0 \bmod X^\sigma$, où $\sigma = (n+1)d - 1$. On dit alors qu'on a deviné un polynôme annulateur $Q$ de $A$. Pour les approximants différentiels, la problématique est la même, la seule différence étant qu'on calcule un approximant de Padé-Hermite du vecteur des dérivées successives $(A, A', \dots, A^{(n)})$.

### Certification d'identités algébriques

Pour un approximant algébrique, il se pose la question de la certification a posteriori. Pour un polynôme annulateur deviné $Q$, cela veut dire qu'on souhaiterait déduire que non seulement $Q(X, A(X))$ est nul modulo $X^\sigma$, mais aussi que $Q(X, A(X)) = 0$.

Le résultat suivant apporte une réponse partielle à la question de la certification. Son avantage est qu'il ne dépend pas de l'algorithme utilisé pour produire l'approximant de Padé-Hermite.

**Théorème 7.15** Supposons que $A \in \mathbb{K}[[X]]$ est racine d'un polynôme irréductible de $\mathbb{K}[X, Y]$ de degré au plus $d$ en $X$ et au plus $n$ en $Y$. Soit

$$Q = (Q_0, Q_1, \dots, Q_n)$$

un approximant de Padé-Hermite de type $(d, \dots, d)$ de

$$F = (1, A, \dots, A^n).$$

Si $\text{val}(Q \cdot F) \ge 2dn$, alors $Q \cdot F = 0$, c'est-à-dire que $A$ est racine du polynôme $Q = \sum_{i=0}^n Q_i Y^i$.

*Démonstration.* Soit $P \in \mathbb{K}[X, Y]$ un polynôme irréductible de degré au plus $d$ en $X$ et au plus $n$ en $Y$ tel que $P(X, A) = 0$. Le polynôme $\text{Res}_Y(P, Q) \in \mathbb{K}[X]$ est de degré borné par $\text{deg}_X(P)\text{deg}_Y(Q) + \text{deg}_Y(P)\text{deg}_X(Q) \le 2dn$. Comme il est de la forme $uP + vQ$ pour $u$ et $v$ deux polynômes, avec $\text{deg}_Y(v) < \text{deg}_Y(P)$, et qu'il ne change pas lorsqu'on remplace la variable $Y$ par la série $A$, c'est un $O(X^{2dn})$. Par conséquent $uP + vQ = 0$, donc $P$ divise $vQ$ et, étant irréductible, il divise $v$ ou $Q$. Or $\text{deg}_Y(v) < \text{deg}_Y(P)$, donc finalement $Q \mid P$ et $Q(X, A)$ est la série nulle. ■
## Reconstruction de récurrences

Soit $(a_n)_n \in \mathbb{K}^\mathbb{N}$ une suite vérifiant une récurrence linéaire à coefficients polynomiaux. Comment, à partir des premiers termes de la suite, retrouver les coefficients de cette récurrence ?

Une possibilité consiste à utiliser la correspondance montrée au Théorème 14.1 en page 253, et à chercher une équation différentielle, à coefficients polynomiaux, portant sur la série génératrice de la suite. Cette équation peut être devinée comme expliqué précédemment, grâce à une approximation de Padé-Hermite. Il suffit ensuite de passer de l'équation différentielle à l'expression de la récurrence, ce qui n'est pas trivial mais purement formel et a été décrit au Chapitre 14.

## Calcul quasi-optimal

Pour le calcul d'approximants de Padé-Hermite, on dispose d'un algorithme quasi-optimal, dont l'efficacité est utilisée en calcul formel comme la base de nombreux autres algorithmes. Cet algorithme peut être vu comme une alternative à l'approche par *matrices structurées*, de même complexité, décrite au Chapitre 10.

**Théorème 7.16** Soient $\mathbf{F} = (f_1, \dots, f_n) \in \mathbb{K}[\mathbf{X}]^n$, $\mathbf{d} = (d_1, \dots, d_n) \in \mathbb{N}^n$ et $\sigma = d_1 + \dots + d_n + n - 1$. Il est possible de calculer un approximant de Padé-Hermite de type $\mathbf{d}$ de $\mathbf{F}$ en $O(\mathrm{MM}(n, \sigma) \log \sigma) = O(n^0 M(\sigma) \log \sigma)$ opérations dans $\mathbb{K}$.

*Esquisse de démonstration.* Comme pour l'algorithmes de complexité quadratique de Derksen, une première idée est de calculer non pas un vecteur, mais toute une matrice d'approximants. Une seconde idée est d'utiliser un « diviser pour régner » tout en exploitant la multiplication rapide de matrices polynomiales. Voici l'esquisse de l'algorithmes dans le cas où $\mathbf{d} = (d, \dots, d)$ :

1. Si $N = 1$, renvoyer une base de l'orthogonal à $\mathbf{F}(0)$.

2. Sinon :

(a) calculer récursivement une matrice polynomiale d'approximants $\mathbf{P}_1$ telle que $\mathbf{P}_1 \cdot \mathbf{F} = \mathrm{O}(X^{\sigma/2})$ et $\mathrm{deg}(\mathbf{P}_1)$ vaille environ $d/2$;

(b) calculer le reste $\mathbf{R}$ tel que $\mathbf{P}_1 \cdot \mathbf{F} = X^{\sigma/2} \cdot (\mathbf{R} + \mathrm{O}(X^{\sigma/2}))$;

(c) calculer récursivement une matrice polynomiale $\mathbf{P}_2$ d'approximants du reste : $\mathbf{P}_2 \cdot \mathbf{R} = \mathrm{O}(X^{\sigma/2})$, avec $\mathrm{deg}(\mathbf{P}_2)$ environ égal à $d/2$;

(d) renvoyer $\mathbf{P} = \mathbf{P}_2 \cdot \mathbf{P}_1$.

Le choix précis des degrés des approximants $\mathbf{P}_1$ et $\mathbf{P}_2$ est une question délicate et ne sera pas détaillé. En admettant ces choix, la correction de l'algorithmes et son analyse de complexité se déduisent aisément. ■

## Exercices

**Exercice 7.3 — Interpolation creuse.** Soit $\mathbb{K}$ un corps effectif de caractéristique zéro. Soit $\mathbf{F}$ une fonction d'arité $n$ sur $\mathbb{K}$ qui calcule les valeurs d'un polynôme $P$ de $\mathbb{K}[X_1, \dots, X_n]$ composé d'au plus $t$ termes non nuls. Le calcul de tous les termes non nuls de $P$ en évaluant seulement $\mathbf{F}$ sur plusieurs points bien choisis est appelé
interpolation creuse.

Soient $p_1, \dots, p_n$ des nombres premiers distincts (en pratique on choisit les plus petits).

1. Soient $X_1^{e_1} \cdots X_n^{e_n}$ et $X_1^{f_1} \cdots X_n^{f_n}$ deux monômes distincts. Expliquer brièvement pourquoi leurs valeurs en $(p_1, \dots, p_n)$ sont distinctes.

2. Supposons que P contienne exactement $\tau \le t$ termes non nuls, et qu'il s'écrive sous la forme $\sum_{i=1}^{\tau} P_i X_1^{e_{i,1}} \cdots X_n^{e_{i,n}}$. Montrer l'identité suivante :

$$ \sum_{i \ge 0} F(p_1^i, \dots, p_n^i) t^i = \sum_{i=1}^{\tau} \frac{P_i}{1 - p_1^{e_{i,1}} \cdots p_n^{e_{i,n}} t} $$

3. Déduire de la formule précédente un algorithme pour obtenir les coefficients du polynôme $D(t) = \prod_{i=1}^{\tau} (1 - p_1^{e_{i,1}} \cdots p_n^{e_{i,n}} t)$ de $\mathbb{Q}[t]$. Estimer le coût de l'algorithme en fonction du nombre d'appels à la fonction F et du nombre d'opérations arithmétiques dans $\mathbb{K}$.

4. En admettant que l'on dispose d'un algorithme pour calculer toutes les racines de D dans $\mathbb{Q}$, expliquer comment obtenir toutes les valeurs $e_{i,j}$, $1 \le i \le \tau$, $1 \le j \le n$.

5. En utilisant les algorithmes rapides du Chapitre 5, et les racines de D, proposer une façon efficace de calculer tous les coefficients $P_i$ de P.

■

## Notes

La présentation de la Section 7.1 est fortement inspirée de l'article de McEliece et Shearer [MS78], celle de la Section 7.2 par l'article de Derksen [Der94]. Les approximants de Padé sont des objets très classiques en analyse numérique; on pourra par exemple consulter le livre de Baker et Graves-Morris [BG96] qui leur est entièrement dédié.

L'algorithme de Berlekamp-Massey a été initialement conçu pour décoder des codes BCH, puis étendu par Zierler [Zie68] et Massey [Mas69], qui ont mis en évidence le lien avec les suites récurrentes linéaires à coefficients constants. Les liens forts entre l'algorithme d'Euclide, l'algorithme de Berlekamp-Massey, l'approximation de Padé et le calcul de fractions continues ont été étudiés par Mills [Mil75], Welch et Scholtz [WS79], Cheng [Che84], et Dornstetter [Dor87]. L'utilisation de ces méthodes pour l'interpolation creuse est à la base de l'algorithme dû à Ben-Or et Tiwari [BT88] adapté dans l'Exercice 7.3.

Les premiers algorithmes rapides pour la reconstruction rationnelle sont dus à Gustavson et Yun [GY79]. Ces algorithmes ont été améliorés par ces deux auteurs en collaboration avec Brent [BGY80] et Gragg-Warner [Gra+82], puis par Cabay et Choi [CC86]. La problématique de la reconstruction rationnelle pour les entiers est classique en théorie des nombres. En calcul formel, elle a été très étudiée en lien avec les algorithmes modulaires, par exemple pour la résolution de systèmes linéaires par remontées $p$-adiques par Dixon [Dix82]. Des algorithmes rapides ont été récemment conçus par Wang et Pan [PW04; WP03].

Les approximants de Padé et de Padé-Hermite ont été introduits et étudiés par
Hermite et Padé [Her73a; Her93; Pad92; Pad94]. Une exposition plus moderne se trouve dans un article de Mahler [Mah68]. Le problème d'approximation de Padé-Hermite a été utilisé par Hermite en 1873 dans sa preuve [Her73b] de la transcendance de *e*, qui utilise le choix très particulier $F = (1, e^X, e^{2X}, \dots)$. Deux autres cas particuliers importants sont : les « approximants algébriques » [Sha74; TD00] avec $F = (1, f, f^2, \dots)$ et les « approximants différentiels » [Kha02] avec $F = (f, f', f'', \dots)$, où $f$ est une série de $\mathbb{K}[[X]]$ donnée.

En calcul formel, les approximants de Padé-Hermite ont été introduits par Della Dora et Dicrescenzo [DD84a; DD84b]. Un premier algorithme de complexité quadratique a été donné par Sergeyev [Ser86]. Cet article est resté méconnu, et d'autres algorithmes quadratiques ont été proposés par Paszkowski [Pas87], Cabay et Labahn [CL89], Van Barel et Bultheel [VB91b], et Cabay, Labahn, Beckermann [CLB92]. L'algorithme de complexité quasi-optimale esquissé en Section 7.2 a été conçu par Beckermann et Labahn [BL94].

Il existe diverses généralisations utiles du problème d'approximation de Padé-Hermite, par exemple les « approximants de Padé-Hermite simultanés et matriciels », ou encore des approximants modulo un polynôme arbitraire de degré $\sigma = \sum_i (d_i + 1) - 1$ au lieu de $X^\sigma$. Des algorithmes rapides existent pour toutes ces généralisations [Bec92; BL00; BL94; GJV03; VB92].

## Bibliographie

Bec92 BECKERMANN, Bernhard (1992). « A reliable method for computing M-Padé approximants on arbitrary staircases ». In : *Journal of Computational and Applied Mathematics*, vol. 40, n°1, p. 19-42.

BG96 BAKER Jr., George A. et Peter GRAVES-MORRIS (1996). *Padé approximants*. 2ᵉ éd. Vol. 59. Encyclopedia of Mathematics and its Applications. Cambridge University Press.

BGY80 BRENT, Richard P., Fred G. GUSTAVSON et David Y. YUN (1980). « Fast solution of Toeplitz systems of equations and computation of Padé approximants ». In : *Journal of Algorithms*, vol. 1, n°3, p. 259-295.

BL00 BECKERMANN, Bernhard et George LABAHN (2000). « Fraction-free computation of matrix rational interpolants and matrix GCDs ». In : *SIAM Journal on Matrix Analysis and Applications*, vol. 22, n°1, p. 114-144.

BL94 BECKERMANN, B. et G. LABAHN (1994). « A uniform approach for the fast computation of matrix-type Padé approximants ». In : *SIAM Journal on Matrix Analysis and Applications*, vol. 15, n°3, p. 804-823.

BT88 BEN-OR, M. et P. TIWARI (1988). « A deterministic algorithm for sparse multivariate polynomial interpolation ». In : *STOC'88 : ACM Symposium on Theory of Computing*. ACM Press, p. 301-309.

CC86 CABAY, Stanley et Dong Koo CHOI (1986). « Algebraic computations of scaled Padé fractions ». In : *SIAM Journal on Computing*, vol. 15, n°1, p. 243-270.

Che84 CHENG, Unjeng (1984). « On the continued fraction and Berlekamp's algorithm ». In : *IEEE Transactions on Information Theory*, vol. 30, n°3, p. 541-544.
CL89 CABAY, S. et G. LABAHN (1989). « A fast, reliable algorithm for calculating Padé-Hermite forms ». In : *ISSAC'89 : International Symposium on Symbolic and Algebraic Computation*. Portland, Oregon, United States : ACM Press, p. 95-100.

CLB92 CABAY, S., G. LABAHN et B. BECKERMANN (1992). « On the theory and computation of nonperfect Padé-Hermite approximants ». In : *Journal of Computational and Applied Mathematics*, vol. 39, n°3, p. 295-313.

DD84a DELLA DORA, J. et C. DICRESCENZO (1984). « Approximants de Padé-Hermite. I. Théorie ». In : *Numerische Mathematik*, vol. 43, n°1, p. 23-39.

DD84b — (1984). « Approximants de Padé-Hermite. II. Programmation ». In : *Numerische Mathematik*, vol. 43, n°1, p. 41-57.

Der94 DERKSEN, Harm (1994). *An algorithm to compute generalized Padé-Hermite forms*. Report n°9403. Dept. of Math., Catholic University Nijmegen.

Dix82 DIXON, John D. (1982). « Exact solution of linear equations using p-adic expansions ». In : *Numerische Mathematik*, vol. 40, n°1, p. 137-141.

Dor87 DORNSTETTER, Jean-Louis (1987). « On the equivalence between Berlekamp's and Euclid's algorithms ». In : *IEEE Transactions on Information Theory*, vol. 33, n°3, p. 428-431.

GJV03 GIORGI, Pascal, Claude-Pierre JEANNEROD et Gilles VILLARD (2003). « On the complexity of polynomial matrix computations ». In : *ISSAC'03 : International Symposium on Symbolic and Algebraic Computation*. Éd. par J. R. SENDRA. New York : ACM Press, p. 135-142.

Gra+82 GRAGG, William B., Fred G. GUSTAVSON, Daniel D. WARNER et David Y. Y. YUN (1982). « On fast computation of superdiagonal Padé fractions ». In : *Mathematical Programming Study*, n°18. Algorithms and theory in filtering and control (Lexington, Ky., 1980), p. 39-42.

GY79 GUSTAVSON, Fred G. et David Y. Y. YUN (1979). « Fast algorithms for rational Hermite approximation and solution of Toeplitz systems ». In : *IEEE Transactions on Circuits and Systems*, vol. 26, n°9, p. 750-755.

Her73a HERMITE, C. (1873). « Extrait d'une lettre de Monsieur Ch. Hermite à Monsieur Paul Gordan ». In : *Journal für die Reine und Angewandte Mathematik*, vol. 78, p. 303-311.

Her73b — (1873). « Sur la fonction exponentielle ». In : *Comptes-rendus de l'Académie des sciences de Paris*, vol. 77, p. 18-24.

Her93 — (1893). « Sur la généralisation des fractions continues algébriques ». In : *Annali di Matematica pura ed applicata*, vol. 21, n°2, p. 289-308.

Kha02 KHAN, M. A. H. (2002). « High-order differential approximants ». In : *Journal of Computational and Applied Mathematics*, vol. 149, n°2, p. 457-468.

Mah68 MAHLER, K. (1968). « Perfect systems ». In : *Compositio Mathematica*, vol. 19, 95-166 (1968).

Mas69 MASSEY, James L. (1969). « Shift-register synthesis and BCH decoding ». In : *IEEE Transactions on Information Theory*, vol. IT-15, p. 122-127.
Mil75 MILLS, W. H. (1975). « Continued fractions and linear recurrences ». In : Mathematics of Computation, vol. 29, p. 173–180.

MS78 MCELIECE, Robert J. et James B. SHEARER (1978). « A property of Euclid's algorithm and an application to Padé approximation ». In : SIAM Journal on Applied Mathematics, vol. 34, n°4, p. 611–615.

Pad92 PADÉ, H. (1892). « Sur la représentation approchée d'une fonction par des fractions rationnelles ». In : Annales scientifiques de l'É.N.S. vol. 9, p. 3–93.

Pad94 — (1894). « Sur la généralisation des fractions continues algé-briques ». In : Journal de mathématiques pures et appliquées, vol. 10, n°4, p. 291–330.

Pas87 PASZKOWSKI, Stefan (1987). « Recurrence relations in Padé–Hermite ap-proximation ». In : Journal of Computational and Applied Mathematics, vol. 19, n°1, p. 99–107.

PW04 PAN, Victor Y. et Xinmao WANG (2004). « On rational number reconstruc-tion and approximation ». In : SIAM Journal on Computing, vol. 33, n°2, p. 502–503.

Ser86 SERGEYEV, A. V. (1986). « A recursive algorithm for Padé–Hermite approxi-mations ». In : USSR Comput. Maths Math. Phys, vol. 26, n°2, p. 17–22.

Sha74 SHAFER, R. E. (1974). « On quadratic approximation ». In : SIAM Journal on Numerical Analysis, vol. 11, p. 447–460.

TD00 TOURIGNY, Y. et Ph. G. DRAZIN (2000). « The asymptotic behaviour of algebraic approximants ». In : The Royal Society of London. Proceedings. Series A. Mathematical, Physical and Engineering Sciences, vol. 456, n°1997, p. 1117–1137.

VB91b VAN BAREL, Marc et Adhemar BULTHEEL (1991). « The computation of nonperfect Padé–Hermite approximants ». In : Numerical Algorithms, vol. 1, n°3, p. 285–304.

VB92 VAN BAREL, M. et A. BULTHEEL (1992). « A general module-theoretic frame-work for vector M-Padé and matrix rational interpolation ». In : Numerical Algorithms, vol. 3, n°1-4, p. 451–461.

WP03 WANG, Xinmao et Victor Y. PAN (2003). « Acceleration of Euclidean algo-rithm and rational number reconstruction ». In : SIAM Journal on Computing, vol. 32, n°2, p. 548–556.

WS79 WELCH, L. R. et R. A. SCHOLTZ (1979). « Continued fractions and Berle-kamp's algorithm ». In : IEEE Transactions on Information Theory, vol. 25, n°1, p. 19–27.

Zie68 ZIERLER, N. (1968). « Linear recurring sequences and error-correcting codes ». In : Error Correcting Codes (Proc. Sympos. Math. Res. Center, Madison, Wis., 1968). Wiley, p. 47–59.
# Matrices

8 Algèbre linéaire dense : de Gauss à Strassen 159

8.1 Introduction

8.2 Multiplication de matrices

8.3 Autres problèmes d'algèbre linéaire

9 Algèbre linéaire creuse : algorithme de Wiedemann 189

9.1 Introduction

9.2 Polynôme minimal et résolution de systèmes

9.3 Calcul du polynôme minimal

9.4 Calcul du déterminant

9.5 Calcul du rang

10 Algèbre linéaire structurée 195

10.1 Introduction

10.2 Le cas quasi-Toeplitz

11 Systèmes linéaires à coefficients polynomiaux 209

11.1 Des séries aux solutions rationnelles

11.2 Développement comme une fraction rationnelle

11.3 L'algorithme de Storjohann

12 Principe de transposition de Tellegen 217

12.1 Introduction

12.2 La version en termes de graphes du principe de Tellegen

12.3 Principe de Tellegen pour les programmes linéaires

12.4 Applications

13 Équations différentielles à coefficients séries 237

13.1 Équations différentielles d'ordre 1

13.2 Équations d'ordre supérieur et systèmes d'ordre 1

13.3 Cas particuliers

13.4 Extensions
<figure><img src="image_27.png" /></figure>
# 8. Algèbre linéaire dense : de Gauss à Strassen

## Résumé

L'algèbre linéaire est au cœur de nombreux problèmes algorithmiques. La multiplication matricielle usuelle et la méthode du pivot de Gauss ont une complexité arithmétique cubique en la taille de l'entrée. Dans ce chapitre, nous montrons que pour la plupart des questions d'algèbre linéaire dense — multiplication, inversion, déterminant, résolution de système — il existe des algorithmes plus efficaces, de complexité strictement sous-cubique.

## 8.1 Introduction

En mathématiques, il est de coutume de considérer qu'un problème est rendu tri- vial lorsqu'il est ramené à une question d'algèbre linéaire. Du point de vue calculatoire se pose cependant la question de l'efficacité des opérations matricielles.

### Classification

Les problèmes algorithmiques en algèbre linéaire cachent des difficultés très subtiles. En général, les premières questions que l'on est amené à se poser en rapport avec la manipulation des matrices concernent le calcul efficace du produit matrice-matrice, du produit matrice-vecteur, de l'inverse, ainsi que la résolution de systèmes.

Les réponses à ces questions varient fortement selon le type de matrices considé- rées. Une classification possible est la suivante.

### Les matrices denses

Ce sont les matrices quelconques, sans aucune structure particulière. Nous verrons que l'algorithme des matrices denses peut se ramener essentiellement au produit matriciel, dont la complexité est une question extrêmement délicate.
## Les matrices creuses

Beaucoup de problèmes linéaires se formulent en termes de matrices possédant un nombre important d'éléments nuls, dites « creuses ». Dans ce cas, l'algorithmique dense est inappropriée; il est possible de tirer parti de la forme creuse en utilisant des outils mieux adaptés, à base de récurrences linéaires.

## Les matrices structurées

Enfin, il existe des familles de matrices particulières, souvent associées à des applications linéaires entre espaces de polynômes, telles que les matrices de type Sylvester pour le résultant, de type Vandermonde pour l'évaluation et l'interpolation, de type Toeplitz et Hankel pour l'approximation de Padé et de Padé-Hermite, etc. Pour ces types de matrices clairement identifiés, on développe soit une algorithmique ad hoc, soit une théorie unifiée reposant sur le concept de « rang de déplacement ».

Schématiquement, l'algorithmique des matrices denses est la plus lente, de complexité située entre $O(n^2)$ et $O(n^3)$, en taille $n$. La complexité du calcul avec les matrices creuses est de l'ordre de $O(n^2)$, alors que celle des matrices structurées est en $O(n)$, à des facteurs logarithmiques près.

L'algorithmique des matrices denses constitue l'objet d'étude de ce chapitre. Les matrices creuses seront abordées brièvement au Chapitre 9 et les matrices structurées feront l'objet d'une étude plus approfondie au Chapitre 10.

Nous allons travailler avec un corps effectif noté $K$ et avec l'algèbre des matrices carrées $M_n(K)$ à coefficients dans $K$. Signalons toutefois que la plupart des résultats de ce chapitre s'étendent au cas où le corps $K$ est remplacé par un anneau effectif $A$, et les matrices carrées par des matrices rectangulaires à coefficients dans $A$.

## Résultat principal

Les questions qui seront étudiées, ou simplement évoquées, dans ce chapitre sont celles de la complexité du calcul du produit matriciel dans $M_n(K)$, de l'inverse, du déterminant, du polynôme caractéristique ou minimal, de la résolution de systèmes linéaires, ou encore de la mise sous diverses « formes canoniques » (échelon, LUP, compagnon par blocs, ...).

Pour ces questions, on dispose d'une première famille d'algorithmes « naïfs », reposant sur l'application directe des définitions ou des propriétés mathématiques. Si l'algorithme naïf de multiplication a une complexité raisonnable, cubique en la taille, dans d'autres cas l'emploi des algorithmes naïfs est vivement déconseillé. Par exemple, si $A = (a_{i,j})_{i,j=1}^n$ est une matrice de $M_n(K)$, on n'utilisera pas la définition

$$ \det(A) = \sum_{\sigma \in S_n} \operatorname{sgn}(\sigma) \prod_{i=1}^{n} a_{i,\sigma(i)} $$

pour le calcul de son déterminant; en effet, cela mènerait à un algorithme de complexité $O(n \cdot n!)$, donc hyper-exponentielle en la taille de $A$. De même, la définition du rang d'une matrice comme la taille maximale d'un mineur non nul ne se prête pas directement, via une recherche exhaustive, à une algorithmisation efficace. Dans un registre légèrement différent, les formules de Cramer pour résoudre un système
linéaire ne sont pas d'une grande utilité pratique (en revanche, elles servent à borner les tailles des solutions, et s'avèrent souvent utiles dans les estimations de complexité).

Une seconde classe d'algorithmes, basés sur la *méthode d'élimination de Gauss*, permet de résoudre de manière raisonnablement efficace la plupart des problèmes évoqués plus haut. Appliquée à une matrice A de $M_n(\mathbb{K})$, cette méthode fournit des algorithmes de complexité $O(n^3)$ pour le calcul du rang $rg(A)$, du déterminant $det(A)$, de l'inverse $A^{-1}$ (si A est inversible), pour le calcul d'une forme échelonnée et d'une décomposition LUP de A, ainsi que pour la résolution d'un système linéaire de matrice A.

Le résultat principal de ce chapitre est que l'on peut faire mieux, et qu'il existe une constante $2 \le \omega < 3$ (dépendante *a priori* du corps de base $\mathbb{K}$), qui contrôle la complexité du produit matriciel et de toutes les opérations d'algèbre linéaire.

Pour formaliser ce point, on associe à chacun des problèmes son *exposant*. Dans le cas du produit, il est défini comme suit.

**Définition 8.1 — Exposant de la multiplication matricielle.** On dit que le réel $\theta$ est un exposant réalisable pour la multiplication de matrices sur le corps $\mathbb{K}$ s'il existe un algorithme de complexité arithmétique $O(n^\theta)$ pour multiplier deux matrices arbitraires de $M_n(\mathbb{K})$. On définit $\omega_{\text{mul}} = \inf\{\theta \mid \text{\textbf{\theta}} \text{ est un exposant réalisable}\}$.

*A priori*, on devrait mettre en indice le corps $\mathbb{K}$ pour les exposants $\theta$ et $\omega_{\text{mul}}$. Les résultats présentés dans la suite sont cependant indépendants de $\mathbb{K}^1$. L'algorithme naïf pour le produit est de complexité $O(n^3)$, donc $\theta = 3$ est un exposant réalisable et $\omega_{\text{mul}} \le 3$. D'un autre côté, $\omega_{\text{mul}}$ est forcément au moins égal à 2 : le nombre d'éléments d'une matrice est $n^2$ et il faut bien autant d'opérations pour l'écrire.

Ensuite, on peut définir de la même manière les exposants de tous les autres problèmes, de la forme $\omega_{\text{inv}}$ pour l'inversion, $\omega_{\text{polcar}}$ pour le polynôme caractéristique, $\omega_{\text{det}}$ pour le déterminant, $\omega_{\text{LUP}}$ pour la décomposition LUP, etc. On a alors le résultat suivant.

**Théorème 8.1** Les exposants $\omega_{\text{mul}}$, $\omega_{\text{inv}}$, $\omega_{\text{polcar}}$, $\omega_{\text{det}}$, $\omega_{\text{LUP}}$, sont tous égaux, et inférieurs à 2,38; on note $\omega$ leur valeur commune. L'ex exposant correspondant à la résolution de systèmes linéaires est inférieur ou égal à $\omega$.

Ce théorème contient deux (familles de) résultats :

- des preuves d'équivalence entre problèmes,

- des bornes supérieures sur leur complexité, c'est-à-dire des algorithmes.

Dans le présent chapitre, il est impossible de démontrer ce théorème dans son intégralité. Nous nous contenterons de prouver (et encore, sous des hypothèses simplificatrices) que le produit et l'inverse ont le même exposant, et que celui-ci est inférieur à $\log_2 7 < 2,81$.

Nous montrerons également que $\omega_{\text{det}} \le \omega_{\text{polcar}} \le \omega_{\text{mul}} = \omega_{\text{inv}}$ et laisserons en exercice (Exercice 8.13) la preuve de l'inégalité $\omega_{\text{inv}} \le \omega_{\text{det}}$ qui permet de conclure que cette suite d'inégalités est en fait une suite d'égalités.

<footnote n="1">1. On peut montrer que la constante $\omega_{\text{mul}}$ ne dépend que de la caractéristique du corps de base $\mathbb{K}$ et il est conjecturé qu'elle est même entièrement indépendante du corps $\mathbb{K}$.</footnote>
## Applications

Les problèmes d'algèbre linéaire abordés ci-dessus sont très fréquemment renforcés en pratique ; les résultats du Théorème 8.1 seront invoqués à plusieurs reprises dans cet ouvrage. Par exemple, l'algorithme de Beckermann-Labahn (Chapitre 7) pour le calcul d'approximants de Padé-Hermite repose ultimement sur des produits de matrices de polynômes. C'est aussi le cas de l'algorithme de Storjohann pour la résolution de systèmes linéaires à coefficients polynomiaux (Chapitre 11), et de l'algorithme de recherche de solutions séries d'équations algébriques (Chapitre 3) ou différentielles linéaires (Chapitre 13).

De nombreux algorithmes de résolution de systèmes polynomiaux reposent aussi sur de l'algèbre linéaire : les méthodes de calcul d'une base de Gröbner (Chapitre 23) peuvent se ramener à des systèmes linéaires (éventuellement creux) en très grande taille ; les résultats de complexité que nous donnerons pour les calculs de bases de Gröbner (Chapitre 26) reposent sur le calcul rapide de forme échelon ; l'algorithme de « résolution géométrique » (Chapitre 28) utilise une itération de Newton faisant intervenir des calculs d'inverses de matrices de séries formelles. Dans un contexte différent, les algorithmes de factorisation des polynômes (Partie IV), et ceux pour les séries différentiellement finies (Chapitre 14) et les algorithmes de sommation et d'intégration (Chapitres 29 et 32) font intervenir comme sous-problème la résolution d'un système linéaire. L'algèbre linéaire est au cœur de nombreux autres problèmes algorithmiques qui ne seront pas abordés dans cet ouvrage ; c'est le cas de la factorisation des entiers, ou encore du logarithme discret en cryptanalyse. Mentionnons enfin qu'une large partie des ordinateurs du monde passent leurs cycles à faire des calculs d'algèbre linéaire, que ce soit pour faire des calculs d'optimisation combinatoire (programmation linéaire par l'algorithme du simplexe), de la simulation numérique (méthode des éléments finis), ou de la recherche de pages web (le système de classement PageRank utilisé par le moteur de recherche Google repose sur la recherche d'un vecteur propre d'une gigantesque matrice creuse).

## 8.2 Multiplication de matrices

Tout comme le produit d'entiers et de polynômes, la multiplication matricielle est une opération fondamentale en calcul formel, dont l'étude de complexité s'avère être un problème hautement non trivial.

L'algorithme naïf pour effectuer le produit d'une matrice de taille $m \times n$ par un vecteur de taille $n$ utilise $O(mn)$ opérations arithmétiques : $mn$ multiplications et $m(n-1)$ additions. Un résultat de Winograd affirme qu'en général on ne peut pas faire mieux : par exemple, pour les matrices réelles dont les coefficients sont algébriquement indépendants sur $\mathbb{Q}$, l'algorithme naïf est optimal. Puisque le produit matriciel peut s'interpréter comme une succession de produits matrice-vecteur, une question naturelle est de savoir si la multiplication naïve de deux matrices est aussi quasi-optimale.

La réponse à cette question, et les résultats les plus importants de cette section, classés par ordre chronologique de découverte, sont regroupés dans le théorème suivant.