**Entrée** A dans $M_{r \times r}(\mathbb{A})$, $v \in M_{r \times 1}(\mathbb{A})$.

**Sortie** $Y = \sum_{i=0}^{N-1} Y_i X^i$ dans $M_{r \times r}(\mathbb{A})[X]$ tel que
$Y' = AY + O(X^{N-1})$.

1. Calculer les vecteurs $v, Av, A^2v, A^3v, \dots, A^{2r}v$.

2. Pour tout $j = 1, \dots, r$ faire :

a. reconstruire la fraction rationnelle $z_j$ ayant pour développement en série $\sum (A^i v)_j X^i$;

b. développer cette fraction en série à précision $X^N$;

c. reconstruire le développement de $y_j$ à partir de celui de $z_j$.

3. Renvoyer le vecteur $(y_1, \dots, y_r)$.

Algorithme 13.4 – Résolution approchée de $Y' = AY$.

## Équations différentielles linéaires à coefficients constants

Dans le cas où les coefficients de l'équation ou du système sont constants, la formule (13.1) s'applique encore et devient :

$$y(X) = \exp(XA)y(0).$$

Pour étudier les propriétés de cette solution, il est usuel de faire intervenir la forme de Jordan de la matrice, mais cette forme ne se prête pas facilement au calcul.

Un algorithme efficace est obtenu en exploitant la *transformée de Laplace formelle*.
Si $S = s_0 + s_1 X + s_2 X^2 + \dots$ est une série, cette transformée est définie par

$$\mathcal{L}S = s_0 + 1!s_1X + 2!s_2X^2 + \dots$$

Les troncatures de $\mathcal{L}S + O(X^N)$ et de la transformée inverse $\mathcal{L}^{-1}S + O(X^N)$ se calculent en N opérations. Cette transformation simplifie les systèmes différentiels linéaires à coefficients constants en les rendant linéaires non différentiels. Ceci découle de

$$\mathcal{L}(y) = \mathcal{L}\left((I + XA + \frac{1}{2!}X^2A^2 + \dots) y(0)\right) = \left(I + XA + X^2A^2 + \dots\right) y(0) = (I - XA)^{-1}y(0).$$

Les vecteurs solutions ont donc des transformées de Laplace dont les coordonnées sont rationnelles. En outre, d'après les formules de Cramer, le numérateur et le dénominateur de ces fractions ont des degrés bornés par l'ordre $r$ du système. L'Algorithme 13.4 s'en déduit. L'étape (1) est effectuée par la méthode naïve en $O(r^3)$ opérations. Une méthode plus efficace, en $O(\mathrm{MM}(r) \log r)$ opérations, est à la base de l'algorithme de Keller-Gehrig pour le calcul du polynôme caractéristique d'une matrice présenté au Chapitre 8. L'étape (2.a) se ramène à un calcul d'algèbre linéaire en posant des coefficients indéterminés pour les numérateurs et dénominateurs. Là aussi, le coût peut être diminué en faisant appel au calcul d'approximants de Padé (Chapitre 7), mais cette partie du calcul ne dépend pas de N. L'étape (2.b) utilise l'algorithme de développement de fractions rationnelles du Chapitre 4 (Théorème 4.9, page 88) en $O(NM(r)/r)$ opérations, et c'est là que se concentre le gain en complexité. Enfin, l'étape (2.c) ne demande que $O(N)$ opérations. L'ensemble de l'algorithme a donc une complexité arithmétique en $O(NM(r))$.
## 13.4 Extensions

### Composer se ramène à résoudre

Le seul algorithme du Chapitre 3 dont la complexité n'est pas quasi-optimale est l'algorithme de composition.

Dans les applications où l'on connaît une équation $\Phi(X, f, f', \dots, f^{(m)}) = 0$ vérifiée par la série $f$ il est parfois possible de calculer la série $h = f \circ g$ efficacement sans passer par l'algorithme de composition en remarquant qu'elle vérifie une équation du même ordre, obtenue en remplaçant $X$ par $g$ dans $\Phi$ et en composant les dérivées. Si $\Phi$ ne fait intervenir que des polynômes ou des fractions rationnelles, le résultat de cette opération a pour coefficients des séries, ce qui mène assez généralement à une complexité en $O(M(N))$ opérations en utilisant les algorithmes de ce chapitre.

**Exemple 13.1** Pour calculer le développement de $h = \tan(f)$, il est possible de calculer en $O(M(N))$ ceux de $\sin(f)$ et $\cos(f)$ (via l'exponentielle) et de diviser, mais il est aussi possible d'observer que $h$ vérifie l'équation $h' = 1 + h^2 f'$, et d'utiliser les méthodes ci-dessus. Il faut ensuite comparer les constantes dans les $O(\cdot)$ (ou deux implantations !) pour déterminer laquelle des deux méthodes est préférable.

**Exercice 13.3** Le $n^e$ polynôme de Legendre $P_n(X)$ est défini par l'équation $P_n(X) = \frac{d^n}{dX^n} (X^2 - 1)^n / (2^n n!)$. Montrer qu'il vérifie l'équation différentielle

$$ \frac{d}{dX} \left[ (1 - X^2)^n \frac{d}{dX} P_n(X) \right] + n(n+1)P_n(X) = 0. $$

En déduire que pour toute série de terme constant nul $F \in \mathbb{Q}[[X]]$, et pour tout $N \ge 0$, la série composée $P_N \circ F \bmod X^N$ peut se calculer en $O(M(N))$ opérations. ■

### Systèmes non linéaires

La linéarisation effectuée dans le cas des équations d'ordre 1 se généralise aux systèmes. L'énoncé devient alors le suivant.

**Théorème 13.5** Soient $\phi_1, \dots, \phi_r$ $r$ séries de $\mathbb{A}[[X, Y_1, \dots, Y_r]]$, telles que pour tout $r$-uplet de séries $(s_1(X), \dots, s_r(X)) \in \mathbb{A}[[X]]^r$, les $N$ premiers termes des $\phi_i(X, s_1, \dots, s_r)$ et de Jac$(\phi)(X, s_1(X), \dots, s_r(X))$ puissent être calculés en utilisant $O(L(N))$ opérations dans $\mathbb{A}$. On suppose en outre que $L(n)/n$ est une suite croissante. Alors, le système différentiel

$$ \left\{ \begin{array}{l} y'_1(t) = \phi_1(t, y_1(t), \dots, y_r(t)), \\ \vdots \\ y'_r(t) = \phi_r(t, y_1(t), \dots, y_r(t)), \end{array} \right. $$

avec conditions initiales $(y_1, \dots, y_r)(0) = v$ admet une solution série formelle, et ses $N$ premiers termes peuvent être calculés en

$$ O(L(N) + \min(MM(r, N), r^2 M(N) \log N)). $$
Le symbole Jac($\phi$) désigne la matrice jacobienne : son coefficient sur la ligne $i$ et la colonne $j$ vaut $\partial\phi_i/\partial y_j$. L'intérêt d'énoncer le théorème à l'aide de la fonction L est le même que dans le cas des équations.

*Démonstration.* Il s'agit essentiellement d'une généralisation de la Proposition 13.1. La linéarisation de $\Phi : Y \mapsto Y' - \varphi(Y)$ au voisinage d'une solution approchée $y$ s'obtient en écrivant :

$$ \Phi(y+h) = \Phi(y) + h' + Jac(\phi)(y)h + O(h^2). $$

L'équation à résoudre pour une itération de Newton est donc le système linéaire inhomogène

$$ h' = Jac(\phi)(y)h - (y' - \phi(y)). $$

Si $h$ est un vecteur solution de ce système à précision $2m$ et $y$ une solution de $\Phi$ à précision $m$, ces équations montrent que $y+h$ est solution à précision $2m$. Le cas non linéaire est ainsi ramené au cas linéaire. ■

## Notes

Les résultats de la Section 13.1 sont tirés de l'article de Brent et Kung de 1978 [BK78]. Une bonne partie des autres résultats de ce chapitre provient de l'article de Bostan, Chyzak, Ollivier, Salvy, Schost, et Sedoglavic de 2007 [Bos+07a]. La technique de la Section 13.2 dans le cas particulier de l'exponentielle d'une série est due à Hanrot, Quercia, et Zimmermann [HQZ04]. L'algorithme « diviser pour régner » de la Section 13.2 se trouve au moins implicitement dans des travaux de van der Hoeven [Hoe02c]. Des algorithmes plus récents peuvent être consultés dans l'article de van der Hoeven de 2010 [Hoe10].

## Bibliographie

* BK78 BRENT, R. P. et H. T. KUNG (1978). « Fast algorithms for manipulating formal power series ». In : *Journal of the ACM*, vol. 25, n°4, p. 581-595.

* Bos+07a BOSTAN, A., F. CHYZAK, F. OLLIVIER, B. SALVY, É. SCHOST et A. SEDOGLAVIC (2007). « Fast computation of power series solutions of systems of differential equations ». In : *SODA'07 : ACM-SIAM Symposium on Discrete Algorithms*. SIAM, p. 1012-1021.

* Hoe02c HOEVEN, Joris van der (2002). « Relax, but don't be too lazy ». In : *Journal of Symbolic Computation*, vol. 34, n°6, p. 479-542.

* Hoe10 — (2010). « Newton's method and FFT trading ». In : *Journal of Symbolic Computation*, vol. 45, n°8, p. 857-878.

* HQZ04 HANROT, Guillaume, Michel QUERCIA et Paul ZIMMERMANN (2004). « The middle product algorithm I ». In : *Applicable Algebra in Engineering, Communication and Computing*, vol. 14, n°6, p. 415-438.
<figure><img src="image_46.png" /></figure>
# Équations différentielles et récurrences linéaires

<table><tbody><tr><td>14</td><td>Séries différentiellement finies</td><td>251</td></tr><tr><td>14.1</td><td>Équations différentielles et récurrences</td><td></td></tr><tr><td>14.2</td><td>Propriétés de clôture</td><td></td></tr><tr><td>14.3</td><td>Séries algébriques</td><td></td></tr><tr><td>14.4</td><td>Au-delà</td><td></td></tr><tr><td>15</td><td>Récurrences linéaires à coefficients polynomiaux</td><td>269</td></tr><tr><td>15.1</td><td>Calcul naïf de N! et de suites polynomialement récurrentes</td><td></td></tr><tr><td>15.2</td><td>Pas de bébés, pas de géants</td><td></td></tr><tr><td>15.3</td><td>Scindage binaire</td><td></td></tr><tr><td>15.4</td><td>Récurrences singulières</td><td></td></tr><tr><td>16</td><td>Résolution de récurrences linéaires</td><td>289</td></tr><tr><td>16.1</td><td>Suites solution</td><td></td></tr><tr><td>16.2</td><td>Solutions à support fini</td><td></td></tr><tr><td>16.3</td><td>Solutions polynomiales</td><td></td></tr><tr><td>16.4</td><td>Solutions rationnelles</td><td></td></tr><tr><td>17</td><td>Résolution d'équations différentielles linéaires</td><td>311</td></tr><tr><td>17.1</td><td>Séries solution</td><td></td></tr><tr><td>17.2</td><td>Solutions polynomiales</td><td></td></tr><tr><td>17.3</td><td>Solutions rationnelles</td><td></td></tr><tr><td>17.4</td><td>Séries différentiellement finies analytiques</td><td></td></tr></tbody></table>
<figure><img src="image_47.png" /></figure>
# 14. Séries différentiellement finies

## Résumé

Les séries différentiellement finies sont les solutions d'équations différentielles linéaires à coefficients polynomiaux. Leurs coefficients vérifient des récurrences linéaires qui permettent de les manipuler plus rapidement que des séries arbitraires. L'équation différentielle linéaire ou la récurrence linéaire devient une structure de données adaptée pour manipuler les solutions. De nombreuses opérations dont l'addition et la multiplication peuvent être effectuées directement sur les équations. Une conséquence immédiate est une multiplication de ces séries différentiellement finies en temps linéaire.

Les séries différentiellement finies (c'est-à-dire solutions d'équations différentielles linéaires à coefficients polynomiaux) ont des coefficients qui satisfont une récurrence linéaire, ce qui permet d'en calculer les N premiers en $O(N)$ opérations, donc plus vite que la plupart des autres séries. Il est de ce fait crucial de reconnaître les séries qui sont différentiellement finies et de disposer des équations différentielles les définissant. De plus, les coefficients des séries différentiellement finies forment des suites qui sont appelées polynomialement récurrentes, dont l'algorithmique est évidemment étroitement liée à celle des séries différentiellement finies. L'importance de ces séries et suites provient en grande partie de leur omniprésence dans les applications. Ainsi, le *Handbook of Mathematical Functions*, référence importante en physique, chimie et mathématiques appliquées, comporte environ 60% de fonctions solutions d'équations différentielles linéaires ; de même, les suites polynomialement récurrentes forment environ un quart des plus de 100 000 suites référencées dans la version en ligne de l'*Encyclopedia of Integer Sequences* de N. Sloane. Enfin, l'algorithmique de ces séries permet la preuve automatique d'identités, présentée dans ce chapitre en présence d'une seule variable. La généralisation à plusieurs variables est traitée dans les chapitres de la cinquième partie.
## 14.1 Équations différentielles et récurrences

### Définitions

Si $\mathbb{K}$ désigne un corps, l'anneau des séries formelles à coefficients dans $\mathbb{K}$ est noté $\mathbb{K}[[X]]$. Ses principales propriétés ont été présentées au Chapitre 3. Son corps des fractions, noté $\mathbb{K}((X))$, est égal à $\mathbb{K}[[X]][1/X]$. Ses éléments sont appelés des séries de Laurent formelles. C'est une algèbre non seulement sur $\mathbb{K}$, mais aussi sur le corps des fractions rationnelles $\mathbb{K}(X)$.

**Définition 14.1** Une série formelle $A(X)$ à coefficients dans un corps $\mathbb{K}$ est dite *différentiellement finie* lorsque ses dérivées successives $A, A', \dots, \text{engendrent dans l'algèbre } \mathbb{K}((X)) \text{ des séries de Laurent formelles un espace vectoriel de dimension finie sur le corps } \mathbb{K}(X) \text{ des fractions rationnelles.}$

De manière équivalente, cette série est solution d'une équation différentielle linéaire à coefficients dans $\mathbb{K}(X)$ : si c'est le cas, alors l'équation différentielle permet de récrire toute dérivée d'ordre supérieur à celui de l'équation en termes des dérivées d'ordre moindre (en nombre borné par l'ordre). À l'inverse, si l'espace est de dimension finie, alors pour $m$ suffisamment grand, $A, A', \dots, A^{(m)}$ sont liées, et une relation de liaison entre ces dérivées est une équation différentielle linéaire.

**Exemple 14.1** De nombreuses séries classiques sont de ce type : $\exp(X)$, $\log(1+X)$, $(1+X)^a$, $\sin X$, $\sinh X$, $\cos X$, $\cosh X$, $\arcsin X$, $\arccos X$, $\arctan X$, ..., les séries hypergéométriques généralisées, et d'autres fonctions classiques de la physique mathématique.

R La Définition 14.1 pourrait être donnée plus généralement pour des séries de Laurent formelles de $\mathbb{K}((X))$, mais est sciemment restreinte ici aux séries formelles de $\mathbb{K}[[X]]$, dans l'objectif notamment de bien fonctionner avec la définition qui suit et avec le Théorème 14.1. Dans la Partie VI, la notion sera étendue à toute sorte d'espaces vectoriels sur $\mathbb{K}(X)$ clos par dérivation.

**Définition 14.2** Une suite $(a_n)_{n \ge 0}$ d'éléments d'un corps $\mathbb{K}$ est appelée suite *poly-nomialement récurrente* si elle satisfait à une récurrence de la forme

$$p_d(n)a_{n+d} + p_{d-1}(n)a_{n+d-1} + \dots + p_0(n)a_n = 0, \quad n \ge 0, \quad (14.1)$$

où les $p_i$ sont des polynômes de $\mathbb{K}[X]$.

Dans la suite, $\mathbb{K}$ aura toujours caractéristique nulle. On peut donc penser sans rien perdre aux idées à $\mathbb{K} = \mathbb{Q}$.

R Comme pour les séries différentiellement finies, on pourrait donner une définition en termes d'espaces vectoriels engendrés par les suites décalées, mais il faut alors faire attention aux zéros entiers du coefficient de tête de la récurrence $p_d$. Une possibilité est de considérer les *germes* de suites à l'infini, c'est-à-dire les classes d'équivalences de suites pour la relation « être égal à partir d'un certain
rang ». Cette approche n'est pas très effective, et nous travaillerons plutôt dans les chapitres suivants avec des conditions initiales généralisées.

**Exemple 14.2** Outre les solutions de récurrences linéaires à coefficients constants comme les nombres de Fibonacci, des suites classiques de ce type sont : $n!$, $1/n!$, $\binom{2n}{n}$, les nombres harmoniques $H_n$, les nombres de Catalan, de Motzkin,...

**Exemple 14.3** Le cas particulier des récurrences d'ordre 1 ($d = 1$) donne lieu aux suites hypergéométriques, qui jouent un rôle important dans les algorithmes de sommation symbolique, abordés au Chapitre 29.

## Équivalence

Le résultat qu'il s'agit de considérer d'un point de vue algorithmique est le suivant. Il s'agit d'un analogue du Lemme 4.7.

**Théorème 14.1** Une série formelle est différentiellement finie si et seulement si la suite de ses coefficients est polynomialement récurrente.

*Démonstration.* Soit $A(X) = a_0 + a_1 X + \cdots$ une série différentiellement finie et

$$q_m(X)A^{(m)}(X) + \cdots + q_0(X)A(X) = 0 \quad (14.2)$$

une équation différentielle qui l'annule. En notant $[X^n]f(X)$ le coefficient de $X^n$ dans la série $f(X)$ avec la convention que ce coefficient est nul pour $n < 0$, on a pour $n \ge 0$

$$[X^n]f'(X) = (n+1)[X^{n+1}]f(X), \quad [X^n]X^k f(X) = [X^{n-k}]f(X). \quad (14.3)$$

Par conséquent, l'extraction du coefficient de $X^n$ de (14.2) fournit une récurrence linéaire sur les $a_n$ valide pour tout $n \ge 0$ avec la convention $a_k = 0$ pour $k < 0$. Pour obtenir une récurrence de la forme (14.1), il faut décaler les indices de $n_0 := \max_{0 \le i \le m}(\deg q_i - i)$ si cet entier est strictement positif. Les équations obtenues alors pour les indices moindres fournissent des contraintes linéaires sur les premiers coefficients $a_n$ pour qu'ils correspondent aux coefficients d'une série solution de (14.2).

À l'inverse, soit $(a_n)$ une suite vérifiant la récurrence (14.1). Les identités analogues à (14.3) sont maintenant

$$\sum_{n \ge 0} n^k a_n X^n = \left( X \frac{d}{dX} \right)^k A(X), \quad \sum_{n \ge 0} a_{n+k} X^n = (A(X) - a_0 - \cdots - a_{k-1} X^{k-1})/X^k,$$

où $A$ est la série génératrice des coefficients $a_n$ et la notation $(X \frac{d}{dX})^k$ signifie que l'opérateur $X \frac{d}{dX}$ est appliqué $k$ fois. En multipliant (14.1) par $X^n$ et en sommant pour $n$ allant de 0 à $\infty$, puis en multipliant par une puissance de $X$, on obtient donc une équation différentielle linéaire de la forme

$$q_d(X)A^{(d)}(X) + \cdots + q_0(X)A(X) = p(X),$$
où le membre droit provient des conditions initiales. Il est alors possible, quitte à augmenter l'ordre de l'équation de 1, de faire disparaître ce membre droit, par une dérivation et une combinaison linéaire.

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

Ces récritures d'équations différentielles en récurrences et réciproquement seront interprétés comme des morphismes d'anneaux d'opérateurs dans la Partie VI.

**Exemple 14.4** L'équation $A' - X^k A = 0 (k \in \mathbb{N})$ donne la récurrence linéaire $(n+1)a_{n+1} - a_{n-k} = 0$ valide pour tout $n \ge 0$, avec la convention que les $a_n$ d'indice négatif sont nuls. On en déduit que $a_0$ est libre, puis les contraintes $a_1 = \dots = a_k = 0$, et les coefficients suivants sont fournis par la récurrence décalée $(n+k+1)a_{n+k+1} - a_n = 0$, valide pour $n \ge 0$. On reconnaît ainsi les coefficients de $a_0 \exp(X^{k+1}/(k+1))$.

**Exercice 14.1** Retrouver la récurrence satisfaite par la suite $1/n!$ à partir de l'équation $A' = A$, et réciproquement. ■

**Exercice 14.2** Calculer une équation différentielle linéaire satisfaite par la série génératrice des nombres de Fibonacci, définis par $F_{n+2} = F_{n+1} + F_n$, $F_0 = 0$, $F_1 = 1$. ■

**Exercice 14.3** Estimer la complexité d'un algorithme direct utilisant cette idée en nombre d'opérations dans $\mathbb{K}$. ■

Un algorithme plus efficace pour passer d'une équation différentielle d'ordre élevé à la récurrence linéaire satisfaite par les coefficients des solutions est donné en Section 14.1.

## Exemples d'applications

**Exemple 14.5** Pour calculer le coefficient de $X^{1000}$ dans

$$(1 + X)^{1000}(1 + X + X^2)^{500},$$

une méthode efficace part de l'observation que ce polynôme vérifie l'équation différentielle linéaire du premier ordre

$$\frac{y'}{y} = 1000 \frac{1}{1+X} + 500 \frac{2X+1}{1+X+X^2}.$$

Il s'ensuit que les coefficients de ce polynôme vérifient une récurrence d'ordre 3 à coefficients de degré 1. Cette récurrence permet d'obtenir le N-ième coefficient en $O(N)$ opérations arithmétiques, voire $\tilde{O}(\sqrt{N})$ en utilisant la technique de « pas de bébés, pas de géants » présentée au Chapitre 15. Pour la complexité binaire, ce même chapitre présente la méthode de scindage binaire qui est quasi-optimale.
**Exemple 14.6** Pour calculer une racine du polynôme $P_N(X)$ défini par la série génératrice

$$ \sum_{n \ge 0} P_n(X) \frac{Z^n}{n!} = \left( \frac{1+Z}{1+Z^2} \right)^X, $$

lorsque N est grand, il n'est pas nécessaire de calculer ce polynôme. Il suffit d'observer que cette série vérifie une équation différentielle linéaire d'ordre 1 (avec X en paramètre), ainsi que la série génératrice des dérivées des $P_n$, et d'utiliser les récurrences que l'on en déduit sur ces polynômes pour en calculer des valeurs. Ces valeurs permettent alors d'appliquer la méthode de Newton, par exemple pour résoudre le polynôme. Cette idée peut aussi être combinée avec la méthode de scindage binaire du Chapitre 15.

## Test d'égalité

**Lemme 14.2** Si $(u_n)$ et $(v_n)$ sont deux suites solutions de (14.1), et s'il existe $k \in \mathbb{N}$ tel que $u_0 = v_0, \dots, u_{k-1} = v_{k-1}$ et $0 \notin p_d(k-d+1+\mathbb{N})$, alors ces suites sont égales.

*Démonstration.* Par récurrence, puisque $0 \notin p_d(k-d+1+\mathbb{N})$ permet d'inverser le coefficient de tête de la récurrence et donc d'exprimer chaque terme à partir de l'indice $d$ en fonction des précédents. ■

**Exemple 14.7** Pour une récurrence à coefficients constants, il suffit de connaître une borne sur l'ordre de la récurrence pour tester l'égalité de deux solutions. Cette observation est utile en liaison avec les opérations de clôture de la section suivante.

**Corollaire 14.3** Si $f(X)$ et $g(X)$ sont deux séries formelles solutions de (14.2), $f(0) = g(0), \dots, f^{(m-1)}(0) = g^{(m-1)}(0)$ et $q_0(0) \neq 0$, alors ces séries sont égales.

*Démonstration.* D'après (14.3), un terme $cX^iA^{(j)}(X)$ intervient dans la récurrence sur les coefficients sous la forme $c(n-i+1)\cdots(n-i+j)a_{n-i+j}$. L'indice maximal est donc atteint pour $j-i$ maximal et donc pour $j=m$ si $i=0$. La récurrence obtenue, valide pour $n \ge 0$ est donc de la forme $q_0(0)(n+1)\cdots(n+m)a_{n+m} + \cdots = 0$. Les contraintes linéaires sur les conditions initiales jusqu'à l'indice $n_0$ introduit dans la preuve du Théorème 14.1 définissent les coefficients d'indice $m$ à $m+n_0$ à partir des précédents et le lemme précédent s'applique pour les valeurs suivantes. ■

## Algorithme rapide pour la conversion

Vue l'efficacité obtenue grâce au passage de l'équation différentielle à la récurrence, il est souhaitable de maîtriser le coût de cette conversion. Il est possible d'obtenir la récurrence plus efficacement que par la méthode naïve, en mettant en œuvre des techniques qui exploitent la multiplication rapide.
Cas d'un opérateur donné en $\theta = X \frac{d}{dX}$

Si l'équation différentielle linéaire de départ est donnée non pas comme un polynôme en X et $d/dX$, mais comme un polynôme en X et $\theta = X \frac{d}{dX}$ ($\theta$ est parfois appelé opérateur d'Euler), alors la conversion en récurrence est assez facile : l'application de l'opérateur

$$ \sum_{\substack{0 \le j \le m \\ 0 \le i \le d}} p_{i,j} X^j \theta^i, $$

à une série A(X) suivie de l'extraction du coefficient de $X^n$ donne au vu des relations (14.3) la récurrence

$$ \sum p_{i,j} (n-j)^i a_{n-j} = 0. $$

De cette conversion découle le résultat de complexité suivant.

**Proposition 14.4** La récurrence satisfaite par les coefficients des séries solutions d'une équation différentielle linéaire de degré $d$ en $\theta = X \frac{d}{dX}$ et $m$ en X se calcule en $O(mM(d))$ opérations sur les coefficients.

Par rapport au nombre $dm$ de coefficients du résultat, cette complexité est quasi-optimale.

La preuve utilise la formule ci-dessus et l'observation que calculer les coefficients du polynôme $P(X-j)$ connaissant ceux du polynôme $P(X)$ de degré $d$ ne requiert que $O(M(d))$ opérations, par exemple en utilisant la somme composée (Chapitre 3).

### Cas général

Quitte à le multiplier au préalable par une puissance de X égale au plus à son degré en $\partial = d/dX$, il est toujours possible de récrire un polynôme en X et $\partial$ en un polynôme en X et $\theta$. Cette récriture peut elle-même être effectuée assez rapidement.

Une première observation est que de la *commutation*

$$ (\theta - i)X^i = X^i \theta $$

se déduit par multiplication à droite par $\partial^i$ la relation

$$ X^{i+1} \partial^{i+1} = (\theta - i)X^i \partial^i = (\theta - i)(\theta - i + 1) \cdots \theta. $$

Étant donnés des polynômes $p_i(X)$ de degré au plus $m+d$, il s'agit donc maintenant de calculer des polynômes $q_i(X)$ tels que

$$ \sum_{i=0}^{d} p_i(X) X^i \partial^i = \sum_{i=0}^{d} p_i(X) (\theta - i + 1) \cdots \theta = \sum_{i=0}^{d} q_i(X) \theta^i. $$

Récrire le polynôme sous la forme

$$ \sum_{j=0}^{m+d} X^j \sum_{i=0}^{d} p_{i,j} X^i \partial^i $$
s'effectue en nombre linéaire d'opérations et montre qu'il suffit de savoir traiter efficacement le cas où les $p_i$ (et donc aussi les $q_i$) sont constants. La transition des uns vers les autres se calcule alors par évaluation-interpolation sur $\theta = 0, 1, 2, \dots$. Soit P le polynôme à calculer. Les premières identités obtenues par évaluation sont

$$p_0 = q_0, \quad p_0 + p_1 = \sum q_i, \quad p_0 + 2p_1 + 2p_2 = \sum 2^i q_i,$$

et plus généralement

$$e^x \sum p_i X^i = \sum \frac{P(i)}{i!} X^i,$$

ce qui montre que les valeurs de P en $0, \dots, d$ peuvent être obtenues en $O(M(d))$ opérations à partir des coefficients $p_i$, et donc les coefficients de P en $O(M(d) \log d)$ opérations, par interpolation rapide, en utilisant les techniques du Chapitre 5. En résumé, nous avons obtenu le résultat suivant.

**Théorème 14.5** Le calcul de la récurrence linéaire satisfaite par les coefficients des séries solution d'une équation différentielle linéaire d'ordre $d$ à coefficients des polynômes de degré au plus $m$ requiert un nombre d'opérations arithmétiques borné par

$$O((m+d)M(d) \log d).$$

## 14.2 Propriétés de clôture

Les séries différentiellement finies jouissent de propriétés de clôture. En ce sens, elles sont un analogue différentiel des nombres algébriques, et une généralisation des suites récurrentes linéaires à coefficients *constants* (srlcc) évoquées en Section 4.5.

### Somme et produit

**Théorème 14.6** L'ensemble des séries différentiellement finies à coefficients dans un corps $K$ est une algèbre sur $K$. L'ensemble des suites polynomialement récurrentes d'éléments de $K$ est aussi une algèbre sur $K$.

Une conséquence algorithmique fondamentale de ce théorème est le fait que les séries différentiellement finies peuvent être multipliées plus rapidement que les séries arbitraires, en complexité linéaire.

**Corollaire 14.7** Soient $A, B \in K[[X]]$ deux séries différentiellement finies, solutions d'équations différentielles d'ordre $O(1)$ à coefficients polynomiaux de degré $O(1)$. On peut calculer le produit $AB$ mod $X^N$ en $O(N)$ opérations dans $K$.

*Démonstration du théorème.* Les preuves pour les suites et les séries sont similaires. Les preuves pour les sommes sont plus faciles que pour les produits, mais dans le même esprit. Nous ne donnons donc que la preuve pour le produit $C = AB$ de deux séries différentiellement finies $A$ et $B$. Par la formule de Leibniz, toutes les dérivées de
C s'écrivent comme combinaisons linéaires de produits entre une dérivée $A^{(i)}$ de A et une dérivée $B^{(j)}$ de B. Les dérivées de A et de B étant engendrées par un nombre fini d'entre elles, il en va de même pour les produits $A^{(i)}B^{(j)}$, ce qui prouve la finitude différentielle de C.

■

**Exercice 14.4** Prouver le cas du produit de suites polynomialement récurrentes. ■

Cette preuve donne également un algorithme pour trouver l'équation différentielle (resp. la récurrence) cherchée : il suffit de calculer les dérivées (resp. les décalées) successives en les récrivant sur un ensemble fini de générateurs. Une fois en nombre suffisant (c'est-à-dire au pire égal à la dimension plus 1), elles satisfont une relation linéaire. À partir de la matrice dont les lignes contiennent les coordonnées des dérivées successives (resp. des décalés successifs) sur cet ensemble fini de générateurs, cette relation s'obtient par un calcul du noyau de la transposée. Le calcul se ramène donc ultimement à la résolution d'un système linéaire à coefficients dans $K[X]$, ce qui peut se faire efficacement grâce à l'algorithme de Storjohann (Chapitre 11).

**Exemple 14.8** Voici comment prouver (et même découvrir) l'identité

$$ \arcsin(X)^2 = \sum_{k \ge 0} \frac{k!}{\left(\frac{1}{2}\right)! \cdots \left(k + \frac{1}{2}\right)!} \frac{X^{2k+2}}{2k+2}. $$

Le calcul consiste à partir d'une équation satisfaite par $\arcsin(X)$, en déduire une équation satisfaite par son carré, traduire cette équation en récurrence sur les coefficients, et conclure en constatant que cette récurrence est satisfaite par les coefficients de la série.

Le point de départ est la propriété :

$$ (\arcsin(X))' = \frac{1}{\sqrt{1-X^2}}, $$

qui permet de représenter $\arcsin$ par l'équation différentielle $(1 - X^2)y'' - Xy' = 0$ avec les conditions initiales $y(0) = 0, y'(0) = 1$.

Ensuite, en posant $h = y^2$, les dérivations et réductions successives donnent :

$$ \begin{aligned} h' &= 2yy', \\ h'' &= 2y'^2 + 2yy'' = 2y'^2 + \frac{2X}{1-X^2}yy', \\ h''' &= 4y'y'' + \frac{2X}{1-X^2}(y'^2 + yy'') + \left(\frac{2}{1-X^2} + \frac{4X^2}{(1-X^2)^2}\right)yy', \\ &= \left(\frac{4X+2}{1-X^2} + \frac{6X^2}{(1-X^2)^2}\right)yy' + \frac{2X}{1-X^2}y'^2. \end{aligned} $$

Les quatre vecteurs $h, h', h'', h'''$ sont combinaisons linéaires des trois vecteurs $y^2, yy', y'^2$. Ils sont donc liés et une relation de liaison s'obtient en calculant le noyau
de la matrice $3 \times 4$ qui découle de ce système. Le résultat est

$$ (1 - X^2)h''' - 3Xh'' - h' = 0. $$

La récurrence qui s'en déduit est

$$ (n+1)(n+2)(n+3)a_{n+3} - (n+1)^3 a_{n+1} = 0. $$

Comme le facteur $(n+1)$ ne s'annule pas sur $\mathbb{N}$, il est possible de simplifier pour obtenir la récurrence équivalente

$$ (n+2)(n+3)a_{n+3} - (n+1)^2 a_{n+1} = 0. $$

Si le membre droit de l'identité n'était pas connu, il serait possible de le déterminer grâce à cette récurrence à deux termes et aux conditions initiales $h(0) = 0$, $h'(0) = 0$, $h''(0) = 2$. Lorsque le membre droit est donné, la vérification que les coefficients de la série ci-dessus vérifient cette identité est facile.

## Bornes

La preuve du Théorème 14.5 permet aussi de borner l'ordre des équations : l'ordre de l'équation satisfaite par une somme est borné par la somme des ordres des équations des sommants, et l'ordre de l'équation satisfaite par un produit est borné par le produit des ordres. On en déduit alors un algorithme simple de preuve d'identités.

**Exemple 14.9** L'identité de Cassini sur les nombres de Fibonacci s'écrit

$$ F_{n+2}F_n - F_{n+1}^2 = (-1)^n. $$

Pour prouver cette égalité, le point de départ est simplement la récurrence définissant les nombres de Fibonacci :

$$ F_{n+2} = F_{n+1} + F_n, $$

qui exprime que tous les décalés de $F_n$ sont des combinaisons linéaires de $F_n$ et $F_{n+1}$. Les produits qui interviennent dans l'identité de Cassini s'expriment donc $a_{n+1}$i comme combinaison linéaire de $F_n^2$, $F_n F_{n+1}$ et $F_{n+1}^2$ et donc le membre de gauche vérifie une récurrence d'ordre borné par 4. Le membre droit vérifiant une récurrence d'ordre 1, leur différence vérifie une récurrence d'ordre au plus 5. Il n'est pas nécessaire de calculer cette récurrence. Il suffit de vérifier que ses 5 conditions initiales sont nulles. Autrement dit, vérifier l'identité pour $n = 0, \dots, 4$ la prouve !

Si le membre droit n'est pas donné, on calcule la récurrence satisfaite par le membre gauche en récrivant ses décalés sur les générateurs $F_{n+1}F_n, F_n^2, F_{n+1}^2$ :

$$ \begin{aligned} u_n &= F_{n+2}F_n - F_{n+1}^2 = F_{n+1}F_n + F_n^2 - F_{n+1}^2, \\ u_{n+1} &= F_{n+2}F_{n+1} + F_{n+1}^2 - F_{n+2}^2 = F_{n+1}^2 - F_n^2 - F_nF_{n+1} = -u_n. \end{aligned} $$

et on conclut en observant que $u_0 = 1$.
**Exercice 14.5** De la même manière, montrer l'identité $\sin^2 X + \cos^2 X = 1$, avec et sans calcul.

■

## Produit d'Hadamard

**Corollaire 14.8** Si $A = \sum_{n \ge 0} a_n X^n$ et $B = \sum_{n \ge 0} b_n X^n$ sont deux séries différentiellement finies, alors leur produit d'Hadamard $A \odot B = \sum_{n \ge 0} a_n b_n X^n$ l'est aussi.

La preuve est également un algorithme : des deux équations différentielles se déduisent deux récurrences satisfaites par les suites $(a_n)$ et $(b_n)$ ; d'après la section précédente, le produit $(a_n b_n)$ vérifie alors une récurrence linéaire, dont se déduit enfin l'équation différentielle satisfaite par sa série génératrice.

**Exemple 14.10** Les polynômes orthogonaux de Hermite ont pour série génératrice

$$ \sum_{n \ge 0} H_n(X) \frac{Z^n}{n!} = \exp(Z(2X - Z)). $$

À partir de là, la détermination du membre droit de l'identité suivante due à Mehler est entièrement algorithmique :

$$ \sum_{n \ge 0} H_n(X) H_n(Y) \frac{Z^n}{n!} = \frac{\exp\left(\frac{4Z(XY-Z(X^2+Y^2))}{1-4Z^2}\right)}{\sqrt{1-4Z^2}}. $$

**Exemple 14.11** Si l'intégrale

$$ I(t) = \int_{-\infty}^{\infty} e^{-xt^2} f(x) dx, $$

peut être calculée à partir du développement de Taylor de $f$ en intervertissant sommation et intégration, alors $I(t)$ est différentiellement finie si $f$ l'est. Il en va de même pour la transformée de Laplace ou la transformée de Fourier.

**Exercice 14.6** Montrer que dans les mêmes conditions, $f$ est différentiellement finie si $I$ l'est, et donner un algorithme pour calculer une équation différentielle linéaire satisfaite par $f$ à partir d'une équation satisfaite par $I$.

■

## 14.3 Séries algébriques

**Définition 14.3** Une série formelle $F(X) \in K[[X]]$ est une série algébrique lorsqu'il existe un polynôme non nul $P(X, Y) \in K[X, Y]$ tel que $P(X, F(X)) = 0$.
**Théorème 14.9** Si la série F(X) annule un polynôme P(X, Y) ∈ K[X, Y] de degré *d* en Y, où K est un corps de caractéristique 0, alors elle est solution d'une équation différentielle linéaire d'ordre au plus *d*.

*Démonstration.* La preuve est algorithmique. Quitte à diviser d'abord P par son pgcd avec sa dérivée P<sub>Y</sub> par rapport à Y, il est possible de le supposer premier avec P<sub>Y</sub> (car la caractéristique est nulle!). En dérivant P(X, Y) = 0 et en isolant Y', il vient

$$Y' = -\frac{P_X}{P_Y}.$$

Par *inversion modulaire* de P<sub>Y</sub> (voir le Chapitre 6), cette identité se récrit via un calcul de pgcd étendu en

$$Y' = R_1(Y) \bmod P,$$

où R₁ est un polynôme en Y de degré au plus *d* et à coefficients dans K(X). Ceci signifie que Y' s'écrit comme combinaison linéaire de 1, Y, Y², ..., Y<sup>*d*−1</sup> à coefficients dans K(X). Dériver à nouveau cette équation, puis récrire Y' et prendre le reste de la division par P mène à nouveau à une telle combinaison linéaire pour Y'' et plus généralement pour les dérivées successives de Y. Les *d* + 1 vecteurs Y, Y', ..., Y<sup>*d*</sup> sont donc linéairement dépendants et la relation de liaison est l'équation cherchée. ■

**Exemple 14.12** Les dénombrements d'arbres mènent naturellement à des équations algébriques sur les séries génératrices. Ainsi, la série génératrice C(X) des nombres de Catalan (nombre d'arbres binaires à *n* sommets internes) vérifie

$$C = 1 + XC^2;$$

la série génératrice M(X) des nombres de Motzkin (le nombre d'arbres unaires-binaires à *n* sommets internes) vérifie

$$M = 1 + XM + XM^2.$$

Dans les deux cas, il est aisé d'obtenir d'abord une équation différentielle puis une récurrence qui permet de calculer efficacement ces nombres par scindage binaire (Chapitre 15). Dans le cas des nombres de Catalan, la récurrence est d'ordre 1, la suite est donc hypergéométrique et s'exprime aisément.

**Exercice 14.7** Trouver une formule explicite des nombres de Catalan. Récrire les fonctions Γ qui pourraient apparaître dans le résultat en termes de factorielles, puis de coefficient binomial. ■

**Exercice 14.8** À l'aide d'un système de calcul formel, calculer une récurrence linéaire satisfaite par les coefficients de la série Y(X) solution de Y = 1 + XY + XY⁷. ■

Les mêmes arguments que ci-dessus mènent à une autre propriété de clôture des séries différentiellement finies.
**Corollaire 14.10** Si F est une série différentiellement finie et A une série algébrique sans terme constant, alors F∘A est différentiellement finie.

La preuve consiste à observer que les dérivées successives de F∘A s'expriment comme combinaisons linéaires des $F^{(i)}(A)A^j$ pour un nombre fini de dérivées de F (par finitude différentielle) et de puissances de A (par la même preuve que pour le Théorème 14.9). Cette preuve fournit encore un algorithme.

**Exemple 14.13** À l'aide d'un système de calcul formel, calculer une récurrence linéaire satisfaite par les coefficients du développement en série de Taylor de

$$ \exp\left(\frac{1 - \sqrt{1 - 4X}}{2}\right). $$

**Corollaire 14.11** Si A est une série algébrique, alors $\exp \circ A$ est différentiellement finie.

**Exercice 14.9** Faire la preuve.

**Exemple 14.14** — La moyenne arithmétique-géométrique. Si $a$ et $b$ sont deux réels tels que $0 \le b \le a$, les deux suites définies par

$$ a_{n+1} = \frac{a_n + b_n}{2}, \quad b_{n+1} = \sqrt{a_n b_n}, \quad a_0 = a, \quad b_0 = b $$

ont une limite commune dont l'existence se déduit de $b_n \le b_{n+1} \le a_{n+1} \le a_n$ et

$$ a_{n+1}^2 - b_{n+1}^2 = \left( \frac{a_n - b_n}{2} \right)^2. $$

Cette limite est notée $M(a, b)$. L'identité suivante

$$ M(a, b) = \frac{a}{F\left(\frac{1}{2}, \frac{1}{2}; 1; 1 - \frac{b^2}{a^2}\right)}. $$

due à Gauss, où la série hypergéométrique est définie comme

$$ F(a, b; c; z) := \sum_{n \ge 0} \frac{(a)_n (b)_n z^n}{(c)_n n!}, \quad \text{où} \quad (x)_n = x(x+1)\cdots(x+n-1), $$

peut se prouver par les algorithmes de ce chapitre.

Tout d'abord, la fonction M est clairement homogène : pour $\lambda > 0$, $M(\lambda a, \lambda b) = \lambda M(a, b)$ ce qui permet de concentrer l'étude sur la fonction de une variable $M(1, x)$,
dont on peut prouver qu'elle est analytique au voisinage de $x = 1$. Ensuite, l'homo-généité et le fait que $M(a_1, b_1) = M(a_0, b_0)$ avec $a_0 = 1 + x$ et $b_0 = 1 - x$ fournissent l'identité

$$M(1, \sqrt{1-x^2}) = (1+x)M\left(1, \frac{1-x}{1+x}\right).$$

En posant $M(1,x) = 1/A(1-x^2)$, cette équation entraine l'équation fonctionnelle

$$(1+x)A(x^2) - A\left(1 - \frac{(1-x)^2}{(1+x)^2}\right) = 0, \quad (14.4)$$

qui admet une seule solution série avec $A(0) = 1$. Il suffit donc de montrer que cette solution est $f(x) = F(\frac{1}{2}, \frac{1}{2}; 1; x)$.

Pour cela, on va construire une équation différentielle satisfaite par le membre gauche de l'équation ci-dessus évalué en $A(x) = f(x)$, puis observer que quelques conditions initiales suffisent à prouver la nullité. La définition de la série hyper-géométrique entraîne clairement une récurrence d'ordre 1 sur les coefficients

$$\frac{u_{n+1}}{u_n} = \frac{(a+n)(b+n)}{(c+n)(1+n)}.$$

En convertissant cette récurrence en équation différentielle, on obtient que $f(x)$ satisfait

$$4x(x-1)f''(x) + 4(x-2)f'(x) + f(x) = 0.$$

À partir de cette équation, on construit une équation pour $g(x) = f(x^2)$ :

$$x(x^2 - 1)g''(x) + (3x^2 - 1)g'(x) + xg(x) = 0$$

et une autre pour $h(x) = g(1 - (1-x)^2/(1+x)^2)$ :

$$x(x-1)(x+1)^2 h''(x) + (x+1)(x^2 + 2x - 1)h'(x) - (x-1)h(x) = 0 \quad (14.5)$$

(il s'agit de cas simples du corollaire ci-dessus). En utilisant la clôture par somme, on peut alors construire une équation pour $(1+x)g(x) + h(x)$, et on regroupe exactement la même que pour $h$ (ceci peut se voir directement : en remplaçant $g(x)$ dans l'équation qui le définit par $u(x)/(1+x)$, on obtient l'équation satisfaite par $h$).

À ce stade, l'équation (14.5) est donc satisfaite par le membre gauche de l'équation (14.4) lorsque $A$ vaut $f$ et il s'agit de prouver que cette solution de (14.5) est nulle. Le Corollaire 14.3 ne s'applique pas ici puisque le coefficient de tête s'annule en 0. En revanche, on peut passer à la récurrence satisfaite par les coefficients des solutions séries de cette équation (14.5) :

$$(n+3)^2 u_{n+3} + (n^2 + 2n - 1)u_{n+2} - (n^2 + 4n + 2)u_{n+1} - n^2 u_n = 0,$$

ce qui montre qu'il suffit de vérifier les conditions initiales $h(0) = h'(0) = h''(0) = 0$. Pour cela, aucun calcul n'est nécessaire puisque la série A est définie comme solution de (14.4), ce qui conclut.
Dans cet exemple, nous sommes partis de l'équation satisfaite par la série hypergéométrique pour prouver qu'elle vérifiait une équation fonctionnelle. Il est également possible d'utiliser une méthode de coefficients indéterminés pour calculer suffisamment de coefficients de la solution $A(x)$ de (14.4), puis conjecturer cette équation différentielle par un calcul d'approximants de Padé-Hermite (Chapitre 7), et enfin de montrer comme nous l'avons fait que la série définie par cette équation conjecturée est bien solution de l'équation fonctionnelle.

## 14.4 Au-delà

En général, la composition de deux séries différentiellement finies n'est pas différentiellement finie. Voici trois résultats plus forts, dont la preuve repose sur la théorie de Galois différentielle et dépasse le cadre de cet ouvrage.

**Théorème 14.12** Soient A et F deux séries de $K[[X]]$.

1. Les séries F et 1/F sont simultanément différentiellement finies si et seulement si F'/F est algébrique.

2. Les séries F et $\exp(\int F)$ sont simultanément différentiellement finies si et seulement si F est algébrique.

3. Soit A algébrique de genre supérieur ou égal à 1, alors F et $A \circ F$ sont différentiellement finies si et seulement si F est algébrique.

**Exercice 14.10** Prouver le sens « si » de ces trois propriétés. ■

### Exercices

**Exercice 14.11** Parmi les suites données par les termes généraux suivants, indiquer celles qui sont polynomialement récurrentes et celles qui ne le sont pas :

$$ \frac{1}{\binom{2n}{n}}, \cos n, 2^{2n}, H_n = 1 + \frac{1}{2} + \dots + \frac{1}{n}, [\log n]. $$

Parmi les développements en série des fonctions suivantes, indiquer ceux qui sont différentiellement finis et ceux qui ne le sont pas :

$$ \exp(\sqrt[3]{1-X}), \prod_{i \ge 1} \frac{1}{1-X^i}, \sum_{i \ge 0} X^{2i}, \sum_{n \ge 0} F_n^3 \frac{X^n}{n!}, \tan(X). $$

$(F_n)$ désigne le n-ième nombre de Fibonacci, dont la suite vérifie $F_{n+2} = F_{n+1} + F_n$ et $F_1 = F_0 = 1$.) ■

**Exercice 14.12**

1. Soit $(a_n)_{n \ge 0}$ une suite polynomialement récurrente d'éléments d'un corps de caractéristique 0. Montrer que les suites extraites $(a_{2n})_{n \ge 0}$ et $(a_{2n+1})_{n \ge 0}$ sont
également polynomialement récurrentes.

2. Soit $(g_n)_{n \ge 0}$ une suite dont la série génératrice exponentielle est

$$ \sum_{n \ge 0} g_n \frac{X^n}{n!} = \frac{\exp(-X/2 - X^2/4)}{\sqrt{1-X}}. $$

Montrer que la suite $(g_n)_{n \ge 0}$ satisfait une récurrence linéaire homogène, d'ordre au plus 3 et à coefficients polynomiaux de degrés au plus 2. La récurrence n'est pas demandée explicitement.

3. Si $N > 0$, quel est le coût binaire, en fonction de $N$, du calcul de $g_N$ ? [Utiliser des résultats du Chapitre 15.]

■

**Exercice 14.13 — Opérations de clôture pour les équations différentielles linéaires à coefficients constants.** Soient $f(X)$ et $g(X)$ solutions des équations différentielles linéaires homogènes

$$ a_m f^{(m)}(X) + \dots + a_0 f(X) = 0, \quad b_n g^{(n)}(X) + \dots + b_0 g(X) = 0, $$

avec $a_0, \dots, a_m, b_0, \dots, b_n$ des coefficients rationnels.

1. Montrer que $f(X)g(X)$ et $f(X) + g(X)$ sont solutions d'équations différentielles du même type.

Si le polynôme $a_m X^m + \dots + a_0$ se factorise sur $\mathbb{C}$ en

$$ a_m (X - \alpha_1)^{d_1} \dots (X - \alpha_k)^{d_k}, $$

on rappelle qu'une base de l'espace des solutions de

$$ a_m f^{(m)}(X) + \dots + a_0 f(X) = 0 $$

est donnée par $\{e^{\alpha_1 X}, Xe^{\alpha_1 X}, \dots, X^{d_1-1}e^{\alpha_1 X}, \dots, e^{\alpha_k X}, Xe^{\alpha_k X}, \dots, X^{d_k-1}e^{\alpha_k X}\}$.

2. Montrer qu'une équation satisfaite par $f(X) + g(X)$ peut être calculée à l'aide de l'algorithme d'Euclide.

3. Montrer qu'une équation satisfaite par $f(X)g(X)$ peut être obtenue par un calcul de résultant.

■

**Exercice 14.14 — Puissance symétrique d'équation différentielle.** Soient $m$ et $d$ deux entiers naturels et soient $a(X), b(X)$ deux polynômes dans $\mathbb{Q}[X]$ de degrés au plus $d$.

1. Montrer qu'il existe une équation différentielle, notée $\mathcal{E}_m$, linéaire homogène d'ordre $m+1$, à coefficients dans $\mathbb{Q}[X]$, qui admet $\phi(X)^m$ comme solution, quelle que soit la solution $\phi(X)$ de l'équation différentielle

$$ (\mathcal{E}_1) \qquad y''(X) + a(X)y'(X) + b(X)y(X) = 0. $$
On admettra que pour toute base $(\phi_1, \phi_2)$ de solutions de $\mathcal{E}_1$, les fonctions $\phi_1^i \phi_2^{m-i}$, $i = 0, \dots, m$ sont linéairement indépendantes.

2. Montrer que si $a = 0$, alors $\mathcal{E}_2 : y'''(X) + 4b(X)y'(X) + 2b'(X)y(X) = 0$.

3. Expliciter un algorithme pour calculer $\mathcal{E}_m$, en ramenant le calcul final à un calcul sur des matrices de polynômes.

4. Estimer la complexité de cet algorithme en nombres d'opérations dans $\mathbb{Q}$ en fonction de $m$ et $d$.

5. Pour $m$ fixé, on associe à une fonction $y$ de $X$ des fonctions $L_k$ par les valeurs initiales

$$L_0(X) = y(X), \quad L_1(X) = y'(X)$$

et la relation de récurrence

$$L_{k+1}(X) = L'_{k}(X) + ka(X)L_{k}(X) + k(m-k+1)b(X)L_{k-1}(X).$$

a. Montrer que lorsque $y = \phi^m$, pour $0 \le k \le m$,

$$L_k(X) = m(m-1)\dots(m-k+1)\phi^{m-k}(X)\phi'(X)^k,$$

et $L_{m+1}(X) = 0$.

b. Montrer que $L_{m+1}(X) = 0$ n'est autre que l'équation $\mathcal{E}_m$.

c. En déduire des bornes sur les degrés des coefficients de $\mathcal{E}_m$ et un algorithme de complexité $O(m^3M(d))$ pour le calcul de $\mathcal{E}_m$.

6. Montrer que $a, b \in \mathbb{Q}$, alors $\mathcal{E}_m$ est à coefficients constants et qu'on peut calculer $\mathcal{E}_m$ en $O(M(m) \log m)$ opérations dans $\mathbb{Q}$.

■

## Notes

La première édition du *Handbook of Mathematical Functions* [AS92] remonte à 1964 et a été suivie de huit rééditions; le NIST Handbook of Mathematical Functions [Olv+10] est vu comme un successeur de l'œuvre d'origine. Le site web [Slo06] est le prolongement en ligne du livre éponyme *Encyclopedia of Integer Sequences* [SP95] et est interrogeable directement depuis le système SAGE.

Les propriétés de clôture des séries différentiellement finies ont été décrites avec leurs applications par Stanley [Sta80; Sta99] et Lipshitz [Lip89].

L'utilisation des séries différentiellement finies pour les séries algébriques en combinatoire est exploitée à de nombreuses reprises par Comtet dans son livre [Com70], où il utilise l'algorithme décrit dans la preuve du Théorème 14.9. L'histoire de cet algorithme est compliquée. Il était connu d'Abel qui l'avait rédigé dans un manuscrit de 1827 qui n'a pas été publié. Ce manuscrit est décrit (page 287) dans les œuvres complètes d'Abel [Abe92]. Ensuite, ce résultat a été retrouvé par Sir James Cockle en 1860 et popularisé par le révérend Harley en 1862 [Har62]. Quelques années plus tard, il est encore retrouvé par Tannery dans sa thèse, dont le manuscrit est remarquablement clair [Tan74]. Pour le calcul du développement de séries algébriques, il faut tenir compte de l'ordre assez élevé des équations obtenues [Bos+07b], ce qui rend cet algorithme utile pour des très grandes précisions, l'itération de Newton du Chapitre 3 étant préférable pour des précisions modérées.
Les limitations de la dernière section ne sont pas très connues. Elles sont dues à Harris et Sibuya pour la première [HS85] et à Singer pour les deux autres [Sin86].

En ce qui concerne les algorithmes et les implantations, la plupart des algorithmes présentés dans ce chapitre sont implantés dans le package gfun de MAPLE [SZ94]. L’algorithme rapide de la Section 14.1 provient essentiellement de la thèse de doctorat de Bostan [Bos03], qui donne une variante légèrement plus efficace (d’un facteur constant).

Bibliographie

<table>
  <tr>
    <td>Abe92</td>
    <td>ABEL, Niels Henrik (1992). Œuvres complètes. Tome II. Réimpression de la seconde édition (1881). Éditions Jacques Gabay. URL : http://www.gallica.fr.</td>
  </tr>
  <tr>
    <td>AS92</td>
    <td>ABRAMOWITZ, Milton et Irene A. STEGUN, éds. (1992). Handbook of mathematical functions with formulas, graphs, and mathematical tables. Réimpression de l'édition de 1972. Dover Publications Inc.</td>
  </tr>
  <tr>
    <td>Bos+07b</td>
    <td>BOSTAN, Alin, Frédéric CHYZAK, Grégoire LECERF, Bruno SALVY et Éric SCHOST (2007). « Differential equations for algebraic functions ». In : IS-SAC'07 : International Symposium on Symbolic and Algebraic Computation. Éd. par Dongming WANG. ACM Press, p. 25-32.</td>
  </tr>
  <tr>
    <td>Bos03</td>
    <td>BOSTAN, Alin (2003). « Algorithmique efficace pour des opérations de base en calcul formel ». Thèse de Doctorat. École polytechnique.</td>
  </tr>
  <tr>
    <td>Com70</td>
    <td>COMTET, L. (1970). Analyse combinatoire. 2 volumes. Presses Universitaires de France.</td>
  </tr>
  <tr>
    <td>Har62</td>
    <td>HARLEY, Rev. Robert (1862). « On the theory of the transcendental solution of algebraic equations ». In : Quarterly Journal of Pure and Applied Mathematics, vol. 5, p. 337-360.</td>
  </tr>
  <tr>
    <td>HS85</td>
    <td>HARRIS, William A. et Yasutaka SIBUYA (1985). « The reciprocals of solutions of linear ordinary differential equations ». In : Advances in Mathematics, vol. 58, n°2, p. 119-132.</td>
  </tr>
  <tr>
    <td>Lip89</td>
    <td>LIPSHITZ, L. (1989). « D-finite power series ». In : Journal of Algebra, vol. 122, n°2, p. 353-373.</td>
  </tr>
  <tr>
    <td>Olv+10</td>
    <td>OLVER, Frank W. J., Daniel W. LOZIER, Ronald F. BOISVERT et Charles W. CLARK, éds. (2010). NIST handbook of mathematical functions. U.S. Department of Commerce, National Institute of Standards and Technology, Washington, DC ; Cambridge University Press, Cambridge.</td>
  </tr>
  <tr>
    <td>Sin86</td>
    <td>SINGER, Michael F. (1986). « Algebraic relations among solutions of linear differential equations ». In : Transactions of the American Mathematical Society, vol. 295, n°2, p. 753-763.</td>
  </tr>
  <tr>
    <td>Slo06</td>
    <td>SLOANE, N. J. A. (2006). The on-line encyclopedia of integer sequences. URL : http://www.research.att.com/~njas/sequences/.</td>
  </tr>
  <tr>
    <td>SP95</td>
    <td>SLOANE, N. J. A. et Simon PLOUFFE (1995). The encyclopedia of integer sequences. Academic Press, Inc., San Diego, CA.</td>
  </tr>
  <tr>
    <td>Sta80</td>
    <td>STANLEY, R. P. (1980). « Differentiably finite power series ». In : European Journal of Combinatorics, vol. 1, n°2, p. 175-188.</td>
  </tr>
</table>
Sta99 STANLEY, Richard P. (1999). *Enumerative combinatorics*. Vol. 2. Cambridge University Press.

SZ94 SALVY, Bruno et Paul ZIMMERMANN (1994). « Gfun : a Maple package for the manipulation of generating and holonomic functions in one variable ». In : *ACM Transactions on Mathematical Software*, vol. 20, n°2, p. 163-177.

Tan74 TANNERY, Jules (1874). « Propriétés des intégrales des équations différentielles linéaires à coefficients variables ». Thèse de doctorat ès sciences mathématiques. Faculté des Sciences de Paris. URL : http://gallica.bnf.fr.
# 15. Récurrences linéaires à coefficients polynomiaux : Nᵉ terme, N premiers termes

## Résumé

Pour calculer les N premiers termes d'une suite polynomialement récurrente, l'algorithme direct par déroulement de la récurrence est quasi-optimal vis-à-vis de N : sa complexité arithmétique (resp. binaire) est linéaire en N (resp. quasi-quadratique en N). Le N-ième terme d'une telle suite peut être calculé plus rapidement que l'ensemble des N premiers termes, en essentiellement $\sqrt{N}$ opérations arithmétiques et N opérations binaires. Ces calculs rapides sont la base de l'évaluation numérique à grande précision en complexité quasi-optimale pour les solutions d'équations différentielles linéaires.

Le Chapitre 4 a montré comment calculer le N-ième terme d'une suite donnée par une récurrence à coefficients constants en complexité arithmétique $O(\log N)$. Ce chapitre attaque le cadre plus général des récurrences à coefficients polynomiaux.

Le calcul du N-ième terme ou des N premiers termes d'une suite polynomialement récurrente, donnée par la récurrence à coefficients polynomiaux

$$p_r(n)u_{n+r} + \cdots + p_0(n)u_n = 0, \quad n \in \mathbb{N}, \qquad (15.1)$$

intervient fréquemment en combinatoire, ou pour calculer des troncatures de séries, ainsi que comme étape clé dans des algorithmes de recherche de solutions polynomiales d'équations fonctionnelles linéaires, dans un algorithme de factorisation déterministe d'entiers, ou encore en cryptologie, dans des questions de comptage de points sur des courbes. Ce calcul intervient aussi dans l'évaluation numérique efficace des séries différentiellement finies au Chapitre 17.
**Définition 15.1** La récurrence (15.1) est dite régulière sur l'intervalle $\{0,...,N\}$ lorsque le coefficient de tête $p_r$ de la récurrence (15.1) ne s'annule sur aucun des entiers $0,..., N-r$. Elle est dite régulière si cette propriété est vraie pour tout N, et singulière sinon.

En dehors de sa dernière Section 15.4, ce chapitre traite de récurrences régulières.

Si les coefficients $p_i(n)$ de la récurrence (15.1) ont degré au plus $d$ en $n$, l'algorithme naïf par déroulement de la récurrence a une complexité arithmétique en $O(rdN)$ et une complexité binaire en $O(rN^2M_Z(d \log N))$ (voir §15.1). Ces complexités sont quasi-optimales vis-à-vis de N pour le calcul de l'ensemble des N premiers termes.

Pour le calcul du N-ième terme seul, les algorithmes présentés dans ce chapitre n'ont pas une aussi bonne complexité que pour des coefficients constants, mais cette fois encore, ils sont plus efficaces que le simple déroulement de la récurrence. Un nouveau phénomène distingue ces algorithmes de ceux qui ont été présentés jusqu'ici : les idées qui permettent le calcul du N-ième terme en bonne complexité arithmétique ne sont pas les mêmes que pour abaisser la complexité binaire. Ainsi, l'algorithme par « pas de bébés, pas de géants » de la Section 15.2 donne un gain en racine carrée de N sur la complexité arithmétique, par rapport à la méthode naïve, mais ne gagne rien sur la complexité binaire. Quant à l'algorithme par scindage binaire de la Section 15.3, il n'abaisse en rien la complexité arithmétique, et améliore pourtant la complexité binaire jusqu'à la rendre quasi-linéaire en la taille de sa sortie.

## 15.1 Calcul naïf de N! et de suites polynomialement récurrentes

### Cas de la factorielle

La définition de la factorielle comme produit,

$$N! = 1 \times 2 \times \dots \times N,$$

montre que N! peut être calculé en $N-1$ opérations arithmétiques, ce qui est optimal pour le calcul des N premières factorielles.

L'estimation de la complexité binaire de cette méthode repose sur la formule de Stirling :

$$\log N! = N \log N - N + \frac{1}{2} \log N + O(1), \quad N \to \infty. \quad (15.2)$$

En particulier, l'équivalence $\log N! \sim N \log N$ donne la taille de N!. La k-ième étape du produit multiplie $k!$ par $k+1$, soit donc un entier de taille $\log k! \sim k \log k$ avec un entier de taille $\log(k+1) \sim \log k$. Ce produit déséquilibré coûte donc moins de $ckM_Z(\log k)$ opérations binaires pour une certaine constante $c$. Le calcul complet de N! par la méthode naïve effectue donc moins de

$$\sum_{k=1}^{N} ckM_{Z}(\log k) = O(N^2 M_{Z}(\log N))$$

opérations binaires. La formule de Stirling montre que cette complexité binaire est quasi-optimale pour le calcul conjoint des nombres 1!, 2!, ..., N!.
Cas général — complexité arithmétique

Pour une suite polynomialement récurrente définie l'équation (15.1), le calcul par
l'algorithme direct de $u_{n+r}$ à partir des $r$ valeurs précédentes demande $O vd$ opérations arithmétiques pour l'évaluation des polynômes (par exemple par le schéma de Horner) puis $O(r)$ opérations pour effectuer la combinaison linéaire des $u_{n+i}$. Le calcul naïf de $u_0, \dots, u_N$ à partir des valeurs initiales $u_0, \dots, u_{r-1}$ demande donc $O vd N)$ opérations arithmétiques.

Une amélioration de cette borne vient de l'observation que les coefficients $p_i(n)$ peuvent être évalués sur les entiers $n = 0, 1, \dots, N$ en utilisant des techniques d'évaluation multipoint rapide, avant de procéder au déroulement de (15.1). Cela permet d'abaisser la complexité arithmétique du calcul de $O vd(N)$ à

$$
O(r \max(N, d)M(d) \log d/d),
$$

voire à O(r max(N, d)M(d)/d) en utilisant l'observation du Corollaire 4.12 page 89, ce
qui revient dans les deux cas à Œ(r max(N, d)) si la multiplication polynomiale est à
base de FFT.

Cas général — complexité binaire

Lorsque les coefficients des polynômes $p_i$ sont entiers, la complexité binaire de l'al-
gorithme naïf découle de nouveau essentiellement de la croissance de $u_N$. Pour estimer
celle-ci, il est commode d'adopter une vision matricielle et de faire apparaître le ra-
tionnel $u_N$ comme un quotient d'entiers. La suite vectorielle des $U_n = {}^t(u_n, \dots, u_{n+r-1})$
vérifie la récurrence du premier ordre

$$
U_{n+1} = \frac{1}{p_r(n)} A(n) U_n \quad \text{avec} \quad A(n) =
\begin{pmatrix}
0 & p_r(n) & & & \\
& 0 & p_r(n) & & \\
& & \ddots & \ddots & \\
& & & 0 & p_r(n) \\
-p_0(n) & & & & -p_{r-1}(n)
\end{pmatrix}
$$

Avec la suite hypergéométrique ($w_n$) définie par la condition initiale $w_0 = 1$ et la
récurrence

$$
w_{n+1} = p_r(n)w_n, \quad (n \in \mathbb{N}),
$$

qui permet de récrire $U_n$ sous la forme

$$
U_n = \frac{1}{w_n} (A(n-1) \cdots A(0)) U_0,
$$

le rationnel $u_N$ s'écrit $u_N = v_N/w_N$, où $v_N$ est le premier coefficient du vecteur $(A(N-1)\cdots A(0))U_0$.

La croissance de $v_N$ est estimée à l'aide de la norme d'indice 1 et de sa norme
subordonnée. Pour un vecteur U de dimension r et une matrice M de dimension r × r,
il s'agit de

$$
\|U\| = \sum_{j=1}^{r} |v_j| \quad \text{et} \quad \|M\| = \max_{1 \le j \le r} \sum_{i=1}^{r} |m_{i,j}|.
$$
Il est alors facile de vérifier (et c'est le sens de la norme subordonnée) que pour tout $u$, $||Mu|| \le ||M|| ||u||$. Soit $\ell$ un majorant commun pour la taille binaire des coefficients des $p_i$, $0 \le i \le r$, et de même $\ell'$ un majorant de la taille des conditions initiales $u_j$, $0 \le j < r$. Il s'ensuit que la norme de $A(n)$ est bornée par

$$||A(n)|| \le r 2^{\ell}(d+1)(n+1)^d,$$

et celle de $U_0$ par $r2^{\ell'}$. Les majorations

$$|v_N| \le ||A(N-1)|| \cdots ||A(0)|| ||U_0|| \le r^{N+1} 2^{\ell N + \ell'} (d+1)^N (N!)^d$$

fournissent la borne $O(dN \log N + N\ell + \ell' + N \log r)$ sur la taille $\log|v_N|$. Par le même raisonnement, la taille du dénominateur $w_N$ est du même ordre (mêmes formules pour $r=1$), si bien que par le même type d'argument que pour la factorielle, la complexité binaire du calcul naïf de $u_N$ est $O(rN^2 M_Z(d \log N + \ell + \log r + \ell'/N))$.

## Conclusion

Les mêmes raisonnements sur la taille et sur la complexité s'appliquent sans modification au cas de récurrences inhomogènes

$$p_r(n)u_{n+r} + \cdots + p_0(n)u_n = q(n), \quad n \in \mathbb{N}, \qquad (15.3)$$

où $q$ est lui aussi un polynôme. Nous avons donc obtenu le résultat suivant. L'algorithme correspondant dans le cas légèrement plus général des récurrences singulières est présenté en Algorithme 15.3 page 281.

**Proposition 15.1** Étant donnée la récurrence (15.3) supposée régulière et des conditions initiales, le calcul des N premiers éléments de la suite peut être effectué en

$$O\left(r \max(N, d) \frac{M(d)}{d}\right)$$

opérations arithmétiques, où $d$ est une borne sur le degré des polynômes coefficients.

Si les coefficients de ces polynômes sont des entiers de taille bornée par $\ell$ et que les conditions initiales sont des rationnels de numérateurs et dénominateurs bornés par $O(dN \log N + \ell N + N \log r)$, alors la complexité binaire de ce calcul est

$$O(rN^2 M_Z(d \log N + \ell + \log r)),$$

la taille du résultat étant $O(N^2(d \log N + \ell + \log r))$.

## 15.2 Pas de bébés, pas de géants

On a vu au Chapitre 4 que le N-ième terme d'une récurrence linéaire à coefficients constants peut être calculé en complexité arithmétique $O(\log N)$. Dans le cas des
coefficients polynomiaux, les choses se compliquent : à ce jour, on ne connaît pas d'algorithmes polynomial en log N (il est fort possible qu'un tel algorithme n'existe pas). L'exemple typique en est le calcul du N-ième terme de la factorielle $u_N = N!$, qui vérifie la récurrence à coefficients polynomiaux $u_{n+1} = (n+1)u_n$ pour $n \ge 0$. Une solution efficace exploite le théorème « évaluation interpolation » (Théorème 5.1). Elle utilise la technique des « pas de bébés, pas de géants » et requiert un nombre d'opérations arithmétiques en $\sqrt{N}$, à des facteurs logarithmiques près. Pour simplifier la présentation, supposons que N est un carré parfait. L'idée de l'algorithmme est de poser

$$P(X) = (X + 1)(X + 2)\cdots(X + N^{1/2}),$$

afin d'obtenir la valeur de $u_N$ à l'aide de l'équation

$$u_N = \prod_{j=0}^{N^{1/2}-1} P(jN^{1/2}). \qquad (15.4)$$

Cette égalité suggère la procédure suivante :

1. **Pas de bébés.** Calculer les coefficients de P. Ceci requiert $O(M(\sqrt{N}) \log N)$ opérations arithmétiques (en construisant un arbre binaire de feuilles $X + i$ comme en Section 5.4).

2. **Pas de géants.** Évaluer P sur les points $0, \sqrt{N}, 2\sqrt{N}, \dots, (\sqrt{N}-1)\sqrt{N}$ et retrouver la valeur de $u_N$ à l'aide de l'équation (15.4). En utilisant les techniques d'évaluation multipoint rapide (Section 5.4), ceci peut se faire également en $O(M(\sqrt{N}) \log N)$ opérations arithmétiques.

Le coût total de cet algorithme est $O(M(\sqrt{N}) \log N)$ opérations arithmétiques. Si la FFT est utilisée pour la multiplication des polynômes, le gain par rapport à la méthode directe est de l'ordre de $\sqrt{N}$, à des facteurs logarithmiques près, typique pour la technique des « pas de bébés, pas de géants ». La technique gagne encore par rapport à l'algorithmme naïf si l'on emploie une multiplication par Karatsuba, mais ne présente aucun intérêt avec la multiplication naïve. Cette technique se généralise aux suites hypergéométriques quelconques ainsi qu'au calcul d'un terme d'une suite récurrente linéaire à coefficients polynomiaux, comme détaillé dans l'Algorithme 15.1.

**Théorème 15.2** Étant donné la récurrence linéaire (15.1) des conditions initiales et un entier N, si la récurrence est régulière sur $\{0, \dots, N\}$, alors l'Algorithme 15.1 est correct. Il calcule le N-ième terme d'une suite solution en

$$O(r^{\theta} M(\sqrt{dN}) \log(dN))$$

opérations arithmétiques, où d est une borne sur le degré des coefficients et $\theta$ est un exposant réalisable pour la complexité du produit de matrices (voir Chapitre 8).

**Démonstration.** L'algorithme calcule les produits $M = A(n-r) \cdots A(0)$ et $w = p_r(N-r) \cdots p_r(0)$ (pour ce dernier la récurrence est simplement $u_{n+1} = p_r(n)u_n$), puis le vecteur $^t(u_{N-r+1}, \dots, u_N)$ d'après la définition de A. La correction est donc claire.
**Entrée** Une récurrence de la forme (15.1) avec $\deg p_i \le d$, $i = 0, \dots, r$, des conditions initiales $u_0, \dots, u_{r-1}$, un entier N.

**Sortie** Les termes $u_{N-r+1}, \dots, u_N$ de la solution.

1. Calculer $P(X) = A(X+m-1)\cdots A(X+1)A(X)$, pour $m = [(N/d)^{1/2}]$, où A est la matrice définie page 271.

2. Évaluer les coefficients des matrices $P(0), P(m), \dots, P(\lfloor \frac{N-r-m}{m} \rfloor m)$.

3. Évaluer ceux de $A(N-r), A(N-r-1), \dots, A(\lfloor \frac{N-r}{m} \rfloor m)$.

4. Calculer le produit de matrices scalaires $M = A(N-r) \cdots A(\lfloor \frac{N-r}{m} \rfloor m) P(\lfloor \frac{N-r}{m} \rfloor m) \cdots P(0)$.

5. Si $p_r(n) \neq 1$, calculer $w = p_r(N-r) \cdots p_r(0)$ par ce même algorithme; sinon $w = 1$.

6. Renvoyer $\frac{1}{w} M \cdot t(u_0, \dots, u_{r-1})$.

**Algorithme 15.1** – Pas de bébés, pas de géants pour les récurrences linéaires.

Le produit de matrices de polynômes à l'étape (1) a degré au plus $md$. Il est effectué récursivement par un arbre des sous-produits comme au chapitre 5. Son coût est donc borné par $O(\text{MM}(r, md) \log m)$ opérations arithmétiques. L'évaluation à l'étape (2) est effectuée en réalisant séparément les évaluations multipoints des $r^2$ coefficients de la matrice P par la technique du Corollaire 4.12 page 89 en $O(r^2 \frac{N}{m} M(md)/(md))$ opérations arithmétiques. L'étape (3) effectue au plus $m$ évaluations de matrices de degré au plus $d$, et est donc dominée par la précédente (qui calcule environ $md$ évaluations de polynômes de degré $md$). L'étape (4) multiple $O(m + N/m)$ matrices scalaires en $O(\text{MM}(r)(m + N/m))$ opérations. L'étape suivante a une complexité au plus égale à celle des précédentes. En résumé, la complexité arithmétique de cet algorithme est donc majorée par

$$ 
\begin{align*}
& O(\text{MM}(r, md) \log m + r^2 \frac{N}{m} \frac{M(md)}{md} + \text{MM}(r) \sqrt{Nd}) \\
&= O\left( r^\theta M(\sqrt{Nd}) \log(Nd) + r^2 \sqrt{Nd} \frac{M(\sqrt{Nd})}{\sqrt{Nd}} + r^\theta \sqrt{Nd} \right), \\
&= O(r^\theta M(\sqrt{Nd}) \log(Nd)). \quad \blacksquare
\end{align*}
$$

**Exercice 15.1** Étant donné un polynôme $f \in \mathbb{K}[X]$ de degré 2 et un entier positif N, borner le nombre d'opérations dans $\mathbb{K}$ suffisantes pour déterminer le coefficient de $X^N$ du polynôme $f^N$. (Indication : La solution consiste à calculer par exponentiation binaire tous les coefficients $u_0, \dots, u_N$ de $f^N$ mod $X^{N+1}$. Une approche encore plus rapiide repose sur le fait que les $u_n$ satisfont à une récurrence à coefficients polynomiaux d'ordre 2.) ■

## Factorisation déterministe des entiers

Pour factoriser un entier N, tester tous les diviseurs plus petits que $\sqrt{N}$ a une complexité binaire au moins linéaire en $\sqrt{N}$. Afin d'accélérer ce calcul, on peut
rassembler tous les entiers plus petits que $\sqrt{N}$ en $\sqrt[4]{N}$ blocs, chacun contenant $\sqrt[4]{N}$ entiers consécutifs. Soit $c$ de l'ordre de $\sqrt[4]{N}$ et $f_0 = 1 \cdots c \bmod N$, $f_1 = (c+1)\cdots(2c) \bmod N$, ..., $f_{c-1} = (c^2-c+1)\cdots(c^2) \bmod N$. Si les valeurs $f_0, \dots, f_{c-1}$ sont connues, alors il devient facile de déterminer un facteur de N, en prenant des pgcd de $f_0, \dots, f_{c-1}$ avec N : la complexité binaire associée est $O(\sqrt[4]{N} M_Z(\log N) \log \log N)$.

Ainsi, la principale difficulté est de calculer les valeurs $f_i$. Pour ce faire, on prend $\mathbb{A} = \mathbb{Z}/N\mathbb{Z}$ et F le polynôme $(X+1)\cdots(X+c) \in \mathbb{A}[X]$, qu'on évalue en les points $0, c, 2c, \dots, c(c-1)$ en $O(M(c) \log c)$ opérations de $\mathbb{A}$. Par le théorème « évaluation interpolation », puisque $c$ est d'ordre $\sqrt[4]{N}$, la complexité pour calculer les $f_i$ est de $O(M(\sqrt[4]{N}) \log N)$ opérations modulo N, soit $O(M(\sqrt[4]{N}) M_Z(\log N) \log N)$ opérations binaires. Ce terme est dominant et est donc aussi la complexité totale du calcul.

## 15.3 Scindage binaire

### Cas de la factorielle

D'après la formule de Stirling (15.2), la factorielle de $N/2$ fait à peu près la moitié de la taille de $N!$. Il en va donc de même pour le produit des entiers de $N/2$ à $N$. En utilisant cette observation récursivement, l'idée de la méthode de scindage binaire consiste à exploiter la multiplication rapide en équilibrant les produits, calculant $P(a, b) = (a+1)(a+2)\cdots b$ récursivement par

$$P(a, b) = P(a, m)P(m, b) \quad \text{où} \quad m = \left\lfloor \frac{a+b}{2} \right\rfloor.$$

Soit $C(a, b)$ le coût binaire du calcul de $P(a, b)$. Il résulte immédiatement de la méthode que ce coût vérifie l'inégalité

$$C(a, b) \le C(a, m) + C(m, b) + M_Z(\max(\log P(a, m), \log P(m, b))),$$

Sous l'hypothèse raisonnable que la complexité de la multiplication d'entiers est croissante avec la taille des entiers, le coût du calcul de $P(a, m)$ est inférieur au coût du calcul de $P(m, b)$, d'où la nouvelle inégalité

$$C(a, b) \le 2C(m, b) + M_Z(\log P(m, b)).$$

À ce stade, le théorème « diviser pour régner » n'est pas suffisant pour conclure, mais l'utilisation de l'inégalité précédente pour $n$ une puissance de 2 donne les inégalités successives

$$
\begin{align*}
C(0, n) &\le 2C(n/2, n) + M_Z(\log P(n/2, n)) \\
&\le 2C(n/2, n) + M_Z\left(\frac{n}{2} \log n\right) \\
&\le 4C(3n/4, n) + 2M_Z\left(\frac{n}{4} \log n\right) + M_Z\left(\frac{n}{2} \log n\right) \\
&\le 4C(3n/4, n) + 2M_Z\left(\frac{n}{2} \log n\right)
\end{align*}
$$
$$ \le \dots \\ \le 2^k C(n - 2^{-k}n, n) + k M_Z\left(\frac{n}{2} \log n\right) $$

pour tout entier positif $k$, où la quatrième ligne découle de la sous-additivité de la fonction $M_Z$. En choisissant finalement d'arrêter la récursion lorsque $n - 2^{-k}n$ et $n$ diffèrent d'au plus un, ce qui impose $k$ de l'ordre de $\log n$, on aboutit à la borne

$$ C(0, n) \le O(n) + M_Z\left(\frac{n}{2} \log n\right) \log n $$

donc à une complexité binaire pour le calcul de $N!$ dans $O(M_Z(N \log N) \log N)$.

Lorsque la multiplication entière utilisée repose sur la FFT, cette complexité devient $O(N \log^3 N \log \log N)$; si la multiplication entière utilisée est d'exposant $\alpha$ strictement plus grand que 1, alors elle s'écrit $O((N \log N)^\alpha)$: au lieu d'utiliser simplement la sous-additivité de la multiplication, on majore la somme des $M_Z$ par une somme géométrique multipliée par le dernier sommant.

## Récurrences d'ordre 1

La factorielle de la section précédente suit la récurrence $u_{n+1} = (n+1)u_n$. La méthode se généralise tout d'abord aux solutions de récurrences de la forme

$$ u_{n+1} = p(n)u_n, \quad p \in \mathbb{Z}[X]. $$

Si $p$ a degré $d$, alors $\log p(N)$ est dans $O(d \log N)$, si bien que la taille de $u_N$ est au plus $O(dN \log N)$.

De même, pour toute récurrence de la forme

$$ u_{n+1} = \frac{p(n)}{q(n)} u_n, \quad p, q \in \mathbb{Z}[X], $$

on peut pour calculer le terme $u_N$ appliquer séparément la même technique de scindage binaire sur le numérateur et le dénominateur. Si $d$ est le maximum des degrés de $p$ et $q$, le calcul produit deux entiers de taille $O(dN \log N)$ en un nombre d'opérations binaires en $O(M_Z(dN \log N) \log N)$.

À ce stade, on dispose d'un rationnel $u_N = P/Q$ dont le numérateur et le dénominateur sont des entiers potentiellement grands. Une évaluation numérique de ce rationnel à $2^{-e}$ près consiste à calculer $2^{-e}[2^e P/Q]$. La méthode de Newton (Théorème 4.2) permet de calculer l'inverse de $Q$ à précision $e$ en $O(M_Z(e))$ opérations. Il est donc possible de terminer le calcul par une approximation de $u_N$ à précision $e$ jusqu'à $O(dN \log N)$ sans affecter la complexité. Cette observation est exploitée dans les sections suivantes. (Le raisonnement est le même dans une base différente de 2; le changement de base a un coût légèrement plus que linéaire.)

## Calcul de $e = \exp(1)$

Le point de départ est la suite $(e_n)$ donnée par $e_n = \sum_{k=0}^{n} \frac{1}{k!}$. Cette suite converge vers $e$. Plus précisément, il est classique que le terme de reste dans $e - e_n$ se majore par
une série géométrique de sorte que

$$0 \le e - e_n < \frac{1}{(n+1)!} \left( 1 + \frac{1}{n+1} + \frac{1}{(n+1)^2} + \cdots \right) = \frac{1}{n n!}.$$

Pour calculer $N$ décimales de $e$ par cette série, il suffit de rendre $\frac{1}{nn!}$ inférieur à $10^{-N}$. Il s'ensuit qu'il suffit de choisir $N$ de sorte que $nn!$ et $10^N$ soient du même ordre, c'est-à-dire d'avoir $N$ proportionnel à $n \log n$. Il suffit donc de $n = O(N/\log N)$ termes de la série.

La suite $(e_n)_{n \ge 0}$ vérifie $e_n - e_{n-1} = 1/n!$ et donc

$$n(e_n - e_{n-1}) = e_{n-1} - e_{n-2}.$$

Cette récurrence se récrit

$$\underbrace{\binom{e_n}{e_{n-1}} = \frac{1}{n} \binom{n+1}{n} \binom{-1}{0}}_{A(n)} \binom{e_{n-1}}{e_{n-2}} = \frac{1}{n!} A(n) A(n-1) \cdots A(1) \binom{1}{0}.$$

Le calcul du produit de matrices peut alors être effectué par scindage binaire. Le calcul de $e_N$ par ce procédé demande donc $O(M_Z(N \log N) \log N)$ opérations binaires.

Pour calculer $n$ décimales de $e$, la première étape ci-dessus construit deux entiers de taille $O(N \log N) = O(n)$ en $O(M_Z(n) \log n)$ opérations binaires. La division finale ne demande que $O(M_Z(n))$ opérations par l'algorithme de Newton. L'ensemble du calcul est donc quasi-optimal. (En outre, le $\log n$ n'est pas présent pour des multiplications comme celle de Karatsuba dont l'exposant de complexité est supérieur à 1.)

La généralisation de ces calculs à l'évaluation numérique de solutions d'équations différentielles linéaires est présentée au Chapitre 17.

## Suites polynomialement récurrentes

L'Algorithmme 15.2 traite le cas général d'une suite polynomialement récurrente vérifiant la récurrence (15.1) avec des coefficients polynomiaux $p_i$ à coefficients entiers.

**Théorème 15.3** Si $N$ est un entier positif, $d$ est une borne sur les degrés des polynômes $p_i \in \mathbb{Z}[n]$ coefficients de la récurrence (15.1) supposée régulière sur $\{0, \dots, N\}$, et si $\ell$ est une borne sur la taille des coefficients des $p_i$ telle la taille des numérateurs et dénominateurs des conditions initiales supposées rationnelles soient bornées par $O(dN \log N + \ell N + N \log r)$, alors la complexité binaire de la méthode du scindage binaire pour calculer le $N$-ième terme d'une suite solution est

$$O(r^{\theta} M_{\mathbb{Z}}(dN \log N + \ell N + N \log r) \log N),$$

où $\theta$ est un exposant réalisable pour la complexité du produit de matrices (voir Chapitre 8).
**Entrée** Une récurrence de la forme (15.1) avec $p_i \in \mathbb{Z}[X]$, $i = 0, \dots, r$, des conditions initiales rationnelles $u_0, \dots, u_{r-1}$, un entier N.

**Sortie** Les termes $u_{N-r+1}, \dots, u_N$ de la solution.

1. Évaluer $A(N-r), \dots, A(1), A(0)$, où A est la matrice définie en page 271.

2. Calculer récursivement le produit $M_{N-r+1,0}$ de ces matrices par la formule

$$ M_{i,j} = \begin{cases} A(j), & \text{si } i = j + 1, \\ M_{i,k} \cdot M_{k,j}, & \text{avec } k = \lfloor(i + j)/2\rfloor, \text{ sinon.} \end{cases} $$

3. Si $p_r(n) \neq 1$, calculer $w = p_r(N-r) \cdots p_r(0)$ par ce même algorithme; sinon $w = 1$.

4. Renvoyer $\frac{1}{w}M^t(u_0, \dots, u_{r-1})$.

Algorithme 15.2 – Algorithme de scindage binaire.

*Démonstration.* La complexité du calcul de la matrice de passage est comme ci-dessus, puis le produit matrice vecteur avec les conditions initiales de taille $O(dN \log N + \ell N + N \log r)$ reste dans la complexité annoncée. ■

Comme précédemment, cette borne de complexité peut être améliorée à

$$ O(r^{\theta} M_{\mathbb{Z}}(dN \log N + \ell N + N \log r)) $$

si l'exposant $\alpha$ de la complexité $M_{\mathbb{Z}}(m) = O(m^{\alpha})$ est supérieur à 1.

## 15.4 Récurrences singulières

Les résultats obtenus dans ce chapitre supposent que les récurrences sont régulières. L'extension de ces algorithmes aux récurrences singulières est possible, mais il faut d'abord préciser le problème en analysant les espaces de solutions qui apparaissent dans ce cas.

### Espace des suites solutions

**Lemme 15.4** L'ensemble des suites solutions de la récurrence

$$ p_r(n)u_{n+r} + \cdots + p_0(n)u_n = 0 \qquad (15.5) $$

pour $n \in \mathbb{N}$, forme un espace vectoriel sur $\mathbb{K}$.

*Démonstration.* La combinaison linéaire de solutions reste une solution, il s'agit donc bien d'un sous-espace vectoriel de $\mathbb{K}^N$. ■
Dans le cas des récurrences à coefficients constants, la dimension de l'espace des solutions est exactement l'ordre de l'équation. Ici, la dimension de cet espace est plus délicate, comme le montre l'exemple suivant.

**Exemple 15.1** Si la suite $(u_n)$ obéit à la récurrence

$$n(n-1)u_{n+3} - (n-1)u_{n+2} + (n+1)u_{n+1} + 2nu_n = 0, \quad n \in \mathbb{N},$$

alors les seules contraintes induites par la récurrence pour $n = 0$ et $1$ sont

$$u_2 + u_1 = 0, \quad 2u_2 + 2u_1 = 0$$

alors qu'à partir de $n = 2$, la récurrence force la valeur de $u_{n+3}$. Autrement dit, les valeurs $u_0, u_1, u_3, u_4$ sont libres, et les autres s'en déduisent. La dimension de l'espace des solutions est donc 4 pour cette récurrence d'ordre 3.

Cette dimension dépend des racines entières du coefficient de tête $p_r$. Si

$$H := \{h \in \mathbb{N} \mid p_r(h) = 0\} \quad \text{et} \quad J := \{0, \dots, r-1\} \cup \{h+r \mid h \in H\} \quad (15.6)$$

alors pour tout $n \in \mathbb{N} \setminus J$, l'équation (15.5) détermine de manière unique $u_n$ en fonction des valeurs précédentes. Les espaces vectoriels

$$\Lambda_K := \{(u_0, \dots, u_n) \in \mathbb{K}^{n+1} \mid (15.5) \text{ est satisfaite pour } n \in \{0, \dots, K-r\} \} \quad (15.7)$$

vérifient donc $\dim \Lambda_K = \dim \Lambda_{K-1}$ lorsque $K \notin J$.

Une façon simple de calculer la dimension de l'espace des solutions de la récurrence (15.5) ainsi qu'une base est donc de calculer la dimension et une base de $\Lambda_K$ pour $K = \max J$. Les valeurs de $u_N$ pour $N$ plus grand peuvent ensuite être calculées efficacement pour des grandes valeurs de $N$ par les Algorithmes 15.1 et 15.2 en observant que la récurrence obtenue en remplaçant $n$ par $n - \max J$ dans la récurrence (15.5) est régulière.

L'espace $\Lambda_K$ est défini par $K - r + 1$ équations linéaires (données par la récurrence) en $K + 1$ inconnues (les valeurs de la suite), et les équations obtenues pour $n \notin J$ sont linéairement indépendantes. Cette discussion permet de conclure :

**Proposition 15.5** La dimension de l'espace des solutions de la récurrence (15.5) est comprise entre $r$ et $r + \text{card}(H)$, où $H := \{h \in \mathbb{N} \mid p_r(h) = 0\}$.

Un calcul d'algèbre linéaire en dimension $K := \max J$ permet donc d'obtenir les K premières valeurs d'une base de l'espace des solutions. Dans de nombreux cas, la valeur de $K$ n'est pas très grande, comme dans l'exemple ci-dessus, et cette approche est suffisante. Mais en général, le maximum de $J$ a une valeur exponentielle en la taille binaire de la récurrence.

**Exemple 15.2** Pour la récurrence

$$(n-100)(n-200)u_{n+1} - u_n = 0, \quad n \in \mathbb{N},$$
la valeur de $\max J$ est 201, alors que le plus grand des coefficients de la récurrence (20000) a une taille binaire de seulement 14 bits. En remplaçant 100 et 200 par N et $2N$ avec N grand, ces estimations deviennent ~ 2log N pour la taille binaire de la récurrence, alors que $\max J$ croît comme $2N^2 + 1$. L'approche suggérée ci-dessus mène alors à la résolution d'un système linéaire de taille exponentielle par rapport à la taille binaire de la récurrence. Dans cet exemple, on peut procéder autrement : les espaces $\Lambda_N$ définis en (15.7) sont de dimension 1 pour $N \in \{0,...,100\}$, puis l'équation $-u_{100} = 0$ obtenue pour $n = 100$ force $\Lambda_{101}$ a être à nouveau de dimension 1, avec un générateur tel que $u_0 = \cdots = u_{100} = 0$ et $u_{101} = 1$. Les espaces suivants $\Lambda_N$ pour $N \in \{102,...,200\}$ sont encore de dimension 1, puis à nouveau l'équation $-u_{200} = 0$ force l'espace $\Lambda_{201}$ à être engendré par la suite $u_0 = \cdots = u_{200} = 0$, $u_{201} = 1$, qui fournit la base de l'espace des suites solution de la récurrence. Pour la déterminer, il a donc suffi de résoudre deux systèmes linéaires de petite dimension en $n = 100, 200$. Cette idée est développée dans la section suivante pour faire chuter la complexité du calcul d'une telle base.

## Conditions initiales généralisées

**Définition 15.2** On appellera *conditions initiales généralisées* de la récurrence

$$p_r(n)u_{n+r} + \cdots + p_0(n)u_n = q(n) \quad (15.8)$$

pour le support $\{0,...,N\}$ les valeurs $u_i$ pour $i \in J \cap \{0,...,N\}$, où J est défini par (15.6). Lorsque le support est $\mathbb{N}$, il n'est pas mentionné.

**Exemple 15.3** Pour une récurrence régulière, les conditions initiales généralisées sont les conditions initiales habituelles $u_0, \dots, u_{r-1}$.

**Exemple 15.4** Pour la récurrence de l'Exemple 15.2, les conditions initiales généralisées sont $u_0, u_{101}, u_{201}$ pour les supports contenant $\{0,...,200\}$.

**Exemple 15.5** Pour la récurrence de l'Exemple 15.1, 0,1,3,4 sont les indices des conditions initiales généralisées pour les supports contenant $\{0,...,4\}$.

Ces conditions initiales généralisées fournissent un test effectif d'égalité de suites, d'après le Lemme 14.2. Un autre intérêt immédiat est qu'elles permettent un calcul des éléments de la suite par l'Algorithme 15.3. À titre de comparaison, le calcul des N premiers coefficients est donné par une variante de la Proposition 15.1, de même preuve.

**Proposition 15.6** Étant donnée la récurrence (15.5) et des conditions initiales généralisées, l'Algorithme 15.3 calcule les N premiers éléments de la suite solution
**Entrée** Une récurrence de la forme (15.8), un entier $N \ge r-1$, des conditions initiales généralisées $u_0, \dots, u_{r-1}, u_{j_1}, \dots, u_{j_M}$ pour le support $\{0, \dots, N\}$.

**Sortie** Les termes $u_0, \dots, u_N$ de la solution.

1. Calculer les évaluations des coefficients $p_0(n), \dots, p_r(n), q(n)$ en $n = 0, \dots, N_r$.

2. Pour $n = 0, \dots, N$ faire :

a. Si $n$ est l'indice d'une condition initiale généralisée alors $u_n$ prend la valeur prescrite par cette condition initiale.

b. Sinon, $u_n$ est donné par (15.8) calculé avec les valeurs de l'étape (1) et les valeurs précédentes des $u_i$.

3. Renvoyer $u_0, \dots, u_N$.

**Algorithme 15.3** – Calcul naïf pour les récurrences singulières.

en

$$O\left(r \max(N, d) \frac{M(d)}{d}\right)$$

opérations arithmétiques, où $d$ est une borne sur le degré des polynômes coefficients.

Si les coefficients de ces polynômes sont des entiers de taille bornée par $\ell$ et que les conditions initiales généralisées sont des rationnels de numérateurs et dénominateurs bornés par $O(dN \log N + \ell N + N \log r)$, alors la complexité binaire de ce calcul est

$$O(rN^2 M_{\mathbb{Z}}(d \log N + \ell + \log r)),$$

la taille du résultat étant $O(N^2(d \log N + \ell + \log r))$.

Dans le cas de récurrences homogènes, les algorithmes rapides par « pas de bébés, pas de géants » et par scindage binaire s'étendent aussi au cas singulier, avec presque les mêmes complexités.

**Proposition 15.7** Étant données la récurrence linéaire (15.5) et des conditions initiales généralisées, l'Algorithme 15.4 calcule le N-ième terme d'une suite solution en

$$O\left(r^{\theta} M\left(\sqrt{d(M+1)N}\right) \log(dN)\right)$$

opérations arithmétiques, où $d$ est une borne sur le degré des coefficients de la récurrence, $M \le d$ est le nombre de racines entières positives de $p_r$ et $\theta$ un exposant réalisable pour la complexité du produit de matrices.

Si les coefficients des polynômes $p_i$ sont entiers de taille bornée par $\ell$ et les conditions initiales généralisées sont des entiers de taille au plus $O(dN \log N + \ell N + \dots)$.
**Entrée** Une récurrence de la forme (15.5), un entier $N \ge r-1$, des conditions initiales généralisées $u_0, \dots, u_{r-1}, u_{j_1}, \dots, u_{j_M}$ pour le support $\{0, \dots, N\}$.

**Sortie** Les termes $u_{N-r+1}, \dots, u_N$ de la solution.

1. Si $N = r-1$ renvoyer $u_0, \dots, u_{r-1}$.

2. Si $N = j_k$ pour un certain $k$, calculer récursivement
   $u_{j_k-r}, \dots, u_{j_k-1}$ par ce même algorithme;
   renvoyer $u_{j_k-r+1}, \dots, u_{j_k}$.

3. Soit $m = \max J \cap \{0, \dots, N\}$.

a. Calculer récursivement $u_{m-r+1}, \dots, u_m$.

b. Appliquer l'un des Algorithmes 15.1 ou 15.2 pour calculer $u_{N-r+1}, \dots, u_N$ avec la récurrence obtenue en changeant $n$ en $n-m+r-1$ dans (15.5), et les conditions initiales $u_{m-r+1}, \dots, u_m$.

c. Renvoyer le résultat.

**Algorithme 15.4** – Calcul rapide pour les récurrences singulières homogènes.

$N \log r)$, alors ce calcul utilise

$$O(M_{\mathbb{Z}}(dN \log N + \ell N + N \log r) (r^{\theta} \log N + r(M+1)))$$

opérations binaires.

**Démonstration.** D'abord l'algorithme est correct par récurrence sur $N$: pour $N = r-1$, c'est évident; pour $N$ un élément de $J$, la valeur renvoyée est correcte si elle l'est pour $N-1$; et pour $N \notin J$, les valeurs $u_{m-r+1}, \dots, u_m$ sont correctes par hypothèse de récurrence, et les conditions sont alors remplies pour invoquer les Algorithmes 15.1 ou 15.2 sur une récurrence régulière dans l'intervalle demandé.

Le résultat de complexité arithmétique est obtenu en invoquant l'Algorithmme 15.1 pour chaque $j_k - 1$ avec $j_k < m$ et pour $m$, c'est-à-dire au plus $1+M$ fois, où $M = \text{card} H$. À chaque appel, la récurrence obtenue en décalant les indices $a$ des coefficients dont le degré $d$ est inchangé. Le Théorème 15.2 donne une complexité arithmétique en $O(r^\theta M(\sqrt{dK}) \log(dK))$ pour une étape calculant une différence d'indices $K$.

La somme de ces complexités est donc de la forme

$$O(r^{\theta} (M(\sqrt{dK_1}) \log(dK_1) + \dots + M(\sqrt{dK_m}) \log(dK_m))) = \\ O(r^{\theta} \log(dN) (M(\sqrt{dK_1}) + \dots + M(\sqrt{dK_m})) = O(r^{\theta} \log(dN) M(\sqrt{dK_1} + \dots + \sqrt{dK_m}))$$

avec $dK_1 + \dots + dK_m \le dN$. La fonction racine carrée n'est pas convexe, mais l'inégalité

$$\sum_{i=1}^{n} \sqrt{x_i} \le \sqrt{n} \sqrt{\sum_{i=1}^{n} x_i}$$
se prouve facilement en développant

$$ \sum_{i=1}^{n} \left( \sqrt{x_i} - \frac{1}{n} \sum_{i=1}^{n} \sqrt{x_i} \right)^2 \geq 0 $$

et permet de borner la complexité par

$$ O(r^{\theta}M(\sqrt{mdN})\log(dN)). $$

De la même manière, le résultat de complexité binaire est obtenu en invoquant l'Algorithme 15.2 pour chaque $j_k - 1$ avec $j_k < m$ et pour $m$. À l'étape $k$, les polynômes coefficients de la récurrence décalée ont un degré inchangé, et des coefficients avec une taille majorée par $\ell + d \log j_{k-1}$ (avec $j_0 = r-1$). La complexité binaire de l'ensemble des produits de matrices dans l'Algorithme 15.2 est alors donnée par le Théorème 15.3 comme

$$ O(r^{\theta}M_{\mathbb{Z}}(dK_k \log K_k + (\ell + d \log j_{k-1})K_k + K_k \log r) \log K_k). $$

En utilisant $j_{k-1} \le N$ et la sous-additivité de la fonction $M_{\mathbb{Z}}$, ces complexités se somment en

$$ O(r^{\theta}M_{\mathbb{Z}}(dN \log N + \ell N + N \log r) \log N). $$

Il reste à considérer les produits matrice-vecteur de l'Algorithme 15.2, qui sont au plus au nombre de $M+1$, et font intervenir des entiers de taille bornée par $O(dN \log N + \ell N + N \log r)$, ce qui conclut la preuve. ■

Ces résultats, en $\tilde{O}(\sqrt{N})$ pour la complexité arithmétique, et en $\tilde{O}(N)$ pour la complexité binaire, généralisent ceux du cas régulier, et sont à comparer à ceux de la Proposition 15.6 pour l'algorithme naïf.

## Homogénéisation des récurrences inhomogènes

Lorsque la récurrence inhomogène (15.8) a un membre droit $q(n)$ lui-même polynomial, alors ses solutions sont aussi solutions de la récurrence homogène

$$ q(n) (p_r(n+1)u_{n+1+r} + \dots + p_0(n+1)u_{n+1}) - q(n+1) (p_r(n)u_{n+r} + \dots + p_0(n)u_n) = 0. \quad (15.9) $$

Cette récurrence a d'autres solutions que celles de la récurrence (15.8). Si $J$ est l'ensemble des indices des conditions initiales généralisées pour l'équation inhomogène, alors celui de la récurrence homogénéisée est

$$ J' := J \cup \{r\} \cup \{h+r+1 \mid q(h) = 0\}. $$

Les valeurs en $J' \setminus J$ permettant de retrouver exactement les mêmes solutions sont calculées facilement à partir de la récurrence inhomogène de départ.

**Exemple 15.6** Voici une variante de l'Exemple 15.2 :

$$ (n-100)(n-200)u_{n+1} - u_n = n-50, \quad n \in \mathbb{N}. $$
**Entrée** Une récurrence de la forme (15.8), un entier $N \ge r-1$, des conditions initiales généralisées $u_0, \dots, u_{r-1}, u_{j_1}, \dots, u_{j_M}$ pour le support $\{0, \dots, N\}$.

**Sortie** Les termes $u_{N-r}, \dots, u_N$ de la solution.

1. Si $N = j_k$ pour un certain $k$, calculer récursivement $u_{j_k-r}, \dots, u_{j_k-1}$ par ce même algorithme; renvoyer $u_{j_k-r+1}, \dots, u_{j_k}$.

2. Sinon, si $N \in J' := \{r\} \cup \{h \in \mathbb{N} \mid h > r \text{ and } q(h-r-1) = 0\}$ alors calculer récursivement $u_{N-r-1}, \dots, u_{j_k-1}$ par ce même algorithme (avec la convention $u_{-1} = 0$); calculer $u_N$ par (15.8); renvoyer $u_{N-r}, \dots, u_N$.

3. Soit $m = \max(J \cup J') \cap \{0, \dots, N\}$.

a. Calculer récursivement $u_{m-r}, \dots, u_m$.

b. Appliquer l'un des Algorithmes 15.1 ou 15.2 pour calculer $u_{N-r+1}, \dots, u_N$ avec la récurrence obtenue en changeant $n$ en $n-m+r$ dans la récurrence homogénéisée (15.9), et les conditions initiales $u_{m-r}, \dots, u_m$.

c. Renvoyer le résultat.

**Algorithme 15.5** – Calcul rapide pour les récurrences singulières inhomogènes.

L'équation homogénéisée correspondante est

$$ (n-50)(n-99)(n-199)u_{n+2} - (n-50 + (n-49)(n-100)(n-200))u_{n+1} + (n-49)u_n = 0. $$

Les indices des conditions initiales généralisées pour cette nouvelle récurrence sont $\{0, 1, 52, 101, 201\}$, c'est-à-dire ceux de la partie homogène, ainsi que $\{1, 52\}$. En ces indices, les valeurs fournies par la récurrence de départ sont :

$$ u_1 = \frac{u_0 - 50}{20000}, \quad u_{52} = \frac{u_{51} + 1}{49 \times 149}, $$

et la valeur de $u_{51}$ peut quant à elle être obtenue par scindage binaire sur l'équation homogénéisée.

Procédant ainsi, l'Algorithmme 15.5 permet de résoudre les équations inhomogènes dans la même complexité que précédemment.

**Théorème 15.8** Étant données la récurrence linéaire inhomogène (15.8) et des conditions initiales généralisées, l'Algorithmme 15.5 calcule le N-ième terme d'une suite solution en

$$ O\left(r^0 M\left(\sqrt{d(M+1)N}\right) \log(dN)\right) $$

opérations arithmétiques, où $d$ est une borne sur le degré des coefficients de la récurrence, $M \le 2d$ est le nombre de racines entières positives de $q$ et $p_r$, et $\theta$ un
exposant réalisable pour la complexité du produit de matrices.

Si les coefficients des polynômes $p_i$ et de $q$ sont entiers de taille bornée par $\ell$ et les conditions initiales généralisées sont des entiers de taille au plus $O(dN \log N + \ell N + N \log r)$, alors ce calcul utilise

$$O(M_Z(dN \log N + \ell N + N \log r)(r^\theta \log N + r(M+1)))$$

opérations binaires.

*Démonstration.* L'algorithme est correct par récurrence sur $N$, la seule différence avec la preuve de la Proposition 15.7 provenant des points de J' qui sont soit des conditions initiales généralisées, soit des valeurs où la récurrence inhomogène permet le calcul.

L'analyse de complexité est exactement la même que pour le cas homogène. ■

## Exercices

**Exercice 15.2** — Calcul rapide de factorielle et de coefficients binomiaux centraux. Cet exercice montre comment calculer certaines suites récurrentes linéaires plus vite que par la méthode de scindage binaire. Soit $N \in \mathbb{N}$ et soit $Q = \sum_{i=0}^{2N} q_i X^i \in \mathbb{Z}[X]$ le polynôme $Q(X) = (1+X)^{2N}$.

1. Montrer que $q_N$ peut être calculé en utilisant uniquement des additions d'entiers du triangle de Pascal, c'est-à-dire l'identité suivante sur les coefficients du binôme :

$$\binom{n}{k} = \binom{n-1}{k-1} + \binom{n-1}{k}, \quad k \ge 1, n \ge 1.$$

Quelle est la complexité binaire de cet algorithme ?

On admet que le calcul de tous les nombres premiers inférieurs à N peut être effectué en $O(N \log N / \log \log N)$ opérations binaires, qu'il existe au plus 3N/ log N tels nombres premiers et que la multiplicité avec laquelle apparaît le nombre premier p dans la factorisation de N! vaut

$$\mathrm{ind}(p, N) = \sum_{i=1}^{\infty} \lfloor \frac{N}{p^i} \rfloor,$$

où la notation $\lfloor x \rfloor$ représente la partie entière de $x$.

2. Montrer que le calcul de $\mathrm{ind}(p, N)$ peut être effectué en

$$O(M_Z(\log N) \log N)$$

opérations binaires.

3. Montrer que la décomposition en facteurs premiers de N! peut être effectuée en $O(M_Z(N) \log N)$ opérations binaires, ainsi que celle de $q_N$.

4. Montrer que N! et $q_N$ peuvent alors être reconstruits en un nombre d'opérations binaires bornés respectivement par $O(M_Z(N \log N) \log \log N)$ et $O(M_Z(N) \log N)$.
5. Comparer avec la complexité du scindage binaire pour ces calculs.

**Exercice 15.3** Soit $N \in \mathbb{N}$ et soit $P = \sum_{i=0}^{2N} p_i X^i \in \mathbb{Z}[X]$ le polynôme $P(X) = (1+X+X^2)^N$.

1. Montrer que $O(M(N))$ opérations binaires suffisent pour déterminer la parité de tous les coefficients de P.

Indication : un entier $n$ est pair si et seulement si $n = 0$ dans $\mathbb{K} = \mathbb{Z}/2\mathbb{Z}$.

2. Montrer que P vérifie une équation différentielle linéaire d'ordre 1 à coefficients polynomiaux. En déduire que les $p_i$ suivent une récurrence d'ordre 2 que l'on précisera.

3. Donner un algorithme qui calcule $p_N$ en $O(M(N \log N) \log N)$ opérations binaires.

## Notes

L'algorithme par pas de bébés et pas de géants a été introduit par Strassen [Str77] et généralisé par les frères Chudnovsky [CC88] au problème du calcul d'un terme d'une suite récurrente linéaire à coefficients polynomiaux, dans un article où ils présentent également le scindage binaire et le doublement de précision.

L'algorithme de factorisation déterministe en Section 15.2 est basé également sur l'article de Strassen [Str77]. Cet algorithme peut à son tour être un peu amélioré [BGS07; CH14]. On ne connaît pas d'algorithme déterministe essentiellement meilleur ; l'existence d'un tel algorithme est un grand problème ouvert.

Une application cryptologique du calcul du N-ième terme d'une suite polynomialement récurrente au problème du comptage de points sur une courbe hyperelliptique sur un corps fini a été développée par Bostan, Gaudry et Schost [BGS07]. Le point de départ de cette application est l'exercice 15.1, inspiré de l'article de Flajolet et Salvy de 1997 [FS97, Pb. 4]. D'autres applications sont données au Chapitre 17.

La partie de l'Exercice 15.2 consacrée au calcul de $N!$ est due à P. Borwein [Bor85].

Une conséquence des résultats de ce chapitre est que les N premiers coefficients d'une série différentiellement finie peuvent être calculés en complexité arithmétique $O(N)$ et en complexité binaire quasi-linéaire en $N^2$, en déroulant la récurrence linéaire à coefficients polynomiaux vérifiée par les coefficients de la série. Ces complexités sont quasi-optimales par rapport à la taille de la sortie.

En exploitant plus finement les mêmes récurrences sur les coefficients, on peut calculer le N-ième terme $f_N$ d'une série différentiellement finie F plus vite que l'ensemble des N + 1 premiers coefficients $f_0, \dots, f_N$. Plus exactement, la technique des pas de bébés et pas de géants permet de calculer $f_N$ en complexité arithmétique quasi-linéaire en $\sqrt{N}$ et la technique du scindage binaire de calculer $f_N$ en complexité binaire quasi-linéaire en N. Cette dernière complexité est bien quasi-optimale en la taille de la sortie, mais celle en $\sqrt{N}$ ne l'est pas. Un grand problème ouvert est de savoir si l'on peut calculer $f_N$ plus rapidement, ou si l'exposant $1/2$ est une barrière intrinsèque. La réponse demeure à ce jour inconnue même dans des cas très particuliers, comme par exemple pour le calcul de $(\frac{p-1}{2})!$ mod $p$, pour un nombre premier $p$.

Un autre cas particulier important concerne les séries algébriques. Peut-on calculer le N-ième coefficient d'une série algébrique plus rapidement que celui d'une série
différentiellement finie quelconque? La question vise les deux modèles de complexité, arithmétique et binaire. Par exemple, un tel coefficient a typiquement une taille binaire O(N), alors que le meilleur algorithme actuellement connu passe par le calcul intermédiaire d'une factorielle de matrices, de taille binaire N log N. Dans le modèle arithmétique, aucun algorithme de complexité plus basse que $\sqrt{N}$ n'est actuellement connu, sauf dans le cas des corps finis, où une complexité en log N a été récemment obtenue [BCD16].

Bibliographie

BCD16 BOSTAN, Alin, Gilles CHRISTOL et Philippe DUMAS (2016). « Fast Computation of the Nth Term of an Algebraic Series over a Finite Prime Field ». In: Proceedings of the ACM on International Symposium on Symbolic and Algebraic Computation. ISSAC'16. ACM, p. 119-126.

BGS07 BOSTAN, Alin, Pierrick GAUDRY et Éric SCHOST (2007). « Linear recurrences with polynomial coefficients and application to integer factorization and Cartier-Manin operator ». In: SIAM Journal on Computing, vol. 36, n°6, p. 1777-1806.

Bor85 BORWEIN, Peter B. (1985). « On the complexity of calculating factorials ». In: Journal of Algorithms, vol. 6, n°3, p. 376-380.

CC88 CHUDNOVSKY, D. V. et G. V. CHUDNOVSKY (1988). « Approximations and complex multiplication according to Ramanujan ». In: Proceedings of the Centenary Conference, University of Illinois at Urbana-Champaign, June 1-5, 1987. Éd. par G. E. ANDREWS, R. A. ASKEY, B. C. BERNDT, K. G. RAMANATHAN et R. A. RANKIN. Academic Press, p. 375-472.

CH14 COSTA, Edgar et David HARVEY (2014). « Faster deterministic integer factorization ». In: Mathematics of Computation, vol. 83, n°285, p. 339-345.

FS97 FLAJOLET, Philippe et Bruno SALVY (1997). « The SIGSAM challenges : symbolic asymptotics in practice ». In: ACM SIGSAM Bulletin, vol. 31, n°4, p. 36-47.

Str77 STRASSEN, V. (1976/77). « Einige Resultate über Berechnungskomplexität ». In: Jahresbericht der Deutschen Mathematiker-Vereinigung, vol. 78, n°1, p. 1-8.
<figure><img src="image_49.png" /></figure>
# 16. Résolution de récurrences linéaires

## Résumé

Les solutions polynomiales ou rationnelles de récurrences linéaires peuvent avoir une taille exponentiellement grande en la taille de l'équation. Néanmoins, elles s'obtiennent efficacement en utilisant la structure de leurs espaces de solutions, en évitant de développer les polynômes concernés, et en exploitant le calcul rapide de factorielles de matrices.

Les chapitres précédents montrent que de nombreuses informations concernant les solutions d'équations différentielles ou de récurrences linéaires peuvent être calculées directement à partir de l'équation sans passer par une « résolution » de celle-ci. Le cas particulier des solutions polynomiales et rationnelles mérite une attention spéciale : ces solutions, lorsqu'elles existent, permettent de nombreux calculs efficaces. D'autre part, décider de leur existence est une brique de base dans les algorithmes de sommation et d'intégration de la Partie VI.

Le cas le plus simple est celui d'équations *homogènes*, c'est-à-dire de la forme

$$p_r(n)u_{n+r} + \cdots + p_0(n)u_n = 0, \quad (16.1)$$

où les coefficients $p_i$ sont des polynômes donnés à coefficients dans un corps $\mathbb{K}$ de caractéristique nulle. La question la plus générale que nous abordons est celle du calcul des solutions rationnelles, qui se ramène à celle des solutions polynomiales, elle-même se réduisant à la question des solutions à support fini (voir ci-dessous).

Un problème légèrement plus difficile, y compris du point de vue de la complexité, est celui des équations *inhomogènes*, de la forme

$$p_r(n)u_{n+r} + \cdots + p_0(n)u_n = q(n), \quad (16.2)$$
c'est-à-dire du même type, mais avec un membre droit qui n'est plus simplement 0, mais une suite $q(n)$ qui peut être à support fini, un polynôme ou une fraction rationnelle. Enfin, des algorithmes de sommation définie traités dans la Partie VI de cet ouvrage mènent à des équations inhomogènes paramétrées, où le membre droit est maintenant de la forme

$$q(n) = \lambda_1 q_1(n) + \cdots + \lambda_k q_k(n),$$

où les $q_i$ sont des fractions rationnelles données, et les coefficients $\lambda_i$ sont inconnus. Il s'agit alors de déterminer l'existence et la valeur des coefficients $\lambda_i$ permettant l'existence d'une solution rationnelle de l'équation (16.2).

## 16.1 Suites solution

### La forme d'une suite polynomialement récurrente

Dans le cas fréquent où les coefficients de la récurrence (16.1) sont des polynômes à coefficients entiers, les tailles des éléments de la suite croissent avec leur indice dans le cas général. Ces questions de taille ont été traitées au Chapitre 15. Il en ressort que la taille en bits du N-ième élément de la suite croît en $O(dN \log N)$ où $d$ est le degré des coefficients de la récurrence, et que la taille totale des N premiers éléments est bornée par $O(dN^2 \log N)$.

Lorsque $N$ est grand, ces tailles ont un impact sur les complexités, et il sera utile de concevoir des algorithmes qui opèrent sur ces suites sans nécessairement en représenter explicitement tous les éléments. En effet, les indices de suites pertinents dans les algorithmes de ce chapitre proviennent de récurrences rendues singulières par l'existence de solutions entières de polynômes. Ces solutions peuvent être exponentiellement grandes en la taille binaire de la récurrence (le polynôme $n-N$ est un exemple frappant, voir aussi l'Exemple 15.2).

### Calcul d'une base des solutions

Les Propositions 15.5, 15.7 et le Théorème 15.8 montrent qu'avec des conditions initiales généralisées en nombre au plus $r+\deg p_r+(\deg q)$, il est possible de calculer, et ce efficacement, tout terme d'une suite solution d'une des récurrences (16.1) et (16.2). Un tel ensemble de conditions initiales généralisées fournit donc une structure de données pour représenter une solution. Tous les choix de conditions initiales généralisées ne correspondent pas nécessairement à une solution, mais il est possible de calculer une base des solutions dans cette représentation, tout en évitant des calculs dont la complexité arithmétique serait au moins linéaire en la plus grande racine entière du polynôme coefficient de tête $p_r$, ou du polynôme $q$ dans le cas inhomogène.

Le principe du calcul est d'exploiter l'Algorithme 15.5 de calcule rapide par « pas de bébés, pas de géants » ou par scindage binaire avec des conditions initiales généralisées arbitraires pour en déduire les contraintes que l'évaluation de la récurrence aux zéros du polynôme $p_r$ et à ceux du polynôme $q$ dans le cas inhomogène fait porter sur ces conditions.
**Entrée** Une récurrence de la forme (15.5).

**Sortie** Les conditions initiales généralisées d'une base de ses solutions.

1. Calculer l'ensemble H des racines entières positives du coefficient de tête $p_r$.

2. En déduire $J := \{0, \dots, r-1\} \cup \{h+r \mid h \in H\}$.

3. Utiliser l'Algorithme 15.4 pour calculer des matrices $A_1, \dots, A_{\text{card}H}$ de dimensions $r \times \text{card}J$, où $A_h$ avec $h \in H$ contient pour colonnes les valeurs de $(u_h, u_{h+1}, \dots, u_{h+r-1})$ pour les conditions initiales généralisées fournies par les colonnes successives de la matrice identité $I_{\text{card}J}$.

4. Former une matrice M de dimensions $\text{card}H \times \text{card}J$ avec pour lignes

$$ (p_0(h), \dots, p_{r-1}(h)) \cdot A_h, \quad h \in H. $$

5. Renvoyer une base du noyau de M et l'ensemble H.

**Algorithme 16.1** – Résolution de récurrences homogènes singulières.

**Exemple 16.1** La récurrence

$$ (n-4)(n-8)u_{n+2} - u_{n+1} + u_n = 0, \quad n \in \mathbb{N}, $$

a $\{0,1,6,10\}$ comme ensemble d'indices pour ses conditions initiales généralisées. Les choix de valeurs $(a,b,c,d)$ pour ces conditions initiales ne sont pas tous possibles, mais même avec des valeurs qui ne correspondent pas à des solutions, l'Algorithme 15.4 peut être exécuté, et permet de calculer

$$ u_4 = \frac{5}{2016}a - \frac{13}{2016}b, \quad u_5 = \frac{1}{1260}a + \frac{1}{126}b, $$

$$ u_8 = -\frac{1}{15120}a - \frac{1}{1512}b + \frac{1}{3}c, \quad u_9 = \frac{1}{9072}a + \frac{5}{4536}b - \frac{2}{9}c. $$

Les équations $u_4 - u_5 = 0, u_8 - u_9 = 0$ obtenues aux racines de $(n-4)(n-8)$ sont les seules contraintes qui pèsent sur les conditions initiales, qui sont donc fournies par le noyau de la matrice

$$ \begin{pmatrix} \frac{17}{10080} & -\frac{29}{2016} & 0 & 0 \\ -\frac{1}{5670} & -\frac{1}{567} & 5/9 & 0 \end{pmatrix}, $$

dont une base est donnée par $(1450, 170, 1, 0), (0, 0, 0, 1)$.

L'Algorithme 16.1 explicite cette idée dans le cas homogène. Lorsque la récurrence est régulière, l'étape (3) ne fait rien, la matrice de l'étape (4) est vide et la dernière étape renvoie la matrice identité $I_r$.
## Récurrences inhomogènes

La récurrence linéaire inhomogène

$$p_r(n)u_{n+r} + \cdots + p_0(n)u_n = q(n), \quad n \in \mathbb{N} \qquad (16.3)$$

a pour ensemble de solutions un espace affine : deux solutions diffèrent par un élément de l'espace vectoriel des suites solutions de son membre gauche. On peut voir ce problème comme un problème d'algèbre linéaire en dimension infinie : le vecteur $(u_n)$ est multiplié par une matrice bande dont les coefficients de $(n, n)$ à $(n, n+r)$ sont fournis par les évaluations des polynômes $p_i$ en $n$; ce produit doit être égal au vecteur $(q(n))$.

Le cas régulier où $p_r$ n'a pas de racines entières signifie qu'il n'y a pas de 0 sur la diagonale de la matrice, et fournit donc une solution pour chaque condition initiale $(u_0, \dots, u_{r-1})$, que l'on obtient simplement en déroulant la récurrence.

Dans le cas général, il suffit de trouver des conditions initiales généralisées permettant à une solution d'exister jusqu'au plus grand zéro entier de $p_r$, puisqu'ensuite cette solution se prolonge de manière unique comme dans le cas régulier, et toutes les solutions sont obtenues en y ajoutant les solutions de l'équation homogène associée.

### Exemple 16.2 La partie homogène de la récurrence

$$(n-4)(n-2)u_{n+1} - nu_n = 1$$

a un espace de solutions de dimension 2, engendré par les suites de conditions initiales généralisées $(u_0 = 1, u_2 = 0, u_4 = 0)$ et $(u_0 = 0, u_2 = 0, u_4 = 1)$. Par contre, la récurrence n'a pas de solution : en partant avec $u_0$ arbitraire, la récurrence produit

$$u_1 = \frac{1}{8}, \quad u_2 = \frac{3}{8}$$

puis la relation $-2u_2 - 1 = 0$, qui est incompatible avec ces équations.

### Exemple 16.3 La récurrence

$$(n-4)(n-2)u_{n+2} + nu_{n+1} - u_n = 1$$

a une partie homogène admettant un espace de solutions de dimension 2. Les indices des conditions initiales généralisées sont $(0, 1, 4, 6)$. Les premières évaluations de la récurrence fournissent les équations

$$8u_2 - u_0 = 1, \quad 3u_3 + u_2 - u_1 = 1, \quad 2u_3 - u_2 = 1, \quad -u_5 + 3u_4 - u_3 = 1, \quad 4u_5 - u_4 = 1,$$

puis des équations qui déterminent uniquement $u_7, u_8, \dots$ en fonction des précédents.

À ce stade, le problème est donc réduit à la résolution d'un système linéaire de dimension 6, et il reste à voir comment exploiter la structure supplémentaire induite par la récurrence. Pour cela on peut tout d'abord utiliser la partie homogène
**Entrée** Une récurrence de la forme (16.3).

**Sortie** Les conditions initiales généralisées d'une base de ses solutions.

1. Calculer la matrice M de l'Algorithmme 16.1.

2. Utiliser l'Algorithmme 15.5 avec des conditions initiales généralisées nulles pour calculer un vecteur B contenant les valeurs de la solution correspondante aux éléments de H.

3. Renvoyer une base des solutions de $MY + B = 0$ et l'ensemble H.

**Algorithme 16.2** – Résolution de récurrences homogènes singulières.

de la récurrence avec des conditions initiales généralisées arbitraires en 0,1,4,6 ce qui donne pour $n = 1,3$ les équations $(u_2, u_3) = (u_0/8, -u_0/24 + u_1/3)$ et $(u_4, u_5) = (u_4, u_0/24 - u_1/3 + 3u_4)$. Ensuite, on déroule la récurrence homogénéisée

$$ (n-3)(n-1)u_{n+3} + (n+1-(n-4)(n-2))u_{n+2} - (n+1)u_{n+1} + u_n = 0 $$

avec des conditions initiales généralisées nulles et la nouvelle condition en $n = 2$ fournie par l'équation inhomogène en $n = 0$, à savoir $u_2 = 1/8$. La récurrence homogène permet de calculer facilement (et efficacement même si ces indices étaient bien plus grands) les valeurs $(u_1, u_2, u_3) = (0, 1/8, 7/24)$, puis $(u_3, u_4, u_5) = (7/24, 0, -31/24)$. Par addition, le cas général vaut donc

$$ (u_2, u_3) = (1/8+u_0/8, 7/24-u_0/24+u_1/3), \quad (u_4, u_5) = (u_4, -31/24+u_0/24-u_1/3+3u_4). $$

Il ne reste plus qu'à introduire les équations pour $n = 2$ et $4$ qui entraînent

$$ -13/24 - 5u_0/24 + 2u_1/3 = -37/6 + u_0/6 - 4u_1/3 + 11u_4 = 0, $$

d'où finalement les conditions initiales

$$ u_1 = 13/16 + 5u_0/16, \quad u_4 = 29/44 + u_0/44, \quad u_0 \text{ arbitraire.} $$

En résumé, la résolution d'un seul système linéaire de petite taille permet d'obtenir que la solution générale de l'équation a pour conditions initiales généralisées $(u_0, u_1, u_4, u_6) = (0, 13/16, 29/44, 0) + \lambda(1, 5/16, 1/44, 245/234344) + \mu(0, 0, 0, 1)$ avec $\lambda, \mu$ arbitraires.

L'Algorithmme 16.2 précise cette méthode dans le cas général.

## Racines entières

Sans information supplémentaire sur le corps $K$ où se trouvent les coefficients des polynômes $p_i$ de la récurrence, la recherche de racines entières peut poser des difficultés. Nous faisons donc l'hypothèse qu'un oracle les fournit sans compter ses opérations arithmétiques. Lorsque les coefficients sont entiers en revanche, les algorithmes de la Partie IV sur la factorisation peuvent être utilisés et donnent une complexité binaire
en $O(d^2\ell)$ pour trouver les racines entières d'un polynôme de degré $d$ à coefficients de taille bornée par $\ell$ (grâce au Théorème 21.6 et aux Propositions 21.16 page 383 et 21.22 page 390).

## Complexité

**Théorème 16.1** L'ensemble H des racines entières positives du coefficient $p_r$ (ainsi que celles de $q$ dans le cas inhomogène) étant donné, l'Algorithme 16.1 (resp. 16.2) calcule les conditions initiales généralisées d'une base de solutions de la récurrence (16.1) (resp. (16.3)) en

$$O\left(r^{\theta}\left(M\left(\sqrt{d(M+1)N}\log(dN) + M(1+M/r)\right)\right)\right)$$

opérations arithmétiques, où $d$ est une borne sur le degré des coefficients de la récurrence, $M \le d$ (resp. $M \le 2d$) est le nombre de racines entières positives de $p_r$ (resp. et de $q$) et $N$ est une borne sur le plus grand élément de H.

Si de plus les coefficients de la récurrence sont dans $\mathbb{Z}[n]$ avec des coefficients de taille au plus $\ell$, alors les conditions initiales généralisées peuvent être prises entières de taille $O(T)$ où $T = r(dN \log N + \ell N + N \log r)$, et ce calcul ainsi que celui de l'ensemble H (resp. et des racines de $q$) peut être effectué en

$$O\left((r+M)^2 M_{\mathbb{Z}}(T) \log(T) + r M d^2 M_{\mathbb{Z}}(\ell + \log N)\right)$$

opérations binaires.

**Démonstration.** D'après la Proposition 15.7, l'étape (3) de l'Algorithme 16.1 requiert $O(r^\theta M(\sqrt{d(M+1)N}\log(dN)))$ opérations arithmétiques pour calculer toutes les matrices de passage $r \times r$ nécessaires, qui sont indépendantes des conditions initiales. Il faut ensuite multiplier les $M$ matrices ainsi construites par $\text{card} J = r + M$ vecteurs de dimension $r$, que l'on peut grouper en $\text{card} J/r$ matrices $r \times r$, ce qui amène le coût de l'ensemble de ces produits à $O(r^\theta M(1 + M/r))$ opérations arithmétiques supplémentaires.

L'étape suivante demande l'évaluation de chacun des polynômes coefficients en les éléments de H, ce qui, par évaluation multi-points, coûte $O(rMM(d) \log d/d)$ et est donc négligeable par rapport aux opérations précédentes. Dans le cas inhomogène, le Théorème 15.8 donne le vecteur B en la complexité annoncée. Enfin, la résolution du système linéaire coûte $O((M+r)r^{\theta-1})$ opérations arithmétiques d'après le Théorème 8.6, ce qui est du même ordre que les multiplications précédentes.

En ce qui concerne la complexité binaire, le calcul de l'étape (3) est donné par $O(r^\theta M_{\mathbb{Z}}(dN \log N + \ell N + N \log r) \log N)$ pour obtenir toutes les matrices de passage. Les vecteurs $(u_h, \dots, u_{h+r-1})$ ont une taille binaire en $O(dN \log N + \ell N + N \log r)$ ce qui fait que le calcul des matrices $A_h$ est borné par celui des matrices de passage. Ces mêmes complexités apparaissent dans les étapes de l'Algorithme 15.5. Par les formules de Cramer et la borne d'Hadamard (Proposition 21.10), les vecteurs de la solution du système linéaire ont pour coordonnées des rationnels que l'on peut prendre de numérateur et dénominateur de taille $T = O(r(dN \log N + \ell N + N \log r))$.
**Entrée** Une récurrence linéaire régulière de la forme (16.1).

**Sortie** Une base des solutions à support fini, données par leurs conditions initiales.

1. Calculer $S := \{h \in \mathbb{N} \mid p_0(h) = 0\}$.

2. Si $S = \emptyset$, renvoyer $\emptyset$; sinon poser $N := \max S$.

3. Utiliser l'un des Algorithmes 15.1 ou 15.2 pour calculer la matrice de transition entre les indices $(0, \dots, r-1)$ et $(N+1, \dots, N+r)$.

4. Calculer une matrice C dont les colonnes forment une base du noyau de cette matrice.

5. Si le noyau est réduit à 0, renvoyer $\emptyset$; sinon renvoyer C.

**Algorithme 16.3** – Calcul d'une base de solutions à support fini, cas régulier.

Leur calcul peut-être effectué en utilisant de l'algèbre linéaire modulo des nombres premiers en nombre de cet ordre en $O(r^\theta(1+M/r)T)$ opérations binaires, puis par reconstruction par restes chinois en $O((r+M)^2M(T)\log T)$ opérations binaires, ce qui domine donc l'étape de construction des matrices de passage. ■

## 16.2 Solutions à support fini

Les algorithmes du reste de ce chapitre, pour calculer des solutions polynomiales ou rationnelles, se ramènent tous au calcul de solutions à support fini de récurrences linéaires.

**Définition 16.1** On appelle *support* d'une suite $(u_n)$, l'ensemble des $n \in \mathbb{N}$ tels que $u_n \neq 0$. Lorsque cet ensemble est fini, son plus grand élément s'appelle le *degré* de la suite, et on dit que le degré est infini sinon.

Les récurrences inhomogènes à coefficients polynomiaux comme (15.3) n'ont pas de solution à support fini : pour $n$ plus grand que le degré, le membre gauche s'annule, et donc le membre droit, un polynôme, devrait avoir une infinité de racines, ce qui le force à être nul. Il suffit donc de considérer le cas homogène.

**Proposition 16.2** Si $(u_n)$ est une solution à support fini de la récurrence (16.1), alors son degré est racine du polynôme $p_0$. Réciproquement, si $p_0$ a K racines dans $\mathbb{N}$, alors la récurrence possède un espace de solutions à support fini de dimension comprise entre 1 et K.

*Démonstration.* Pour le sens direct, il suffit d'évaluer la récurrence en $n = N$ où N est le degré, ce qui donne $p_0(N)u_N = 0$. La réciproque est similaire à la Proposition 15.5 : l'évaluation de la récurrence en une racine entière $h$ de $p_0$ ne contraint pas $u_h$, mais crée une contrainte linéaire sur les valeurs suivantes. Selon que cette contrainte est linéairement dépendante ou non des autres, elle permet ou non d'augmenter la dimension. ■
**Entrée** Une récurrence linéaire de la forme (16.1).

**Sortie** Une base des solutions à support fini, données par leurs conditions initiales généralisées.

1. Calculer $S := \{h \in \mathbb{N} | p_0(h) = 0\}$.

2. Si $S = \emptyset$, renvoyer $\emptyset$; sinon poser $N := \max S$.

3. Utiliser l'Algorithme 16.1 pour calculer une matrice B de dimension $\mathbf{J} \times k$ dont les colonnes sont les conditions initiales généralisées d'une base de solutions.

4. Utiliser l'Algorithme 15.4 pour calculer une matrice de dimensions $\mathbf{r} \times k$ dont les colonnes contiennent les valeurs aux indices $(N+1, \dots, N+r)$ de chacun des éléments de cette base.

5. Calculer une matrice C dont les colonnes forment une base du noyau de cette matrice.

6. Si le noyau est réduit à 0, renvoyer $\emptyset$; sinon renvoyer $\mathbf{B} \cdot \mathbf{C}$.

**Algorithme 16.4** – Calcul d'une base de solutions à support fini, cas singulier.

**Exemple 16.4** Voici un cas extrême où toutes les solutions sont à support fini :

$$u_{n+2} + (n-7)u_{n+1} + (n-10)(n-5)u_n = 0,$$

(exercice : le vérifier) mais changer le coefficient 7 en 6 par exemple fait chuter la dimension des solutions à support fini à 1, et le degré à 5.

Le principe de l'Algorithme 16.3 pour le cas de récurrences régulières est très simple : une fois trouvé $N$ la plus grande racine entière positive du coefficient $p_0$, il s'agit de trouver les conditions initiales de la récurrence qui permettent d'obtenir $u_{N+1} = \dots = u_{N+r} = 0$, ce qui donne l'ensemble des solutions à support fini dans ce cas. Dans le cas singulier, où le coefficient de tête de la récurrence a des zéros entiers qui peuvent à la fois augmenter la dimension de l'espace des solutions et ajouter des contraintes sur les conditions initiales, l'Algorithme 16.4, qui généralise le précédent, procède de façon similaire, à partir de l'Algorithme 16.1.

**Théorème 16.3** Les racines entières positives de $p_0$ et $p_r$ étant données, les Algorithmes 16.3 et 16.4 calculent une base de l'espace des solutions à support fini de la récurrence (16.1). Leur complexité arithmétique est

$$O(r^\theta M(\sqrt{dN}) \log(dN))$$

où $d$ est une borne sur le degré des coefficients et $N$ la plus grande racine entière positive de $p_0 p_r$.

Lorsque les coefficients des polynômes $p_i$ sont des entiers, la complexité binaire du calcul des racines entières positives de $p_0$ et $p_r$ ainsi que du reste de l'algorithme
est

$$O((r + M)^2 M_Z(T) \log(T) + r M d^2 M_Z(\ell + \log N)),$$

où $T = r(dN \log N + \ell N + N \log r)$ est alors une borne sur la taille des conditions initiales généralisées trouvées, $\ell$ étant une borne sur la taille des coefficients des polynômes $p_i$.

*Démonstration.* Les bornes de complexité des différentes étapes sont les mêmes que celles de l'Algorithme 16.1, la seule différence étant la borne N sur les indices. ■

R

En pratique, même dans le cas où les coefficients sont entiers, si N est grand, il vaut donc mieux commencer par tester l'existence d'une racine modulo un nombre premier de taille modérée en un nombre d'opérations quasi-linéaire en $\sqrt{N}$, puis, s'il en existe, les calculer en un nombre quasi-linéaire d'opérations en N. Si l'on souhaite calculer tous les coefficients, alors avec ces conditions initiales généralisées, l'Algorithme 15.5 les fournit en complexité quasi-linéaire en $N^2$, ce qui est quasi-optimal. Pour bien des opérations cependant, ce calcul n'est pas nécessaire et la récurrence ainsi que ses conditions initiales généralisées fournissent une bonne structure de données grâce aux algorithmes de cette partie.

## 16.3 Solutions polynomiales

Le principe est simple : on commence par trouver une borne sur le degré des solutions, puis la solution se ramène par coefficients indéterminés à un système linéaire. À nouveau, ce système a une structure induite par la récurrence linéaire, qu'il s'agit d'exploiter.

### Borne sur le degré

**Exemple 16.5** Si la récurrence

$$(n-3)u_{n+2} - (2n-3)u_{n+1} + nu_n = 0$$

a une solution qui se comporte pour $n$ grand comme $u_n \sim n^k$, alors le membre gauche est un polynôme de degré $k+1$. L'extraction des coefficients de degré $k+1, k, k-1, \dots$ fournit des équations qui doivent être satisfaites par une solution polynomiale. Il s'agit de

$$\text{degré } k+1: 1-2+1=0,$$

$$\text{degré } k: 2k-3-2k+3=0,$$

$$\text{degré } k-1: 2k(k-1)-6k-k(k-1)+3k=k(k-4)=0,$$

donc une solution polynomiale ne peut avoir que degré 0 ou 4.

Un algorithme pour ce calcul s'obtient un peu plus facilement en récrivant les décalages $u(n+k)$ de la suite initiale en termes de différences finies. Ainsi, on note $u(n)$ pour
$u_n$, $(\Delta \cdot u)(n) = u(n+1) - u(n)$ et, par récurrence, $(\Delta^{k+1} \cdot u)(n) = (\Delta^k \cdot u)(n+1) - (\Delta^k \cdot u)(n)$. La récurrence à annuler prend la forme

$$L \cdot u = \sum_{k=0}^{m} b_k(n) \Delta^k u = 0$$

pour de nouveaux polynômes $b_k$. L'opérateur $\Delta$ fait décroître de 1 exactement le degré des polynômes. Ainsi,

$$\deg(\Delta^k P) \leq \max(\deg P - k, 0),$$

avec égalité tant que $\Delta^k P$ n'est pas nul, et donc $\deg(L \cdot P) \leq \deg P + \max_k\{\deg(b_k) - k\}$. Soient alors les quantités

$$b := \max_k\{\deg(b_k) - k\}, \quad E := \{k \mid \deg(b_k) - k = b\},$$

qui vont servir à fournir une borne sur le degré des solutions.

**Exemple 16.6** La récurrence de l'exemple précédent se récrit

$$(n-3)(u_{n+2} - 2u_{n+1} + u_n) + (2(n-3) - (2n-3))u_{n+1} + (3 - n - 3 + n)u_n = 0,$$

qui se simplifie en

$$(n-3)\Delta^2 - 3\Delta)(u_n) = 0$$

l'entier $b$ vaut 0 et l'ensemble $E$ est $\{1, 2\}$.

Soit $D$ le degré d'une solution. La discussion distingue deux cas :

- soit $D + b < 0$, et alors $-(b+1)$ est une borne sur $D$;

- sinon le coefficient de degré $D + b$ dans $L(n^D + \dots)$ vaut

$$\sum_{k \in E} \text{ct}(b_k)D(D-1)\cdots(D-k+1),$$

où $\text{ct}$ désigne le coefficient de tête. Cette expression, vue comme un polynôme en $D$, s'appelle le *polynôme indiciel* de la récurrence, lequel est non nul.

Cette discussion mène au résultat suivant.

**Proposition 16.4** Une borne sur le degré des solutions polynomiales de l'opérateur $L = \sum b_k(n)\Delta^k$ est donnée par le maximum de $-(b+1)$ et de la plus grande racine entière positive du polynôme indiciel de $L$.

**Exemple 16.7** Dans l'exemple précédent, le polynôme indiciel vaut

$$D(D-1) - 3D = D(D-4),$$

ce qui redonne bien le résultat espéré.
## Algorithme naïf

Un algorithme simple consiste alors à calculer cette borne et à rechercher ensuite les solutions par un calcul d'algèbre linéaire.

**Exercice 16.1** Trouver les solutions polynomiales de la récurrence

$$3u(n + 2) - nu(n + 1) + (n - 1)u(n) = 0.$$

## Base binomiale

La borne sur le degré peut être exponentiellement grande en la taille binaire de la récurrence.

**Exemple 16.8** La récurrence

$$nu(n + 1) - (n + 100)u(n) = 0$$

a pour solution le polynôme $u(n) = n(n + 1) \cdots (n + 99)$ de degré 100.

Comme dans les sections précédentes, il est alors souhaitable de disposer d'algorithmes permettant de détecter l'existence de solutions polynomiales en complexité maîtrisée par rapport à N. De plus, lorsqu'existent des solutions polynomiales de degré N, les développer utilise O(N) coefficients, alors qu'elles peuvent être représentées de manière compacte par la récurrence et les conditions initiales généralisées correspondantes.

Une difficulté provient de ce que les solutions polynomiales de récurrences linéaires ont des coefficients dans la base $1, n, n^2, \dots$ qui ne vérifient pas une récurrence linéaire d'ordre borné a priori. En revanche, cette propriété apparaît si l'on développe les polynômes dans la base binomiale.

**Proposition 16.5** Si $(u_n)$ est une solution polynomiale de la récurrence (16.1) avec des coefficients de degré au plus $d$, alors ses coefficients $(c_k)$ dans la base binomiale $\{(\binom{n}{k}), k \in \mathbb{N}\}$ sont une suite solution à support fini d'une récurrence linéaire

$$a_m(k)c_{k+m} + \cdots + a_0(k)c_k = 0 \qquad (16.4)$$

qui ne présente pas de nouvelles singularités, dont l'ordre est au plus $d+r$, les $a_i$ ayant degré au plus $d$. Si les coefficients des $p_i$ sont entiers de taille bornée par $\ell$, ceux des $a_i$ sont du même ordre.

*Démonstration.* Le cœur de la preuve se résume au calcul suivant.
**Lemme 16.6** Si $u_n = \sum_{k=0}^{d} c_k \binom{n}{k}$ alors

$$nu_n = \sum_{k=0}^{d+1} k(c_k + c_{k-1}) \binom{n}{k}, \quad u_{n+1} = \sum_{k=0}^{d} (c_k + c_{k+1}) \binom{n}{k}.$$

*Démonstration.* La relation

$$\binom{n}{k+1} = \frac{n-k}{k+1} \binom{n}{k}$$

se récrit en

$$n \binom{n}{k} = (k+1) \binom{n}{k+1} + k \binom{n}{k}.$$

Le premier résultat s'obtient en injectant cette récriture dans la somme $nu_n$ et en extrayant le coefficient de $\binom{n}{k}$.

La seconde propriété découle directement du triangle de Pascal, qui s'obtient lui-même en extrayant le coefficient de $X^k$ dans l'identité $(1+X)^{n+1} = (1+X)^n(1+X)$, ce qui donne

$$\binom{n+1}{k} = \binom{n}{k} + \binom{n}{k-1}.$$

En utilisant répétitivement ces deux propriétés, tout terme $n^i u_{n+j}$ apparaissant dans une récurrence se traduit en une portion de récurrence pour les coefficients, les multiplications par $n$ pouvant augmenter l'ordre de 1, mais conservant le degré en $k$, tandis que les décalages en $n$ ne changent pas le degré, et mènent à des décalages en $k$, d'où le résultat de la proposition en ce qui concerne les degrés et ordre. Enfin, on voit que le coefficient de tête de la récurrence est conservé dans la récurrence transformée, qui a donc les mêmes singularités. ■

**Exercice 16.2** Montrer que l'existence d'une telle récurrence persiste pour toute base binomiale $\{(\binom{n+a}{k}), k \in \mathbb{N}\}$ avec $a \in \mathbb{K}$ quelconque. ■

**Exemple 16.9** Si $(c_k)$ est la suite de coefficients de $\binom{n}{k}$ dans une solution polynomiale de la récurrence de l'Exemple 16.8, alors la traduction de cette récurrence montre que cette suite vérifie

$$k(c_k + c_{k+1} + c_{k-1} + c_k) - k(c_k + c_{k-1}) - 100c_k = k c_{k+1} + (k-100)c_k = 0,$$

ce qui donne une récurrence d'ordre 1 pour ces coefficients, et le degré 100 de la solution.
**Exemple 16.10** La récurrence de l'Exemple 16.5 se traduit en

$$k(c_k + c_{k+1} + c_{k+1} + c_{k+2} + c_{k-1} + c_k + c_k + c_{k+1}) - 3(c_k + c_{k+1} + c_{k+1} + c_{k+2}) - 2k(c_k + c_{k+1} + c_{k-1} + c_k) + 3(c_k + c_{k+1}) + k(c_k + c_{k-1}) = (k-3)c_{k+2} + (k-3)c_{k+1} = 0,$$

équivalent à $c_0$ libre et $(k-4)c_{k+1} + (k-4)c_k = 0$ pour $k \ge 1$, où il apparaît clairement que $k = 4$ mène à une solution à support fini, de degré 4.

### Conditions initiales

Comme pour les récurrences satisfaites par les coefficients de séries solutions d'équations différentielles linéaires (§14.1 page 253), l'ordre de la récurrence obtenue n'est pas nécessairement le même que celui de la récurrence initiale, et la traduction fournit également des contraintes linéaires sur les conditions initiales, qui précisent le sous-espace des solutions de la nouvelle récurrence correspondant aux solutions de celle de départ.

**Exemple 16.11** La récurrence

$$v_{n+2} - v_{n+1} - (2n + 1)v_n = 0$$

se traduit en

$$(c_k + 2c_{k+1} + c_{k+2}) - (c_k + c_{k+1}) - 2k(c_k + c_{k-1}) - c_k = c_{k+2} + c_{k+1} - (2k + 1)c_k - 2kc_{k-1} = 0$$

pour $k \ge 1$, et une contrainte obtenue en faisant $k = 0$ et $c_{-1} = 0$ dans cette équation, ce qui donne $c_2 + c_1 - c_0 = 0$. Ainsi la récurrence d'ordre 2 mène à une récurrence d'ordre 3, plus une contrainte linéaire sur les conditions initiales.

### Conditions initiales généralisées

Si $m$ est indice d'une condition initiale généralisée de la récurrence dont la suite $(u_n)$ est solution, alors le lien entre les deux suites $(u_n)$ et $(c_k)$ donne

$$u_m = \sum_{k=0}^{m} c_k \binom{m}{k}$$

puisque les binomiaux $\binom{m}{k}$ sont nuls pour $k > m$. La condition initiale généralisée est donc donnée par

$$c_m = u_m - \sum_{k=0}^{m-1} c_k \binom{m}{k}.$$

Pour calculer cette valeur efficacement lorsque $m$ est grand, il suffit de calculer d'abord une récurrence pour la suite $(c_k(\binom{m}{k}))_k$, par clôture (Thm. 14.6), puis pour la somme $(\sum_{k=0}^i c_k(\binom{m}{k}))_i$ en augmentant l'ordre de 1. La somme est alors calculée par l'Algorithmme 15.4. Ainsi les conditions initiales généralisées sont calculées efficacement de proche en proche.
**Entrée** Une récurrence linéaire de la forme

$$p_r(n)u_{n+r} + \cdots + p_0(n)u_n = q(n).$$

**Sortie** Une base des solutions polynomiales, sous la forme d'une récurrence sur les coefficients des solutions dans la base binomiale, et de conditions initiales généralisées.

1. Calculer les conditions initiales généralisées d'une base des solutions par l'Algorithme 16.2.

2. Si la récurrence est inhomogène ($q \neq 0$), l'homogénéiser (§15.4 page 283).

3. Calculer la récurrence

$$a_m(k)c_{k+m} + \cdots + a_0(k)c_k = 0$$

donnée par la Proposition 16.5, ainsi que les équations linéaires sur les conditions initiales généralisées correspondant aux suites solutions de l'équation de départ.

4. Utiliser l'Algorithme 16.4 sur cette récurrence, en y adjointant ces équations linéaires avant la résolution.

5. Renvoyer la récurrence et les conditions initiales généralisées de la base trouvée.

**Algorithme 16.5** – Solutions polynomiales dans la base binomiale.

## Solutions polynomiales

Grâce à cette traduction, la recherche de solutions polynomiales se ramène donc à la recherche de solutions à support fini. L'ensemble du calcul est résumé dans l'Algorithmme 16.5. Les complexités arithmétiques et binaires découlent des considérations précédentes et sont résumées dans le théorème suivant.

**Théorème 16.7** Les racines entières positives des coefficients $a_0$ et $a_m$ de la récurrence (16.4) étant données, si N est la plus grande racine entière positive de $a_0$ et $M \le d$ une borne sur le nombre de racines entières positives de $a_m$ inférieures à N, alors on peut déterminer les solutions polynomiales de la récurrence (16.1) en

$$O((r+d)^{\theta} (M(\sqrt{d(M+1)N}))\log(dN) + M(1+M/m))$$

opérations arithmétiques, sous la forme des conditions initiales généralisées d'une base de solutions exprimées dans la base binomiale.

Si de plus les coefficients $p_i$ sont des entiers de taille bornée par $\ell$, la complexité binaire du calcul de ces racines entières et du reste de l'algorithme est

$$O((r+d+M)^2 M_Z(T) \log T + (r+d)M d^2 M_Z(\ell + \log N)),$$

où $T = (r+d)(dN \log N + \ell N + N \log(r+d))$ est une borne sur les conditions initiales généralisées de la solution.
## 16.4 Solutions rationnelles

La recherche de solutions rationnelles se déroule en deux étapes : trouver d'abord un multiple du dénominateur, puis effectuer un changement de suite inconnue et chercher une solution polynomiale à la nouvelle équation ainsi obtenue pour le numérateur.

### L'algorithme d'Abramov

Pour déterminer un multiple du dénominateur, il faut d'abord repérer les positions des pôles éventuels d'une solution rationnelle. L'intuition de départ s'obtient en considérant le cas de fractions dans le plan complexe. Si $u(n) = P(n)/Q(n)$ est solution de la récurrence

$$p_r(n) \frac{P(n+r)}{Q(n+r)} + \cdots + p_0(n) \frac{P(n)}{Q(n)} = q(n) \quad (16.5)$$

alors chaque racine $a$ de $Q(n)$ donne lieu à une racine $a-1$ de $Q(n+1)$, une racine $a-2$ de $Q(n+2)$ et ainsi de suite jusqu'à une racine $a-r$ de $Q(n+r)$. Pourtant, le membre gauche de la récurrence doit rester un polynôme. Cela signifie que ces racines doivent correspondre à des racines des coefficients correspondants. Si aucune paire de racines de $Q$ ne différait d'un entier, alors le dénominateur serait nécessairement un multiple du pgcd des polynômes $p_r(n-r), \dots, p_1(n-1), p_0(n)$. La situation est plus complexe si deux racines de $Q$ diffèrent d'un entier. Cependant, si $\alpha$ et $\beta$ sont deux racines de $Q$ telles que $\alpha - \beta = k \in \mathbb{N}$ est maximal, alors nécessairement $p_0(\alpha) = 0$ et $p_r(\beta - r) = 0$. La plus grande différence entière $H$ entre les racines de $Q$ peut donc être bornée étant donnés $p_0$ et $p_r$.

Cette observation peut être raffinée et abstraite du contexte du plan complexe. Soit

$$M(n) = \text{ppcm}(Q(n+1), \dots, Q(n+r)).$$

La récurrence multipliée par $M(n)$ s'écrit

$$p_r(n)P(n+r) \frac{M(n)}{Q(n+r)} + \cdots + p_0(n)P(n) \frac{M(n)}{Q(n)} = 0.$$

Dans cette récurrence, tous les coefficients $M(n)/Q(n+k)$, $k=1, \dots, r$ sont des polynômes. Par conséquent $Q$ divise $p_0M$, ce qui se récrit

$$Q(n) | p_0(n) \text{ppcm}(Q(n+1), \dots, Q(n+r)).$$

En décalant $n$ et en reportant, il vient

$$Q(n) \begin{aligned} &| p_0(n) \text{ppcm}(p_0(n+1) \text{ppcm}(Q(n+2), \dots, Q(n+r+1)), Q(n+2), \dots, Q(n+r)) \\ &| p_0(n)p_0(n+1) \text{ppcm}(Q(n+2), \dots, Q(n+r+1)) \\ &| p_0(n) \cdots p_0(n+j) \text{ppcm}(Q(n+j+1), \dots, Q(n+j+r)). \end{aligned}$$

Cependant, pour $j > H$, $\text{pgcd}(Q(n), Q(n+j)) = 1$ et donc nous avons obtenu

$$Q(n) | p_0(n) \cdots p_0(n+H).$$
**Entrée** Une récurrence linéaire de la forme

$$p_r(n)u_{n+r} + \cdots + p_0(n)u_n = q(n).$$

**Sortie** Un multiple du dénominateur des solutions rationnelles.

1. Calculer le polynôme

$$R(h) = \mathrm{Res}_n(p_0(n+h), p_r(n-r)).$$

2. Si R n'a pas de racines dans $\mathbb{N}$, alors renvoyer le polynôme 1 ;
sinon, soit $h_1 > h_2 > \cdots > h_m \ge 0$ ses racines entières positives.
Initialiser Q à 1, A à $p_0(n)$, B à $p_r(n-r)$.

3. Pour $i = 1, \dots, m$ faire

$$g(n) := \mathrm{gcd}(A(n+h_i), B(n));$$

$$Q(n) := g(n)g(n-1)\cdots g(n-h_i)Q(n);$$

$$A(n) := A(n)/g(n-h_i);$$

$$B(n) := B(n)/g(n).$$

4. Renvoyer Q.

**Algorithme 16.6** – Algorithme d'Abramov pour les multiples du dénominateur.

De la même manière, multiplier le polynôme ppcm($Q(n)$, ..., $Q(n+r-1)$) par la récurrence permet de voir que

$$Q(n) \mid p_r(n-r) \cdots p_r(n-r-H)$$

et finalement

$$Q(n) \mid \mathrm{gcd}(p_0(n) \cdots p_0(n+H), p_r(n-r) \cdots p_r(n-r-H))$$

fournit un multiple de Q. Nous avons ainsi prouvé la première partie du résultat suivant.

**Proposition 16.8** L'algorithme 16.6 calcule un multiple du dénominateur des éventuelles solutions rationnelles de la récurrence (16.1). Si les racines entières positives du polynôme $R(h)$ sont données, sa complexité arithmétique est bornée par $O(M(d^2) + dM(d)\log d)$ où $d$ est une borne sur le degré des coefficients de la récurrence, à condition de ne pas développer les produits intervenant dans l'expression du polynôme renvoyé.

*Démonstration.* Le polynôme dans l'étape (1) est un résultant particulier qui peut se calculer efficacement comme une somme composée. Ce calcul est présenté dans le Chapitre 3 comme application du calcul rapide de l'exponentielle des séries. Sa complexité arithmétique est en $O(M(d^2))$.

Si on conserve le polynôme Q sous la forme

$$Q(n) = \prod_{i} \prod_{j=0}^{h_i} g_i(n-j),$$
où $g_i$ désigne le pgcd à l'étape $i$, le coût de la boucle est dominé par celui du pgcd. Tout d'abord, il faut décaler le polynôme $A$ dont le degré reste borné par $d$, ce qui se fait à nouveau par une somme composée en $O(M(d))$ opérations arithmétiques. Puis chaque pgcd a une complexité arithmétique en $O(M(d) \log d)$ et le degré de $B$ décroissant à chaque étape, le nombre d'étapes est au plus $d$, d'où le résultat. ■

L'estimation de la complexité binaire n'est pas compliquée, mais plus technique. Il faut borner les tailles des coefficients des différents polynômes et recourir au théorème des restes chinois pour les reconstruire efficacement.

**Exercice 16.3** Trouver les dénominateurs possibles des solutions rationnelles des récurrences suivantes :

$$ 
\begin{aligned}
(n+1)u_{n+1} - nu_n &= 0, \\
(n+1)(n+3)u_{n+2} - 2(n+2)nu_{n+1} + (n-1)(n+1) &= 0, \\
(n+3)(n+2)(n^2 + 6n + 4)u(n+2) - (n+1)(3n^3 + 27n^2 + 64n + 48)u(n+1) \\
&\quad + 2n^2(n^2 + 8n + 11)u(n) = 0.
\end{aligned}
$$

**Exemple 16.12** Pour trouver les dénominateurs possibles pour les solutions rationnelles de la récurrence

$$ (2n-1)(4n+4K+5)(4n+4K+3)u_{n+1} + 2(4n+5)(4n+3)(n+1)u_n = q(n), $$

où $K$ est un entier strictement positif et $q$ un polynôme, on commence par chercher les différences entières positives entre les racines de $(4n+5)(4n+3)(n+1)$ et $(2n-3)(4n+4K+1)(4n+4K-1)$, et il est facile de voir qu'il n'y a que $K-1$. Le calcul par résultant peut se faire aussi, il s'agit de considérer

$$ \mathrm{Res}_n(2(4n+4h+5)(4n+4h+3)(n+h+1), (2n-3)(4n+4K+1)(4n+4K-1)), $$

pour lequel on peut exploiter la multiplicativité du résultant (Cor. 6.6), ramenant ce calcul à des résultats de polynômes unitaires de degré 1, exhibant directement les différences de racines.

Il n'y a ensuite qu'un calcul de pgcd, qui extrait les facteurs menant à la racine commune $K-1 : g(n) = (4n + 4K + 1)(4n + 4K - 1)$. Un multiple des dénominateurs possibles est donc

$$ \prod_{j=0}^{K-1} (4n + 4K - 4j + 1)(4n + 4K - 4j - 1). $$

## Numérateur

Une fois calculé ce multiple $Q(n)$ du dénominateur, le changement d'inconnue $u(n) = P(n)/Q(n)$ dans l'équation (16.5) peut être réalisé efficacement, même si les $h_i$
sont très grands. C'est clair dans le cas où l'équation est homogène ($q(n) = 0$) puisqu'alors il suffit de multiplier le membre gauche par $Q(n)$ et d'observer que

$$ \frac{Q(n)}{Q(n+k)} = \prod_{i} \prod_{j=1}^{h_i} \frac{g_i(n-j)}{g_i(n+k-j)} = \prod_{i} \frac{g_i(n+k-h_i-1)g_i(n+k-h_i-2)\cdots g_i(n-h_i)}{g_i(n+k-1)g_i(n+k-2)\cdots g_i(n)} $$

est une fraction dont le numérateur et le dénominateur ont degré borné par

$$ k \sum \deg(g_i) \leq k \max(\deg p_0, \deg p_r), $$

c'est-à-dire qui ne souffre pas de la taille potentiellement grande des $h_i$.

**Exercice 16.4** Finir de trouver les solutions rationnelles pour les récurrences de l'Exercice 16.3. ■

### Le cas inhomogène

Dans le cas inhomogène, ce changement d'inconnue mène à un membre droit de grand degré dépendant des $h_i$, mais celui-ci peut-être réduit : si le changement d'inconnue donne

$$ A(n, P(n)) = Q(n) \sum p_i(n) \frac{P(n+i)}{Q(n+i)} = Q(n)q(n), $$

l'homogénéisation vue en §15.4 fournit une équation homogène linéaire

$$ q(n)A(n+1, P(n+1)) - q(n+1) \frac{Q(n+1)}{Q(n)} A(n, P(n)) = 0 \quad (16.6) $$

dont les degrés des coefficients restent modérés. L'homogénéisation fournit aussi des équations sur les conditions initiales généralisées supplémentaires. Le problème est ainsi ramené à la recherche d'une base de solutions polynomiales, par l'Algorithmie 16.5.

### Le cas inhomogène paramétré

La récurrence (16.6) obtenue par homogénéisation d'une récurrence dont la partie inhomogène est $\lambda_1 q_1(n) + \cdots + \lambda_k q_k(n)$ produit a priori des divisions par des polynômes en les $\lambda_i$ puisque $q(n)$ divise son coefficient de tête, mais les solutions restent linéaires en les $\lambda_i$ et les conditions initiales en prenant en compte les contraintes linéaires sur les conditions initiales obtenues lors de l'homogénéisation.

**Exemple 16.13** Pour trouver les valeurs de $\lambda, \mu$ permettant l'existence de solutions rationnelles à la récurrence

$$ (2n-1)(2n-3)(4n+4K+5)(4n+4K+3)u_{n+1} + 2(4n+5)(4n+3)(2n-3)(n+1)u_n = 2\lambda(4n+5)(4n+3)(2n-3)(n+1) + 6\mu(4n+4K+3)(4n+5)(n+1), $$

où K est un entier strictement positif, la première étape est de chercher un multiple des dénominateurs possibles, qui est fourni par la partie homogène. Une variante
simple de l'Exemple 16.12 donne le polynôme

$$
Q(n) := (2n-3) \prod_{j=0}^{K-1} ((4n + 4K - 4j)^2 - 1).
$$

Le facteur supplémentaire (2n − 3) peut soit être déduit du fait qu’il divise le
membre gauche de la récurrence, soit obtenu automatiquement par l’Algorithme
d’Abramov, qui fait apparaître 0 comme différence entière entre les coefficients.

Le changement d’inconnue $u_n = v_n/Q(n)$ suivi par une multiplication par Q donne la récurrence

$$
(2n-3)(4n+4K+5)(4n+4K+3) \frac{\prod_{j=0}^{K-1} ((4n+4K-4j)^2 - 1)}{\prod_{j=0}^{K-1} ((4n+4K+4-4j)^2 - 1)} v_{n+1} \\
+ 2(4n+5)(4n+3)(n+1)v_n = \\
(4n+5)(n+1)(2\lambda(4n+3)(2n-3) + 6\mu(4n+4K+3)) \prod_{j=0}^{K-1} ((4n+4K-4j)^2 - 1),
$$

où le facteur de $v_{n+1}$ se simplifie en $(2n-3)(4n+3)(4n+5)$ et il ne reste donc des
polynômes dont le degré dépend de $K$ (donc potentiellement grands) que dans le
membre droit. L’homogénéisation fournit alors une récurrence à coefficients de
faible degré :

$$
\begin{align*}
& (4n+5)(4n+7)(-1+2n)(n+1)((4n+3)(2n-3)\lambda + (12n+12K+9)\mu)v_{n+2} \\
& -(n+2)((4n+7)(2n-3)(16(2n-1)K^2 + 32(2n^2 + n-1)K - 48n^2 - 96n - 45)\lambda \\
& \quad +3(4n+4K+3)(16(2n-3)K^2 + 8(2n^2 - 13)K - 80n^2 - 240n - 175)\mu)v_{n+1} \\
& -2(n+2)(n+1)(4n+4K+5)(4n+4K+3)((4n+7)(2n-1)\lambda + 3(4n+4K+7)\mu)v_n = 0,
\end{align*}
$$

ainsi que la contrainte

$$
-9v_{1} + 6v_{0} = (6\lambda - 2(4K + 3)\mu)Q(0)
$$

sur les conditions initiales. Cette dernière est très importante, puisqu'elle permet
de rester dans un sous-espace où les suites ont pour éléments des combinaisons
linéaires de λ, μ, 1 et non des fractions rationnelles en λ, μ. Par ailleurs, à ce stade, on
observe que quitte à multiplier Q par un facteur constant, on peut prendre Q(0) =
1.

Il faut ensuite chercher les solutions polynomiales de cette récurrence sur (v_n).
La traduction en récurrence sur les coefficients dans la base binomiale mène
à une récurrence d'ordre 7 un peu grosse pour être reproduite ici, ainsi que
des contraintes linéaires sur les conditions initiales, auxquelles il faut ajouter la
contrainte

$$
-9(c_0 + c_1) + 6c_0 = (6\lambda - 2(4K + 3)\mu)Q(0)
$$

déduite de la précédente. Le coefficient de $c_k$ dans la récurrence est le polynôme
indiciel

$$
512\lambda(k+5)(k+4)(k+3)(k+2)(k+1)(k-2K-1)
$$

dont la seule racine entière positive est 2K + 1 qui est donc le seul degré possible
pour une solution polynomiale.

La suite du calcul consiste donc à obtenir les équations sur $c_0, \lambda, \mu$ permettant l'existence d'une solution à support fini. Ces équations sont produites en utilisant la récurrence d'ordre 7 et selon le cas un calcul naïf ou de « pas de bébés-pas de géants » ou de scindage binaire pour obtenir (avec des choix $(c_0, \lambda, \mu) \in \{(0,1,0), (1,1,0), (0,0,1)\}$), les valeurs de $c_{2K+2}, c_{2K+3}, \dots, c_{2K+8}$ comme combinaisons linéaires de $c_0, \lambda, \mu$. L'existence de solutions à ce système linéaire est équivalente à celle de solutions rationnelles à l'équation de départ.

Cette partie du calcul ne peut pas être effectuée avec K symbolique, mais par
exemple avec K = 10, on obtient l'existence d'une solution avec

$$
c_0 = 0, \quad \lambda = \frac{4790112341450}{13996552017} \mu.
$$

À ce stade, tout dépend de ce que l'on souhaite calculer à propos de cette fraction rationnelle. Si l'objectif est d'obtenir les coefficients du numérateur, alors on peut utiliser la récurrence sur les $(c_k)$ avec ces contraintes. On constate alors que le numérateur est un polynôme de degré 21, comme prédit par la borne, qui a un facteur de degré 1, le polynôme $n$, correspondant à $c_0 = 0$, et dont l'autre facteur est irréductible.

Exercice 16.5 Résoudre en (u, λ, μ) la récurrence

$$
3u(n+2) - nu(n+1) + (n-1)u(n) = \lambda n^3 + \mu n^2.
$$

Pour conclure, l'ensemble du calcul de solutions rationnelles, résumé par l'Al-
gorithme 16.7, peut ainsi être réalisé efficacement, même si certains polynômes ont
des degrés exponentiels en la taille de l'équation, à condition de ne pas développer
ces polynômes, sous la forme d'un multiple du dénominateur exprimé comme un
produit, et d'une récurrence avec ses conditions initiales généralisées pour une base
des numérateurs exprimés dans la base binomiale.

L'ordre 1

Lorsque la récurrence est d'ordre 1, il est en outre possible de prédire un facteur
du numérateur, ce qui gagne en efficacité pratique. En effet, dans la récurrence

$$
a(n)u(n+1) + b(n)u(n) = c(n),
$$
**Entrée** Une récurrence linéaire de la forme

$p_r(n)u_{n+r} + \cdots + p_0(n)u_n = q(n)$, où $q(n)$ peut être nul, ou un polynôme, ou une combinaison linéaire de polynômes à coefficients $\lambda_i$ inconnus.

**Sortie** Une base des solutions rationnelles, et éventuellement les valeurs des $\lambda_i$ en permettant l'existence.

1. Calculer un multiple Q du dénominateur par l'Algorithme 16.6.

2. Poser $u_n = v_n/Q(n)$ et former la récurrence qui en découle sur $v_n$ ; éventuellement l'homogénéiser.

3. Calculer les solutions polynomiales de cette récurrence en procédant comme dans l'Algorithme 16.5 à partir de sa troisième étape.

4. Renvoyer Q, ces solutions et la récurrence correspondante.

**Algorithme 16.7** – Solutions rationnelles de récurrences.

si $b$ et $c$ ont une racine commune qui n'est ni un pôle de $u$ ni une racine de $a$, alors celle-ci est nécessairement racine de $u(n+1)$. De même, si $a$ et $c$ ont une racine commune qui n'est ni un pôle de $u(n+1)$ ni une racine de $b$, alors elle est racine de $u$.

Ces calculs préalables permettent de réduire le degré des coefficients de l'équation dont on recherche ensuite les solutions polynomiales. Ils sont utilisés dans l'algorithme de Gosper présenté au Chapitre 29.

## Exercices

Ce chapitre montre qu'il est possible de déterminer les solutions à support fini de récurrences du type (16.1), puis par extension les solutions polynomiales de récurrences sous une forme compacte donnée par : la récurrence elle-même à coefficients dans $\mathbb{Z}[X]$ de degré borné par $d$ ; une borne $N$ sur le degré ; des conditions initiales généralisées (voir page 280) rationnelles de taille $O(dN \log N)$. Ceci n'est utile que si cette forme compacte se prête bien au calcul. L'exercice suivant et l'Exercice 17.4 du chapitre suivant explorent cette question.

**Exercice 16.6** (Forme compacte en base binomiale.)

1. Montrer que si $(u_k)$ est solution d'une récurrence linéaire de la forme (16.1) et si $\alpha$ n'est pas un entier négatif, la suite $(\sum_{i=0}^k u_i(\alpha))$ est solution d'une récurrence linéaire à coefficients polynomiaux que l'on peut calculer efficacement, ainsi que les $r$ premières conditions initiales, dont le terme de tête n'a pas de nouvelles racines entières positives.

2. Montrer que si $P$ est donné sous forme compacte dans la base binomiale, on peut calculer $P(a)$ efficacement.

3. Même question pour le polynôme $\Delta^l P(a)$ (ici $\Delta$ est l'opérateur aux différences est le polynôme $\Delta : P(X) \mapsto P(X+1) - P(X)$, et $\Delta^l$ représente l'itérée $\ell$-ième de cet opérateur).
## Notes

L'existence d'une récurrence pour les coefficients dans la base binomiale pour les solutions de récurrences linéaires est connue au moins depuis le 19e siècle [Boo72, Chap. XIII, art. 5]. Elle a été utilisée pour donner un algorithme calculant les solutions polynomiales par Abramov, Bronstein, Petkovšek [ABP95a]. Ensuite, ces questions ont été reprises du point de vue de l'efficacité en tenant compte de la taille potentiellement exponentielle des degrés en jeu, par l'introduction des algorithmes rapides pour la factorielle de matrice dans ce contexte, ainsi que les adaptations aux cas inhomogène et inhomogène paramétré. Ces idées, ainsi que les Exercices 17.4 et 16.6 proviennent d'articles de Bostan, Chyzak, Cluzeau, et Salvy [BCS05; Bos+06a].

Notre présentation de la preuve de l'algorithme d'Abramov est due à Chen, Paule et Saad [CPS08]. Abramov a proposé plusieurs raffinements de son algorithme pour obtenir des multiples du dénominateur de degré moindre [Abr95] et aussi dans certains cas en tenant compte de tous les coefficients $p_i$ [Abr89]. Une manière plus directe d'aboutir à ces raffinements a été décrite par van Hoeij [Hoe98].

## Bibliographie

ABP95a ABRAMOV, Sergei A., Manuel BRONSTEIN et Marko PETKOVŠEK (1995). « On Polynomial Solutions of Linear Operator Equations ». In: Symbolic and Algebraic Computation. Éd. par A. H. M. LEVELT. Proceedings of ISSAC'95, July 1995, Montreal, Canada. New York: ACM Press, p. 290-296.

Abr89 ABRAMOV, S. A. (1989). « Rational solutions of linear differential and difference equations with polynomial coefficients ». In: USSR Computational Mathematics and Mathematical Physics, vol. 29, n°11, p. 1611-1620.

Abr95 — (1995). « Rational solutions of linear difference and q-difference equations with polynomial coefficients ». In: ISSAC'95: International Symposium on Symbolic and Algebraic Computation. Éd. par A. H. M. LEVELT. New York: ACM Press, p. 285-289.

BCS05 BOSTAN, Alin, Thomas CLUZEAU et Bruno SALVY (2005). « Fast algorithms for polynomial solutions of linear differential equations ». In: ISSAC'05: International Symposium on Symbolic and Algebraic Computation. ACM Press, p. 45-52.

Boo72 BOOLE, George (1972). A treatise on the calculus of finite differences. 2nd. Macmillan.

Bos+06a BOSTAN, A., F. CHYZAK, T. CLUZEAU et B. SALVY (2006). « Low complexity algorithms for linear recurrences ». In: ISSAC'06: International Symposium on Symbolic and Algebraic Computation. ACM Press, p. 31-38.

CPS08 CHEN, William Y. C., Peter PAULE et Husam L. SAAD (2008). « Converging to Gosper's algorithm ». In: Advances in Applied Mathematics, vol. 41, n°3, p. 351-364.

Hoe98 HOEIJ, Mark van (1998). « Rational solutions of linear difference equations ». In: ISSAC'98: International Symposium on Symbolic and Algebraic Computation. New York: ACM Press, p. 120-123.
# 17. Résolution d'équations différentielles linéaires

## Résumé

Comme pour les récurrences, les solutions polynomiales ou rationnelles d'équations différentielles linéaires peuvent avoir une taille exponentiellement grande par rapport à la taille de l'équation. Néanmoins, elles s'obtiennent efficacement en exploitant les algorithmes rapides sur les solutions de récurrences linéaires. Par ailleurs, la méthode de scindage binaire permet l'évaluation numérique garantie et rapide de solutions d'équations différentielles linéaires en complexité quasi-optimale vis-à-vis de la précision.

Comme au Chapitre 16, trois types d'équations sont considérés ici :

- homogènes, c'est-à-dire de la forme

$$c_m(x)y^{(m)}(x) + \dots + c_0(x)y(x) = 0, \quad (17.1)$$

où les coefficients $c_i$ sont des polynômes donnés à coefficients dans un corps $K$;

- inhomogènes, c'est-à-dire du même type, mais avec un membre droit qui n'est plus simplement 0, mais un polynôme $b(x)$;

- inhomogènes paramétrées, où le membre droit est maintenant de la forme

$$\lambda_1 b_1 + \dots + \lambda_k b_k,$$

où les $b_i$ sont des polynômes donnés, et les coefficients $\lambda_i$ sont inconnus.

Les solutions étudiées sont polynomiales, rationnelles, ou en séries. Les questions à aborder concernant ces solutions portent sur leur existence, leur valeur (et la structure de données à utiliser pour représenter cette valeur), la complexité arithmétique et binaire des calculs correspondants. En outre, dans le cas inhomogène paramétré, il s'agira de déterminer l'existence et la valeur des coefficients $\lambda_i$ permettant l'existence d'une solution rationnelle de l'équation, cette opération intervenant dans les algorithmes de la Partie VI sur la sommation et l'intégration.
## 17.1 Séries solution

### Traduction en récurrence

Le Théorème 14.1 page 253 montre que la suite des coefficients d'une série formelle solution de l'équation différentielle linéaire homogène (17.1) vérifie une récurrence linéaire homogène de type (16.1). Si les coefficients de l'équation différentielle sont

$$c_i(x) = \sum_j c_{i,j} x^j,$$

le calcul explicite de la récurrence vérifiée par la suite $(y_n)$ des coefficients de la série $y(x)$ donne

$$\sum_{i,j} c_{i,j} (n-j+1) \cdots (n-j+i) y_{n+i-j} = 0, \quad (17.2)$$

pour tout $n \in \mathbb{Z}$, avec $y_\ell = 0$ pour $\ell < 0$. Quitte à décaler les indices, cette récurrence se récrit sous la forme

$$p_r(n)y_{n+r} + \cdots + p_0(n)y_n = 0, \quad (17.3)$$

à laquelle il faut ajouter les équations

$$p_r(-1)y_{r-1} + \cdots + p_0(-1)y_0 = 0, \quad p_r(-2)y_{r-2} + \cdots = 0, \dots, \quad p_r(-r)y_0 = 0, \quad (17.4)$$

qui tiennent compte des contraintes $y_\ell = 0$ pour $\ell < 0$. (Voir l'Exemple 14.4 page 254.)

**Définition 17.1** Le polynôme $p_r(n-r)$ est appelé polynôme indiciel en 0 de l'équation (17.1).

**Exercice 17.1** Le polynôme $p_0(n)$ est appelé polynôme indiciel à l'infini. En changeant la variable $x$ en $x + \alpha$ dans l'équation, le même calcul fournit deux polynômes. Le premier s'appelle polynôme indiciel en $\alpha$. Montrer que le polynôme indiciel à l'infini est inchangé. ■

### Séries solutions

**Proposition 17.1** Si $F \in \mathbb{K}[[x]]$ est solution de l'équation différentielle (17.1) alors sa valuation valF annule le polynôme indiciel en 0 de (17.1).

*Démonstration.* Il suffit d'évaluer la récurrence (17.3) en $n = \text{val}(F) - r$. ■

**Corollaire 17.2** Si le coefficient de tête $c_m$ de l'équation différentielle (17.1) est non nul en 0, alors cette équation admet une base de $m$ séries solutions ayant pour valuations $\{0, 1, \dots, m-1\}$.

*Démonstration.* L'équation (17.2) montre que lorsque $c_{m,0} \neq 0$, le polynôme indiciel en 0 est donné par

$$p_r(n-r) = c_{m,0}(n-m+1)\cdots(n),$$

et donc la récurrence est régulière au sens du Chapitre 15. ■
Pour l'application importante au calcul des N premiers termes d'une série solution, on obtient donc.

**Corollaire 17.3** Si le coefficient de tête $c_m$ de l'équation différentielle (17.1) est non nul en 0, alors étant données ses $m$ conditions initiales, on peut calculer les N premiers coefficients du développement en série d'une solution en

$$O\left((m+d)\left(M(d)\log d + \frac{M(m)}{m}N\right)\right)$$

opérations arithmétiques.

Si de plus les coefficients $c_i$ de l'équation différentielle sont dans $\mathbb{Z}[x]$ avec des coefficients de taille bornée par $\ell$, et les conditions initiales sont des rationnels de numérateurs et dénominateurs bornés par $\tilde{O}(N(m \log N + d + \ell))$, alors la taille totale de ces N premiers coefficients est bornée par $\tilde{O}(N^2(m \log N + d + \ell))$ et leur calcul peut s'effectuer en

$$\tilde{O}\left((m+d)N^2 M_{\mathbb{Z}}(m \log N + d + \ell)\right)$$

opérations binaires.

*Démonstration.* La traduction d'équation différentielle en récurrence coûte

$$O((m+d)M(d)\log d)$$

opérations arithmétiques d'après le Théorème 14.5. D'après l'équation (17.3), la récurrence obtenue a ordre au plus $d+m$, degré au plus $m$, et ses coefficients, des combinaisons linéaires des coefficients des $c_i$ avec des facteurs bornés par $(d+m)!$, ont donc taille en $\tilde{O}(d+m+\ell)$ lorsque les $c_i \in \mathbb{Z}[x]$ ont des coefficients de taille au plus $\ell$. La conclusion s'obtient par application de la Proposition 15.1. ■

## Séries généralisées solutions

Les calculs ci-dessus s'étendent facilement à une classe un peu plus grande de solutions.

**Définition 17.2** On appelle série généralisée une série de la forme $x^\alpha F$, où $F \in \mathbb{K}[[x]]$ et $\alpha \in \mathbb{K}$.

**Proposition 17.4** Si l'équation différentielle (17.1) admet une solution série généralisée alors $\alpha$ est racine du polynôme indiciel en 0 de (17.1). Si de plus, $\alpha \notin \{0, 1, \dots, m-1\}$, alors $c_m(0) = 0$.

*Démonstration.* L'équation

$$x^{-\alpha} (c_m(x)(x^{\alpha} F)^{(m)} + \dots + c_0(x) x^{\alpha} F) = 0$$
est de même type que (17.1). La récurrence (17.3) qui lui correspond s'obtient en remplaçant $n$ par $n - \alpha$ dans (17.3). Le résultat se déduit alors des deux résultats précédents. ■

### Translation

Les solutions séries généralisées de la forme

$$y(t) = (t-a)^{\alpha} \sum_{i=0}^{\infty} y_i (t-a)^i \quad \text{avec } y_0 \neq 0$$

se ramènent au cas précédent en effectuant le changement d'inconnue $y(a+t) = \tilde{y}(t)$ qui fournit l'équation

$$c_m(a+t)\tilde{y}(t)^{(m)} + \dots + c_0(a+t)\tilde{y}(t) = 0. \quad (17.5)$$

Les résultats obtenus jusqu'ici établissent donc la proposition suivante.

**Proposition 17.5** Si l'équation différentielle (17.1) admet une solution de la forme $(t-a)^{\alpha} \sum_{i=0}^{\infty} y_i (t-a)^i$ avec $y_0 \neq 0$ et $\alpha \notin \{0, \dots, m-1\}$ alors $c_m(a) = 0$ et $\alpha$ est racine du polynôme indiciel de (17.5).

**Définition 17.3** On appelle polynôme indiciel de l'équation différentielle (17.1) en $a$ celui de l'équation (17.5) en 0.

**Définition 17.4** On dit que $\alpha \in \overline{\mathbb{K}}$ est un point ordinaire de l'équation (17.1) si $c_m(\alpha) \neq 0$. On dit qu'il est singulier dans le cas contraire. On dit que l'infini est un point singulier de l'équation si 0 est singulier pour l'équation satisfaites par $y(1/x)$.

**Exemple 17.1** La fraction rationnelle $y = 1/(1-x)$ est solution de l'équation

$$(1-x)y'(x) - y(x) = 0.$$

Le nombre complexe 1 est singularité de la solution et donc nécessairement point singulier de l'équation, ce qui se traduit par l'annulation du coefficient de tête.

**Exemple 17.2** Le polynôme $x^{10}$ est solution de l'équation

$$10xy'(x) - y(x) = 0.$$

La solution n'a pas de point singulier à distance finie, mais le nombre complexe 0 est un point singulier de l'équation; le théorème de Cauchy ne s'y applique pas.
## Preuves de non-finitude différentielle

Une conséquence utile de la Proposition 17.5 est que les fonctions différentiellement finies ne peuvent avoir qu'un nombre fini de singularités (les racines du coefficient de tête). Il s'en déduit des résultats négatifs.

**Exemple 17.3** Les fonctions $1/\sin x$ et $\sqrt{\sin x}$ ne satisfont pas d'équation différentielle linéaire à coefficients polynomiaux.

**Exemple 17.4** La suite des nombres de Bernoulli, qui interviennent en particulier dans la formule d'Euler-Maclaurin, ne peut être solution d'une récurrence linéaire à coefficients polynomiaux, puisque la série génératrice exponentielle

$$ \sum_{n=0}^{\infty} B_n \frac{z^n}{n!} = \frac{z}{\exp(z) - 1} $$

a une infinité de pôles (aux points $2ik\pi$, pour $k \in \mathbb{Z} \setminus \{0\}$).

D'autres types de preuves de non-finitude différentielle sont présentés page 264.

## 17.2 Solutions polynomiales

Les solutions polynomiales de l'équation différentielle linéaire (17.1) ont des suites de coefficients à support fini solutions de la récurrence linéaire (17.3) qui s'en déduit par traduction, avec les contraintes sur les conditions initiales (17.4). Elles peuvent donc être calculées efficacement par les techniques du chapitre précédent. Les cas inhomogènes et inhomogènes paramétrés avec des membres droits polynomiaux passent à leurs analogues sur les suites dans cette traduction.

## 17.3 Solutions rationnelles

En un pôle *a* d'une solution rationnelle (c'est-à-dire un zéro de son dénominateur), la solution admet un développement en série généralisée de la forme $(t-a)^{\alpha} \sum y_i (t-a)^i$ avec $\alpha \in \mathbb{Z} \setminus \mathbb{N}$. D'après la Proposition 17.5, ces pôles sont donc faciles à localiser : ce sont des racines du coefficient de tête de l'équation différentielle, et leur multiplicité est bornée par l'opposé de la plus petite racine entière négative du polynôme indiciel en ces points.

Ceci conduit à l'Algorithmme 17.1, essentiellement dû à Liouville, pour calculer les solutions rationnelles de (17.1). La preuve de l'algorithme se réduit à observer que le polynôme P qu'il calcule est un multiple du dénominateur de toute solution rationnelle. Du point de vue de la complexité, il faut noter que bien que les exposants $N_a$ puissent avoir une valeur exponentielle en la taille de l'équation différentielle, le polynôme P n'a pas besoin d'être développé, et l'équation obtenue par changement de fonction inconnue grossit indépendamment des $N_a$, puisque la dérivée logarithmique de P est une fraction rationnelle avec un numérateur et un dénominateur de degré au plus *d*.
**Entrée** Une équation différentielle linéaire de la forme
$$c_m(x)y^{(m)}(x) + \cdots + c_0(x)y(x) = 0.$$

**Sortie** Une base de ses solutions rationnelles.

1. En toute racine $a$ de $c_m$ :

- calculer le polynôme indiciel $p_a(n)$;

- calculer la plus petite racine entière négative $N_a$ de $p_a$, s'il n'en existe pas, faire $N_a := 0$.

2. Former le polynôme $P = \prod_{c_m(a)=0} (x-a)^{-N_a}$.

3. Effectuer le changement de fonction inconnue $y = Y/P$ et réduire au même dénominateur.

4. Chercher une base $B$ des solutions polynomiales de cette nouvelle équation. Renvoyer la base des solutions rationnelles formée des fractions $\{b/P, b \in B\}$.

**Algorithme 17.1** – Algorithme de Liouville pour les solutions rationnelles des équations différentielles linéaires.

**Exercice 17.2 — Système et équation.** Ramener la résolution du système

$$Y'(x) = A(x)Y(x), \qquad (17.6)$$

où $A(x)$ est une matrice de fractions rationnelles de $K(x)$ et $Y$ un vecteur, à la résolution d'équations par la méthode de Liouville. ■

**Exercice 17.3** Trouver une équation différentielle linéaire satisfaite par $y_1$ solution de $y'_1 = xy_1 - y_2$, $y'_2 = y_1 - xy_2$. ■

### Le cas inhomogène

Lorsque le membre droit est un polynôme, le polynôme $P$ trouvé par l'Algorithmme 17.1 reste un multiple du dénominateur des solutions rationnelles. En revanche, l'étape de réduction au même dénominateur s'écrit

$$A(x, Y(x)) = P(x) \sum_{i} c_i \frac{d^i Y(x)}{dx^i} = P(x)b(x)$$

et risque de faire exploser le degré des coefficients. Pour éviter cette explosion, il suffit de dériver encore une fois l'équation, et d'effectuer une combinaison linéaire :

$$\begin{align*} & b(x) \frac{d}{dx} A(x, Y(x)) - \left( b(x) \frac{P'(x)}{P(x)} + b'(x) \right) A(x, Y(x)) \\ & \qquad = b(x) (P(x)b(x))' - \left( b(x) \frac{P'(x)}{P(x)} + b'(x) \right) P(x)b(x) = 0. \end{align*}$$

Cette méthode s'applique aussi au cas inhomogène paramétré, et l'équation demeure linéaire en les paramètres.
## 17.4 Séries différentiellement finies analytiques

Le calcul de $\exp(1)$ présenté au chapitre précédent (§15.3) se généralise, jusqu'à l'évaluation numérique en bonne complexité de toutes les séries différentiellement finies convergentes dans une partie de $\mathbb{C}$.

### Évaluation dans le disque de convergence

Soit $y(z)$ une solution de

$$a_r(z)y^{(r)}(z) + \dots + a_0(z)y(z) = 0. \quad (17.7)$$

L'objectif est de calculer numériquement $y$ à partir de conditions initiales. Le point de départ est un résultat d'existence, avec une preuve effective utile pour la suite.

**Proposition 17.6** — Cauchy. Si $z_0 \in \mathbb{C}$ est tel que $a_0$ ne s'annule pas dans le disque $D(z_0, R)$ centré en $z_0$ et de rayon $R$, alors pour tout $(y_0, \dots, y_{r-1}) \in \mathbb{C}^r$, il existe une solution de (17.7) telle que $y^{(i)}(z_0) = y_i$ pour $0 \le i < r$ et qui est analytique dans $D(z_0, R)$.

*Démonstration.* La preuve utilise une méthode de séries majorantes. L'équation (17.7) se récrit matriciellement

$$Y' = AY, \quad Y = \begin{pmatrix} y \\ y' \\ \vdots \\ y^{(r-1)} \end{pmatrix}.$$

Comme le coefficient $a_r$ ne s'annule pas dans $D(z_0, R)$, les coefficients de la matrice A ont un développement en série convergeant dans ce disque. De manière équivalente, la matrice admet un développement en série

$$A(z) = \sum_{k \ge 0} A_k (z - z_0)^k,$$

tel que $0 < \rho < r$, il existe $\alpha > 0$ vérifiant $\|A_k\| \le \alpha/\rho^{k+1}$ pour tout $k \ge 0$. On considère alors la série formelle définie par

$$Y(0) = \begin{pmatrix} y_0 \\ \vdots \\ y_{r-1} \end{pmatrix} \quad \text{et} \quad Y' = AY,$$

c'est-à-dire dont les coefficients sont des vecteurs $Y_n$ qui vérifient

$$(n+1)Y_{n+1} = \sum_{k=0}^{n} A_k Y_{n-k}.$$

Des majorations directes donnent

$$(n+1)\|Y_{n+1}\| \le \sum_{k=0}^{n} \|A_k\| \|Y_{n-k}\| \le \sum_{k=0}^{n} \frac{\alpha}{\rho^{k+1}} \|Y_{n-k}\|.$$
Par récurrence, la suite $(u_n)$ définie par $u_0 = \|Y_0\|$ et

$$(n+1)u_{n+1} = \sum_{k=0}^{n} \frac{\alpha}{\rho^{k+1}} u_{n-k}$$

est telle que $\|Y_n\| \le u_n$ pour tout $n$. Mais la série génératrice $u(z) = \sum u_n z^n$ vérifie

$$u'(z) = \frac{\alpha}{\rho - z} u(z),$$

d'où

$$u(z) = u_0 \left(1 - \frac{z}{\rho}\right)^{-\alpha}.$$

Comme les coefficients de la série Y sont majorés par ceux de $u$, qui converge, il s'ensuit que la série Y converge pour $|z - z_0| < \rho$. Et comme cette propriété est vraie pour tout $\rho < R$, la fonction Y est bien analytique dans $\mathcal{D}(z_0, R)$. ■

Cette preuve fournit aussi des bornes sur les queues du développement en série de $Y(z)$ :

$$\begin{aligned} \left\| \sum_{n \ge N} Y_n (z - z_0)^n \right\| &\le \sum_{n \ge N} u_n |z - z_0|^n, \quad \text{où } u_n = u_0 \rho^{-n} \binom{\alpha + n - 1}{n}, \\ &\le u_0 \rho^{-N} \binom{\alpha + N - 1}{N} |z - z_0|^N \left( 1 + \frac{\alpha + N}{1 + N} \frac{|z - z_0|}{\rho} + \dots \right), \\ &\le C N^{\alpha} \left( \frac{|z - z_0|}{\rho} \right)^N, \end{aligned}$$

pour une certaine constante C. Ainsi, pour garantir

$$\left\| \sum_{n \ge N} Y_n (z - z_0)^n \right\| \le 10^{-K},$$

il suffit de prendre

$$N \log_{10} \left( \frac{|z - z_0|}{\rho} \right) + \alpha \log_{10} N + \log_{10} C \le -K \quad (17.8)$$

ou encore

$$N \ge \frac{K}{\log_{10} \frac{\rho}{|z-z_0|}} + C' \log K \quad (17.9)$$

pour une certaine constante $C'$.

La combinaison de cette borne avec la méthode de scindage binaire du Chapitre 15 fournit donc une évaluation efficace de toutes les séries différentiellement finies à l'intérieur de leur disque de convergence.
**Théorème 17.7** Si $y_i \in \mathbb{Q}$, $a_i \in \mathbb{Q}[z]$, $\zeta \in \mathbb{Q} \cap \mathcal{D}(z_0, \mathbb{R})$ ont tous des numérateurs et dénominateurs bornés par $2^\ell$, alors $y(\zeta)$ peut être calculé à précision $10^{-K}$ en

$$O\left( M_{\mathbb{Z}} \left( \frac{K(\ell + \log K)}{\log \frac{\rho}{|z-z_0|}} \right) \log K \right)$$

opérations binaires, ainsi que $y'(\zeta), \dots, y^{(r-1)}(\zeta)$.

Des informations plus précises sur la constante dans le $O(\cdot)$ et sa dépendance en le degré des coefficients et l'ordre de l'équation différentielle (17.7) s'obtiennent de la même manière sans difficulté autre que la taille des expressions à manipuler.

**Démonstration.** L'algorithme consiste dans un premier temps à obtenir une borne $C$ comme dans le développement ci-dessus, puis à l'utiliser pour borner $N$ grâce à l'équation (17.8). Cette borne croît au plus en $O(K/\log(\rho/|z-z_0|))$ d'après (17.9). Ensuite, il s'agit de calculer la somme du développement en série $y(z)$ au voisinage de $z_0$ jusqu'au $N^e$ terme. Les sommants forment une suite polynomialement récurrente, donc la suite des sommes tronquées aussi. Comme $z_0$ est rationnel, les coefficients de la récurrence satisfaite par les coefficients sont des polynômes de $\mathbb{Q}[z]$. Cette récurrence est régulière par le Corollaire 17.2. La conclusion découle alors du Théorème 15.3. Le même raisonnement s'applique à une dérivée $k^e$, en remplaçant la borne $u_n$ sur les coefficients par $n^k u_n$. À nouveau, l'inégalité (17.9) en découle, avec une constante $C'$ éventuellement plus grande. ■

R L'énoncé ci-dessus est donné pour $\zeta \in \mathbb{Q}$ par simplicité, mais le même raisonnement s'applique à des points complexes du disque de convergence dont la partie réelle et la partie imaginaire sont rationnels. Plus généralement, il s'applique à des nombres algébriques dans le disque de convergence.

## Prolongement analytique

L'équation différentielle (17.7) et ses conditions initiales étant données, il est possible d'étendre le raisonnement précédent pour calculer des valeurs de la solution en dehors du disque de convergence de son développement en série. Si les conditions initiales sont données en $z_0$, le principe est d'utiliser un point intermédiaire $z_1$ à l'intérieur du disque de convergence (borné par la racine la plus proche du coefficient de tête de l'équation). La méthode précédente permet de calculer de nouvelles conditions initiales en $z_1$. On calcule donc ainsi une matrice de transition permettant de passer de conditions initiales en $z_0$ à des conditions initiales en $z_1$. L'intérêt de ce calcul est que le disque centré en $z_1$ et borné par la racine la plus proche du coefficient $a_r$ n'est pas nécessairement inclus dans le disque de départ, et permet donc d'atteindre de nouveaux points, où recommencer le calcul. Ensuite, en bornant les normes des matrices de transition, on peut borner les propagations d'erreurs et ainsi estimer les nombres de termes à prendre à chaque étape.
Doublement de précision

Le Théorème 17.7, combiné au prolongement analytique, montre comment calculer
la valeur d'une fonction différentiellement finie à un point arbitraire de coordonnées
rationnelles de taille modérée. Si le point ζ lui-même est donné à précision 10⁻⁵,
alors l'estimation obtenue par ce théorème avec ℓ = K log₁₀ 2 donne une complexité
quadratique. Il est heureusement possible d'obtenir une meilleure complexité par une
technique de doublement de précision.

**Théorème 17.8** Pour tout point ζ ∈ D(z₀, R), la valeur y(ζ) peut être calculée à précision 10⁻⁴ en

$$
O(K \log^3 K \log \log K)
$$

opérations binaires étant donnés O(K) bits de ζ.

Démonstration. Soient z₀, z₁, ..., z_m = ζ définis avec

$$
z_i := \lceil 2^{2^i} \zeta \rceil 2^{-2^i},
$$

c'est-à-dire que zᵢ est le nombre dyadique ayant les 2ᵢ premiers bits de ζ. (On pourrait
aussi raisonner en base 10, en considérant des nombres décimaux.) Le calcul de y(zᵢ+₁)
utilise donc des entiers plus gros que celui de y(zᵢ), mais plus i est grand, plus zᵢ+₁ - zᵢ
est petit, et donc moins de termes de la série sont nécessaires. Ainsi, dans l'estimation
de complexité du théorème, le terme

$$
\frac{\ell + \log K}{\log \frac{\rho}{|z-z_0|}}
$$

devient

$$
\frac{2^i + \log K}{2^i + \log \rho},
$$

et leur somme est

$$
\sum_{i=1}^{\log K} \left( \frac{2^i}{2^i + \log \rho} + \frac{\log K}{2^i + \log \rho} \right) \le \left( \sum_{i=1}^{\log K} 1 \right) + \sum_{i=1}^{\log K} \frac{\log K}{2^i} = O(\log K).
$$

La formule finale est obtenue en utilisant pour simplifier la complexité de la multipli-
cation d'entiers par FFT.
■

Exercice

Cet exercice est le pendant pour les équations différentielles de l'Exercice 16.6 du chapitre précédent sur l'utilisation de formes compactes pour le calcul efficace.

**Exercice 17.4** (Forme compacte dans la base monomiale.)

1. Si le polynôme P est donné en forme compacte dans la base (X^k), c'est-à-dire au moyen d'une récurrence sur ses coefficients, et de conditions initiales généralisées, montrer comment obtenir efficacement son pgcd avec X^N.
2. Montrer que si $(u_k)$ est solution d'une récurrence linéaire de la forme (16.1) et $\alpha \in \mathbb{Q} \setminus \{0\}$, alors la suite $(u_k \alpha^k)$ est solution d'une récurrence linéaire de la même forme, qui se calcule efficacement, dont le coefficient de tête a les mêmes racines entières positives que $p_0$, et dont les conditions initiales généralisées se calculent efficacement à partir de celles de $(u_k)$.

3. Pour la suite $s_k := \sum_{i=0}^{k} u_i a^i$, montrer qu'on a encore facilement une récurrence linéaire, ainsi que les r premières conditions initiales, et que les racines du coefficient de tête sont encore inchangées.

4. Si $(u_k)$ est donnée sous forme compacte, montrer comment obtenir les conditions initiales généralisées pour cette suite $(s_k)$.

5. Montrer que si P est un polynôme donné sous forme compacte dans la base monomiale alors $P(\alpha)$ peut être calculé en un nombre d'opérations binaires quasi-linéaire en N, pour a comme ci-dessus.

6. Même question pour la dérivée $P^{(\ell)}(\alpha)$, pour $\ell = O(N)$.

7. Montrer comment obtenir les mêmes résultats pour $(\sum_{i=0}^{k} u_i \alpha^i)$ lorsque $\alpha$ est un nombre algébrique, donné par un polynôme $Q \in \mathbb{Q}[X]$ tel que $Q(\alpha) = 0$. (On admettra que les racines du terme de tête sont inchangées.)

8. En déduire que si P est donné en forme compacte dans la base $(X^k)$, on peut en obtenir une forme compacte dans la base $((X - \alpha)^k)$ pour $\alpha$ donné par un polynôme $Q \in \mathbb{Q}[X]$ tel que $Q(\alpha) = 0$.

9. Décrire un algorithme permettant de diviser efficacement le numérateur obtenu sous forme compacte dans l'algorithme de Liouville par son pgcd avec le polynôme P multiple du dénominateur. Estimer sa complexité. ■

## Notes

Les idées de base de l'algorithme de recherche de solutions rationnelles sont dues à Liouville [Lio33] qui donne également une méthode par coefficients indéterminés pour trouver les solutions polynomiales. La présentation qui utilise les récurrences donne un algorithme de même complexité mais fait mieux ressortir la structure du calcul [ABP95b]. Une bonne introduction aux équations différentielles dans le plan complexe est l'ouvrage classique de Ince [Inc56].

La recherche de solutions d'équations différentielles linéaires ne s'arrête pas aux solutions rationnelles. En utilisant la théorie de Galois différentielle, on peut concevoir une algorithmique sophistiquée de recherche de solutions *liouvilliennes* (c'est-à-dire formées par l'application répétée d'exponentielles, d'intégrales et de prise de racines de polynômes). Les calculs se ramènent à la recherche présentée ici de solutions rationnelles pour des équations (les puissances symétriques) formées à partir de l'équation de départ [Mar98; Sin81; SU97].

L'une des premières références mentionnant le scindage binaire [BGS72, §178] suggère déjà (sans preuve) qu'il permet l'évaluation numérique rapide des fonctions différentiellement finies. L'article de synthèse de Bernstein [Ber08, §12] est une bonne référence sur l'historique de la méthode.

La méthode de scindage binaire est un grand classique en calcul formel. Elle a été redécouverte à plusieurs reprises dans la littérature, dans différents contextes :
conversion d'une base à l'autre, calcul de décimales de constantes comme $\pi$ ou $e$ [Bre76], calcul avec des récurrences linéaires [KS73], [CC88, §6]. En particulier, c'est ainsi que les systèmes de calcul formel calculent rapidement $\pi$ par une formule découverte en 1989 par les frères Chudnovsky :

$$ \frac{1}{\pi} = \frac{12}{C^{3/2}} \sum_{n=0}^{\infty} \frac{(-1)^n (6n)!(A+nB)}{(3n)!n!^3 C^{3n}} $$

où $A = 13591409$, $B = 545140134$ et $C = 640320$. Cette série donne environ 14 décimales nouvelles par terme. Ceci n'est pas suffisant pour parvenir à une bonne efficacité, que l'on obtient en observant que les sommants vérifient une récurrence linéaire d'ordre 1 et en appliquant la méthode de scindage binaire.

Les algorithmes numériques présentés dans les dernières sections requièrent le calcul de bornes explicites pour le passage à l'implantation. C'est un travail un peu technique pour lequel nous renvoyons à la littérature [MS10].

## Bibliographie

ABP95b ABRAMOV, Sergei A., Manuel BRONSTEIN et Marko PETKOVŠEK (1995). « On polynomial solutions of linear operator equations ». In : *ISSAC'95 : International Symposium on Symbolic and Algebraic Computation*. Éd. par A. H. M. LEVELT. New York : ACM Press, p. 290–296.

Ber08 BERNSTEIN, Daniel J. (2008). « Fast multiplication and its applications ». In : *Algorithmic number theory : lattices, number fields, curves and cryptography*. Vol. 44. Publications of the Research Institute for Mathematical Sciences. Cambridge University Press, p. 325–384.

BGS72 BEELER, Michael, R. William GOSPER et Richard SCHROEPPEL (1972). *HAK-MEM*. Artificial Intelligence Memo No. 239. Massachusetts Institute of Technology.

Bre76 BRENT, Richard P. (1976). « Multiple-precision zero-finding methods and the complexity of elementary function evaluation ». In : *Analytic computational complexity*. Proceedings of a Symposium held at Carnegie-Mellon University, Pittsburgh, PA, 1975. Academic Press, p. 151–176.

CC88 CHUDNOVSKY, D. V. et G. V. CHUDNOVSKY (1988). « Approximations and complex multiplication according to Ramanujan ». In : *Proceedings of the Centenary Conference, University of Illinois at Urbana-Champaign, June 1–5, 1987*. Éd. par G. E. ANDREWS, R. A. ASKEY, B. C. BERNDT, K. G. RAMANATHAN et R. A. RANKIN. Academic Press, p. 375–472.

Inc56 INCE, E. L. (1956). *Ordinary differential equations*. Reprint of the 1926 edition. Dover Publications.

KS73 KOGGE, P. et H. STONE (1973). « A parallel algorithm for the efficient solution of a general class of recurrence equations ». In : *IEEE Transactions on Computers*, vol. C-22, p. 786–793.

Lio33 LIOUVILLE, Joseph (1833). « Second mémoire sur la détermination des intégrales dont la valeur est algébrique ». In : *Journal de l'École polytechnique*, vol. 14, p. 149–193.
Mar98 MAROTTE, F. (1898). « Les équations différentielles linéaires et la théorie des groupes ». In : *Annales de la Faculté des Sciences de Toulouse*, vol. 12, n°4, p. 1-92.

MS10 MEZZAROBBA, Marc et Bruno SALVY (2010). « Effective Bounds for P-Recursive Sequences ». In : *Journal of Symbolic Computation*, vol. 45, n°10, p. 1075-1096.

Sin81 SINGER, Michael F. (1981). « Liouvillian solutions of *n*-th order homogeneous linear differential equations ». In : *American Journal of Mathematics*, vol. 103, n°4, p. 661-682.

SU97 SINGER, Michael F. et Felix ULMER (1997). « Linear differential equations and products of linear forms ». In : *Journal of Pure and Applied Algebra*, vol. 117/118, p. 549-563.
<figure><img src="image_50.png" /></figure>
# IV Factorisation des polynômes

<table>
  <tbody>
    <tr>
      <td>18</td>
      <td>Factorisations sans carré et séparable</td>
      <td>327</td>
    </tr>
    <tr>
      <td>18.1</td>
      <td>Introduction</td>
      <td></td>
    </tr>
    <tr>
      <td>18.2</td>
      <td>Factorisation sans carré</td>
      <td></td>
    </tr>
    <tr>
      <td>18.3</td>
      <td>Séparabilité et déflation</td>
      <td></td>
    </tr>
    <tr>
      <td>18.4</td>
      <td>Factorisation séparable</td>
      <td></td>
    </tr>
    <tr>
      <td>18.5</td>
      <td>Réduction à la factorisation des polynômes séparables</td>
      <td></td>
    </tr>
    <tr>
      <td>19</td>
      <td>Factorisation des polynômes à une variable sur un corps fini</td>
      <td>345</td>
    </tr>
    <tr>
      <td>19.1</td>
      <td>Corps finis, quelques rappels</td>
      <td></td>
    </tr>
    <tr>
      <td>19.2</td>
      <td>Algorithme de Berlekamp</td>
      <td></td>
    </tr>
    <tr>
      <td>19.3</td>
      <td>Algorithme de Cantor et Zassenhaus</td>
      <td></td>
    </tr>
    <tr>
      <td>20</td>
      <td>Réduction de réseaux et algorithme LLL</td>
      <td>359</td>
    </tr>
    <tr>
      <td>20.1</td>
      <td>Réseaux et vecteurs courts</td>
      <td></td>
    </tr>
    <tr>
      <td>20.2</td>
      <td>Le procédé d'orthogonalisation de Gram-Schmidt</td>
      <td></td>
    </tr>
    <tr>
      <td>20.3</td>
      <td>L'algorithme LLL</td>
      <td></td>
    </tr>
    <tr>
      <td>20.4</td>
      <td>Preuve de l'algorithme LLL</td>
      <td></td>
    </tr>
    <tr>
      <td>21</td>
      <td>Factorisation des polynômes à une variable sur les rationnels</td>
      <td>373</td>
    </tr>
    <tr>
      <td>21.1</td>
      <td>Bornes sur les coefficients des facteurs</td>
      <td></td>
    </tr>
    <tr>
      <td>21.2</td>
      <td>Factorisation sans carré</td>
      <td></td>
    </tr>
    <tr>
      <td>21.3</td>
      <td>Algorithme par remontée et recombinaison</td>
      <td></td>
    </tr>
    <tr>
      <td>21.4</td>
      <td>Algorithme de van Hoeij</td>
      <td></td>
    </tr>
    <tr>
      <td>22</td>
      <td>Factorisation des polynômes à plusieurs variables</td>
      <td>399</td>
    </tr>
    <tr>
      <td>22.1</td>
      <td>Algorithme naïf par remontée</td>
      <td></td>
    </tr>
    <tr>
      <td>22.2</td>
      <td>Recombinaison des dérivées logarithmiques</td>
      <td></td>
    </tr>
    <tr>
      <td>22.3</td>
      <td>Algorithme de Lecerf</td>
      <td></td>
    </tr>
    <tr>
      <td>22.4</td>
      <td>Réduction probabiliste au cas de deux variables</td>
      <td></td>
    </tr>
    <tr>
      <td>22.5</td>
      <td>Cas particulier des corps finis</td>
      <td></td>
    </tr>
  </tbody>
</table>
<figure><img src="image_51.png" /></figure>
# 18. Factorisations sans carré et séparable

## Résumé

Les algorithmes de factorisation des polynômes sont constitués de plusieurs étapes de natures bien différentes. Pour des polynômes à une variable, la première d'entre elles consiste à décomposer le polynôme en facteurs dont les racines sont de même multiplicité. On parle alors de factorisation séparable. Lorsque les coefficients appartiennent à un corps de caractéristique zéro ou plus grande que le degré, alors la factorisation séparable coïncide avec la factorisation sans carré. Le présent chapitre présente des algorithmes pour ces tâches.

## 18.1 Introduction

Soit $\mathbb{K}$ un corps, et soit $F$ un polynôme unitaire de $\mathbb{K}[X]$. La factorisation en irréductibles de $F$, notée $\text{firr}(F) = \{(F_1, m_1), \dots, (F_r, m_r)\}$ avec les $F_i$ irréductibles unitaires, correspond à l'identité

$$F = F_1^{m_1} \cdots F_r^{m_r}.$$

L'entier $m_i$ est appelé la *multiplicité du facteur* $F_i$.

On dit que $F$ est *sans carré* s'il n'existe pas de polynôme $P$ de $\mathbb{K}[X] \setminus \mathbb{K}$ tel que $P^2$ divise $F$. Autrement dit, toutes les multiplicités des facteurs irréductibles de $F$ sont égales à 1. Si l'on note $\{\mu_1, \dots, \mu_s\}$ l'ensemble des multiplicités des facteurs irréductibles (avec $s \le r$), et $Q_i$ le produit des $F_i$ de multiplicité $\mu_i = m_i$, alors la factorisation sans carré de $F$, notée $\text{fsc}(F)$, est l'ensemble des paires $(Q, m)$ où $Q$ est le produit des facteurs irréductibles unitaires de $F$ de multiplicité exactement $m \ge 1$. La partie sans carré de $F$ est alors définie comme le produit $\prod_{(Q,m) \in \text{fsc}(F)} Q$; il s'agit manifestement d'un polynôme sans carré.
Par exemple, si $K = Q$, et

$$F = (X + 1)(X + 2)(X + 3)^2(X + 4)^2(X + 5)^3(X + 6)^3,$$

sa partie sans carré est le produit de ses facteurs irréductibles, c'est-à-dire

$$(X + 1)(X + 2)(X + 3)(X + 4)(X + 5)(X + 6).$$

La factorisation sans carré de F est ici :

$$fsc(F) = \{((X + 1)(X + 2), 1), ((X + 3)(X + 4), 2), ((X + 5)(X + 6), 3)\}.$$

On voit bien que le calcul de factorisation sans carré est l'amorce de la factorisation irréductible; en outre, c'est un calcul sensiblement plus facile. En effet, nous verrons plus loin un algorithme de calcul de la partie sans carré, de complexité quasi-linéaire en le degré, qui repose uniquement sur des calculs de pgcd bien choisis lorsque $K$ est de caractéristique zéro ou suffisamment grande.

D'une façon générale, la factorisation irréductible d'un polynôme est l'une des premières questions difficiles du calcul formel. Nous présenterons ces difficultés et les solutions connues pour traiter les cas usuels, rencontrés en pratique, dans le Chapitre 22. Dans un premier temps nous exposerons les idées classiques en caractéristique zéro pour calculer la factorisation sans carré. En caractéristique non nulle la situation est nettement plus complexe et il existe des corps effectifs où l'extraction des racines $p$-ièmes n'est pas calculable. La factorisation sans carré n'est alors pas calculable en général et nous verrons alors l'intérêt d'introduire la factorisation séparable, qui nous sera utile dans les chapitres suivants.

Tout au long de ce chapitre, $K$ représente un corps effectif, $p$ sa caractéristique, et $F$ un polynôme unitaire de $K[X]$ de degré $n \ge 1$.

## 18.2 Factorisation sans carré

Dans cette section nous considérons le cas de la caractéristique zéro, et nous nous intéressons d'abord au calcul de la partie sans carré de $F$. Un algorithme pour la factorisation sans carré s'en déduit : on divise $F$ par sa partie sans carré, on calcule la partie sans carré du résultat, etc. Ce faisant, on n'obtient pas une borne de complexité quasi-linéaire. L'algorithme rapide, dû à Yun, est présenté en fin de section.

### Partie sans carré

Rappelons, qu'en caractéristique zéro un polynôme $F$ de $K[X]$ est sans carré si, et seulement si, $pgcd(F, F') = 1$. Commençons par le lemme suivant, utilisé plusieurs fois dans la suite :

**Lemme 18.1** Supposons $K$ de caractéristique $p = 0$. Soient $G_1, \dots, G_s$ des polynômes sans carré de $K[X] \setminus K$, premiers deux à deux, soit $G = G_1 \cdots G_s$, et soit
$(c_1, \dots, c_s) \in \mathbb{K}^s$. Alors nous avons l'égalité suivante :

$$ \text{pgcd}\left(G, \sum_{i=1}^{s} c_i \frac{G}{G_i} G_i'\right) = \prod_{c_i=0} G_i. $$

*Démonstration.* La preuve découle du calcul suivant :

$$ \text{pgcd}\left(G, \sum_{i=1}^{s} c_i \frac{G}{G_i} G_i'\right) = \prod_{j=1}^{s} \text{pgcd}\left(G_j, \sum_{i=1}^{s} c_i \frac{G}{G_i} G_i'\right) = \prod_{j=1}^{s} \text{pgcd}\left(G_j, c_j \frac{G}{G_j} G_j'\right). $$

Ce dernier pgcd est égal à $G_j$ si $c_j = 0$ et 1 si $c_j \neq 0$ puisque $G_j'$ est premier avec $G_j$. ■

**Proposition 18.2** Si $K$ est de caractéristique zéro, la partie sans carré de $F$ est donnée par la formule $F/\text{pgcd}(F, F')$ et peut être calculée au moyen de $O(M(n) \log n)$ opérations dans $K$.

*Démonstration.* Écrivons la factorisation en irréductibles unitaires de $F$ :

$$ F = F_1^{m_1} \cdots F_r^{m_r}, $$

où tous les $m_i$ sont des entiers $\geq 1$. En dérivant, on obtient

$$ F' = \sum_{i=1}^{r} m_i F_i' \frac{F}{F_i}. $$

En posant $E = \prod_{i=1}^{r} F_i$ la partie sans carré de $F$, on en déduit que

$$ \text{pgcd}(F, F') = \text{pgcd}\left(\prod_{i=1}^{r} F_i^{m_i}, \sum_{i=1}^{r} m_i \frac{F}{F_i} F_i'\right) = \text{pgcd}\left(E, \sum_{i=1}^{r} m_i \frac{E}{F_i} F_i'\right) \prod_{i=1}^{r} F_i^{m_i-1} = \prod_{i=1}^{r} F_i^{m_i-1}, $$

où la dernière égalité provient du Lemme 18.1. La partie sans carré $E$ de $F$ s'obtient alors comme $F/\text{pgcd}(F, F')$. Le coût du calcul est essentiellement une conséquence de la Proposition 6.14 qui donne le coût du pgcd. ■

## Algorithme de Musser

La première méthode de factorisation sans carré, que nous présentons dans l'Algorithmme 18.1, commence par calculer la partie sans carré.

**Proposition 18.3** L'Algorithmme 18.1 de factorisation sans carré est correct et son coût arithmétique est en $O(nM(n) \log n)$.

*Démonstration.* Notons $(G_1, m_1), \dots, (G_s, m_s)$ la factorisation sans carré de $F$ normalisée de sorte que

$$ F' = \sum_{i=1}^{s} m_i G_i' G_i^{m_i-1} \frac{F}{G_i^{m_i}}. $$
**Entrée** F dans $K[X]$, avec $K$ de caractéristique zéro.

**Sortie** fsc(F).

1. Initialiser L avec la liste vide et $i$ à l'entier 1.

2. Calculer $C = \text{pgcd}(F, F')$ et $G = F/C$.

3. Tant que $\text{deg}(G) \ge 1$ faire :

a. calculer $P = \text{pgcd}(C, G)$ et $C = C/P$;

b. si $\text{deg}(G) > \text{deg}(P)$ alors adjoindre $(G/P, i)$ à $L$;

c. $G = P$, et $i = i + 1$.

4. Renvoyer L.

**Algorithme 18.1** – Algorithme de Musser en caractéristique zéro.

**Entrée** F dans $K[X]$, avec $K$ de caractéristique zéro.

**Sortie** fsc(F).

1. Poser $l = 1$ et initialiser L avec la liste vide.

2. Calculer $U = \text{pgcd}(F, F')$, $V = F/U$, et $W = F'/U$.

3. Tant que $\text{deg}(V) \ge 1$ faire :

a. calculer $H = \text{pgcd}(V, W - V')$, $W = (W - V')/H$, et $V = V/H$;

b. si $\text{deg}(H) \ge 1$ alors adjoindre $(H, l)$ à $L$;

c. incrémenter $l$ de 1.

4. Renvoyer L.

**Algorithme 18.2** – Algorithme de Yun en caractéristique zéro.

Alors, à l'étape (2), nous avons $C = \prod_{i=1}^{s} G_i^{m_i-1}$, et $G = \prod_{i=1}^{s} G_i$. Dans la première étape de la boucle nous obtenons $P = \prod_{i=1, m_i \ge 2}^{s} G_i$, $C = \prod_{i=1, m_i \ge 2}^{s} G_i^{m_i-2}$, et $G/P$ correspond aux $G_i$ pour lesquels $m_i = 1$. À la deuxième étape de la boucle nous avons

$$P = \prod_{i=1, m_i \ge 3}^{s} G_i, \quad C = \prod_{i=1, m_i \ge 3}^{s} G_i^{m_i-3},$$

et $G/P$ correspond aux $G_i$ pour lesquels $m_i = 2$. Par récurrence nous déduisons que l'algorithme est correct.

Le nombre d'étapes de l'algorithme est borné par $n$, le coût de chaque étape est dominé par celui du pgcd qui est, rappelons-le, en $O(M(n) \log n)$. $\blacksquare$

## Algorithme de Yun

La méthode suivante de factorisation sans carré, présentée dans l'Algorithme 18.2, permet d'atteindre un coût arithmétique essentiellement linéaire.

**Proposition 18.4** L'Algorithme 18.2 de factorisation sans carré est correct et son coût arithmétique est en $O(M(n) \log n)$ opérations dans $K$.
Démonstration. Notons $V_l$ et $W_l$ les valeurs respectives de V et W au départ de l'étape $l$ de la boucle « tant que », et notons $H_l$ la valeur de H calculée durant l'étape $l$. Notons $\ell$ la valeur de $l$ à la sortie de la boucle. Définissons aussi $V_\ell$ et $W_\ell$ comme les valeurs respectives de V et W à la sortie de la boucle. Nous allons prouver la propriété suivante par récurrence sur $l$ de 1 à $\ell$ :

$$V_l = \prod_{(G,m) \in \text{fsc}(F) \atop m \ge l} G, \quad W_l = \sum_{(G,m) \in \text{fsc}(F) \atop m \ge l} (m-l+1) \frac{V_l}{G} G'$$

À l'étape (2) nous avons :

$$U = \text{pgcd}(F, F'), \quad V = \frac{F}{\text{pgcd}(F, F')}, \quad W = \frac{F'}{\text{pgcd}(F, F')}.$$

Comme nous avons vu précédemment que $\text{pgcd}(F, F') = \prod_{(G,m) \in \text{fsc}(F)} G^{m-1}$, l'hypothèse de récurrence est satisfaite pour $l=1$. Supposons l'hypothèse de récurrence satisfaite pour une certaine valeur de $l$ dans $\{1, \dots, \ell-1\}$. En utilisant le Lemme 18.1 nous déduisons :

$$H_l = \text{pgcd}(V_l, W_l - V_l') = \text{pgcd} \left( \prod_{(G,m) \in \text{fsc}(F) \atop m \ge l} G, \sum_{(G,m) \in \text{fsc}(F) \atop m \ge l} (m-l) \frac{V_l}{G} G' \right) = \prod_{(G,m) \in \text{fsc}(F) \atop m=l} G,$$

ce qui conduit aux formules attendues pour $V_{l+1}$ et $W_{l+1}$. Puisque $\ell$ est le premier entier tel que $\text{deg}(V_\ell) = 0$, nous obtenons que $\ell-1$ est la plus grande des multiplicités des facteurs de F. Finalement L contient tous les facteurs sans carré de F à la fin de la boucle.

L'étape (2) utilise $O(M(n) \log n)$ opérations dans $\mathbb{K}$. L'étape $l$ de la boucle coûte au plus $O(M(\text{deg}(V_l)) \log(\text{deg}(V_l)))$ puisque $\text{deg}(W_l) \le \text{deg}(V_l) - 1$. Par la super-additivité de M, l'étape (3) totalise au plus $O(M(\sum_{l=1}^{\ell-1} \text{deg}(V_l)) \log n)$ opérations dans $\mathbb{K}$, ce qui conduit à la borne de complexité attendue puisque $\prod_{l=1}^{\ell-1} V_l = F$. ■

## 18.3 Séparabilité et déflation

Rappelons que le corps fini à $q$ éléments est noté $\mathbb{F}_q$. Si $p$ est premier alors $\mathbb{F}_p$ n'est autre que $\mathbb{Z}/p\mathbb{Z}$. Les rappels mathématiques sur les corps finis sont regroupés au début du chapitre suivant.

Si $\mathbb{K}$ est le corps des fractions $\mathbb{F}_2(T)$, et si $F(X) = T^2 + X^2 = (T+X)^2$, alors $F'(X) = 0$, et ni l'algorithme de Musser ni celui de Yun ne calcule la factorisation sans carré de F. Nous allons voir que le « bon concept » en caractéristique non nulle n'est pas celui d'être sans carré, mais plutôt d'être séparable :

**Définition 18.1** Un polynôme $F \in \mathbb{K}[X]$ est dit séparable s'il n'a aucune racine multiple dans une clôture algébrique $\bar{\mathbb{K}}$ de $\mathbb{K}$.
Nous rappelons la définition du discriminant, déjà abordée dans l'Exemple 6.7 page 118.

**Définition 18.2** Le discriminant de $F \in \mathbb{K}[X]$ de degré $n$, noté disc$(F)$, est défini par

$$ \mathrm{Res}(F, F') = (-1)^{n(n-1)/2} \mathrm{ct}(F) \mathrm{disc}(F), $$

où ct$(F)$ représente le coefficient de degré $n$ de $F$.

**Proposition 18.5** Un polynôme $F$ qui n'est pas constant est séparable si, et seulement si, son discriminant est non nul.

*Démonstration.* Notons $\alpha_1, \dots, \alpha_r$ les racines distinctes de $F$ dans $\bar{\mathbb{K}}$ de multiplicités respectives $m_1, \dots, m_r$. On peut supposer $F$ unitaire. Si toutes les multiplicités valent 1, alors $F'(\alpha_i) = \prod_{j \neq i}(\alpha_i - \alpha_j) \neq 0$. Par la formule de Poisson (Théorème 6.5), le discriminant est égal au produit des $F'(\alpha_i)$, il est donc non nul. Si l'une des multiplicités, disons $m_1 \ge 2$ alors $F'(\alpha_1) = 0$ et le discriminant s'annule. ■

Par exemple, si $\mathbb{K}$ est le corps des fractions $\mathbb{F}_2(T)$, alors $F(X) = T + X^2$ est sans carré mais n'est pas séparable pour autant puisque $F'(X)$ est identiquement nul. Néanmoins si $p = 0$, alors la condition de séparabilité est équivalente à être sans carré.

Dans la suite nous noterons $\mathcal{B} = \{1, p, p^2, p^3, \dots\}$ l'ensemble des puissances de $p$ si $p > 0$, et $\mathcal{B} = \{1\}$ si $p = 0$. Les problèmes liés à la séparabilité concernent les cas de caractéristique non nulle. Nous aurons besoin des définitions et propositions suivantes.

**Définition 18.3** Le degré d'inséparabilité de $F$, noté $\mathrm{deg}^i(F)$, est le plus grand élément $q$ de $\mathcal{B}$ tel que $F \in \mathbb{K}[X^q] \setminus \mathbb{K}[X^{pq}]$.

**Proposition 18.6** Soient $F$ et $G$ deux polynômes de $\mathbb{K}[X] \setminus \mathbb{K}$.

1. Le degré d'inséparabilité de $F$ est la plus grande puissance de $p$ qui divise la multiplicité de chaque racine de $F$.

2. $\mathrm{deg}^i(\mathrm{FG}) \ge \min(\mathrm{deg}^i(F), \mathrm{deg}^i(\mathrm{G}))$, l'inégalité devient une égalité dès lors que $F$ et $G$ sont premiers entre eux.

*Démonstration.* En utilisant $(X - a)^q = X^q - a^q$ pour tout $a$ dans une clôture algébrique de $\mathbb{K}$, si $q \in \mathcal{B}$ divise la multiplicité de chaque racine de $F$, alors $F$ appartient bien à $\mathbb{K}[X^q]$. Réciproquement, si $F$ peut s'écrire $G(X^q)$ avec $q \in \mathcal{B}$, alors la multiplicité d'une racine de $F$ est bien un multiple de $q$. La partie (2) est une conséquence directe de la partie (1). ■

**Définition 18.4** Le degré de séparabilité de $F \in \mathbb{K}[X]$, noté $\mathrm{deg}^s(F)$, est défini comme
suit :

$$ \deg^s(F) = \sum_{(G,m) \in \text{firr}(F)} \deg(G)/\deg^i(G) = \sum_{(G,m) \in \text{firr}(F)} \deg^s(G). $$

**Exemple 18.1** Soient $K = \mathbb{F}_2(T)$ et $F(X) = X^4 + X^2 + 1 + T$. On a $\deg^i(F) = 2$, et comme $F$ est irréductible, $\deg^s(F) = 2$.

**Proposition 18.7** Soient $F$ et $G$ des polynômes de $K[X] \setminus K$.

1. $\deg^s(F)$ est égal au nombre de racines de $F$ dans $\tilde{K}$ sans compter les multipli-cités.

2. $F$ est séparable si, et seulement si, $\deg^s(F) = \deg(F)$.

3. Pour tout $q \in \mathcal{B}$, $\deg^s(F(X^q)) = \deg^s(F)$.

4. $\deg^s(FG) \le \deg^s(F) + \deg^s(G)$, avec égalité si, et seulement si, $F$ et $G$ sont premiers entre eux.

*Démonstration.* Ces propriétés découlent presque immédiatement des définitions et nous laissons au lecteur le soin de les vérifier à titre d'exercice. ■

Lorsque qu'un polynôme $F$ appartient à $K[X^p]$, il est utile de le voir comme un polynôme en la variable $X^p$. Ceci motive la définition suivante :

**Définition 18.5** Si $F \in K[X] \setminus K$, le polynôme déflaté de $F$ est l'unique polynôme $\tilde{F}$ défini par

$$ \tilde{F}(X^{\deg^i(F)}) = F(X). $$

Les multiplicités des racines de $\tilde{F}$ sont celles de $F$ divisées par $\deg^i(F)$, d'où la termi-nologie « déflation ». Puisque $\tilde{F}$ appartient à $K[X] \setminus K[X^p]$, la fonction suivante $\Phi$ de déflation est bien définie :

$$ \Phi : K[X] \setminus K \to (K[X] \setminus K[X^p]) \times \mathcal{B} \\ F \mapsto (\tilde{F}, \deg^i(F)). $$

**Exemple 18.2** Avec l'exemple précédent où $K = \mathbb{F}_2(T)$ et $F(X) = X^4 + X^2 + 1 + T$, on a $\Phi(F) = (X^2 + X + 1 + T, 2)$.

Notons $\mathcal{Q}$ l'ensemble des puissances $q$-ièmes des polynômes irréductibles de $K[X] \setminus K$ pour $q \in \mathcal{B}$, et $\mathcal{S}$ le sous-ensemble des polynômes séparables irréductibles.

**Proposition 18.8** $\Phi$ est une bijection qui envoie $\mathcal{Q}$ surjectivement sur $\mathcal{S} \times \mathcal{B}$.

*Démonstration.* Le fait que $\Phi$ est une bijection découle clairement des définitions. Si $F \in K[X]$ est irréductible, alors $\tilde{F}$ est nécessairement irréductible et donc séparable
car $\tilde{F} \in \mathbb{K}[X^p]$. Puisque le polynôme déflaté de $F^q$ coincide avec celui de $\tilde{F}^q$ pour tout $q \in \mathcal{B}$, le Lemme 18.9 ci-dessous implique que $\Phi(\mathcal{Q})$ appartient bien à $\mathcal{S} \times \mathcal{B}$.

Réciproquement, si $(G, q) \in \mathcal{S} \times \mathcal{B}$. Nous devons montrer que $F(X) = G(X^q)$ appartient à $\mathcal{Q}$, et que $\deg^i(F) = q$. Cette dernière égalité est claire puisque $G$ est séparable. Soit $h \le q$ le plus grand élément de $\mathcal{B}$ tel que $G(X^h)$ est une puissance $h$-ième, et soit $\tilde{G}$ la racine $h$-ième correspondante, de sorte que $\tilde{G}^h(X) = G(X^h)$. Par le Lemme 18.9 ci-dessous, $\tilde{G}$ et donc $\tilde{G}(X^{q/h})$ sont irréductibles et séparables, ce qui prouve que $F(X) = \tilde{G}(X^{q/h})^h \in \mathcal{Q}$. ■

**Lemme 18.9** Soient $F$ et $G$ deux polynômes de $\mathbb{K}[X] \setminus \mathbb{K}$ tels que $G(X^p) = F(X)^p$.

1. $G$ est séparable si, et seulement si, $F$ est séparable.

2. Si $G$ est irréductible alors $F$ est irréductible. La réciproque est vraie si $F$ ou $G$ est séparable.

*Démonstration.* La preuve de la partie (a) est immédiate. Les multiplicités de $G(X^p)$ (resp. de $F(X)^p$) sont toutes $p$ si, et seulement si, $G$ (resp. $F$) est séparable. La partie (b) est donc une conséquence de la Proposition 18.7. Concernant la partie (c), si $G$ est irréductible, alors pour n'importe quel diviseur strict $A$ de $F$, le polynôme $A^p$ divise $F^p$, et donc $B$ divise $G$, si nous définissons $B$ par $B(X^p) = A(X)^p$. Par conséquent $B$ et donc $A$ sont dans $\mathbb{K}$. Réciproquement, supposons $F$ irréductible, et soit $A$ un diviseur strict de $G$. Puisque $A(X^p)$ divise $F(X)^p$, le polynôme $A(X^p)$ peut s'écrire $F(X)^\alpha$ pour un certain $\alpha \in \{0, \dots, p-1\}$. En dérivant, nous obtenons $\alpha F'F^{\alpha-1} = 0$, d'où $\alpha = 0$ dès que $F$ est séparable. ■

## 18.4 Factorisation séparable

Nous pouvons maintenant définir un nouveau type de factorisation correspondant à écrire $F$ essentiellement comme un produit de facteurs séparables :

**Définition 18.6** La factorisation séparable d'un polynôme $F \in \mathbb{K}[X] \setminus \mathbb{K}$, notée fsep(F), est l'ensemble des triplets $(G_1, q_1, m_1), \dots, (G_s, q_s, m_s)$ de $\mathbb{K}[X] \times \mathcal{B} \times \mathbb{N}^*$ vérifiant les propriétés suivantes :

$$ (S_1) \quad F(X) = \prod_{i=1}^{s} G_i^{m_i}(X^{q_i}), $$

(S2) les $G_i(X^{q_i})$ sont premiers entre eux deux à deux,

(S3) les $m_i$ ne sont pas divisibles par $p$,

(S4) les $G_i$ sont séparables de degré au moins 1,

(S5) les couples $(q_i, m_i)$ sont deux à deux distincts.

Notons que dans le cas où $p = 0$ la factorisation séparable n'est rien d'autre que la factorisation sans carré.

**Théorème 18.10** Tout polynôme $F$ de $\mathbb{K}[X]$ admet une unique factorisation séparable (à permutation et unités près dans $\mathbb{K}$).
Démonstration. Soient $(F_1, e_1), \dots, (F_r, e_r)$ la factorisation irréductible de $F$ dans $K[X]$. Pour chaque $i \in \{1, \dots, r\}$, soit $a_i$ la plus grande puissance de $p$ qui divise $e_i$, et soit $m_i = e_i/a_i$. En définissant $(G_i, q_i) = \Phi(F_i^{a_i})$, nous obtenons des triplets

$$(G_1, q_1, m_1), \dots, (G_s, q_s, m_s)$$

qui satisfont $(S_1), (S_2)$ et $(S_3)$. La propriété $(S_4)$ est aussi satisfaite par la Proposition 18.8.

Pour obtenir $(S_5)$ il suffit de regrouper les triplets qui partagent les mêmes deuxièmes et troisièmes coordonnées comme suit : si $(A_1, q, m), \dots, (A_t, q, m)$ sont de tels triplets alors nous les fusionnons en le triplet $(A_1 \cdots A_r, q, m)$. Les propriétés de $(S_1)$ à $(S_4)$ sont préservées, ce qui prouve l'existence de la factorisation séparable.

Supposons maintenant que nous disposons de triplets

$$(G_1, q_1, m_1), \dots, (G_s, q_s, m_s)$$

qui satisfont les propriétés de $(S_1)$ à $(S_5)$. Pour chaque $i \in \{1, \dots, s\}$, n'importe quelle racine de $G_i(X^{q_i})$ dans $\bar{K}$ admet $q_i$ comme multiplicité, et donc $q_i m_i$ comme multiplicité dans $F$. La propriété $(S_5)$ implique alors que les racines de $G_i(X^{q_i})$ sont exactement celles de $F$ de multiplicité $q_i m_i$, ce qui conduit à l'unicité de la factorisation séparable. ■

**Exemple 18.3** Lorsque $K = \mathbb{F}_3$ et $F = X^2(X+1)^3(X+2)^4 = X^9 + 2X^8 + 2X^3 + X^2$, nous avons $fsep(F) = \{(X, 1, 2), (X+1, 3, 1), (X+2, 1, 4)\}$.

**Exemple 18.4** Lorsque $K = \mathbb{F}_3(T)$ et $F = (X+2T)^7(X^3+2T)^3(X^6+T)$, nous avons $fsep(F) = \{(X+2T, 1, 7), (X+2T^3, 9, 1), (X^2+T, 3, 1)\}$.

## Algorithme de Gianni et Trager

Les algorithmes de factorisation séparable connus utilisent ceux de Musser ou de Yun vus précédemment. Dans le cas de l'algorithme de Musser, le point de départ est la propriété suivante :

**Proposition 18.11** Si l'Algorithme 18.1 est exécuté en caractéristique non nulle, et si $L = \{(G_1, m_1), \dots, (G_s, m_s)\}$ est la liste renvoyée, alors

$$(G_1, 1, m_1), \dots, (G_s, 1, m_s)$$

sont exactement les facteurs séparables de $F$ de la forme $(G, 1, m)$.

Démonstration. Notons $F_0$ le facteur de $F$ composé des racines de multiplicité qui ne sont pas multiples de $p$ et posons $F_1 = F/F_0$. Notons justement $(G_1, 1, m_1), \dots, (G_s, 1, m_s)$ la factorisation séparable de $F_0$. Nous allons montrer que l'algorithme renvoie bien $\{(G_1, m_1), \dots, (G_s, m_s)\}$.