# Algorithmes Efficaces en Calcul Formel

Alin Bostan, Frédéric Chyzak, Marc Giusti, Romain Lebreton, Grégoire Lecerf, Bruno Salvy, Eric Schost

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

## To cite this version:

Alin Bostan, Frédéric Chyzak, Marc Giusti, Romain Lebreton, Grégoire Lecerf, et al. Algorithmes Efficaces en Calcul Formel. published by the Authors, 2017. hal-01431717

HAL Id: hal-01431717

https://inria.hal.science/hal-01431717v1

Submitted on 11 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY-SA 4.0 - Attribution - ShareAlike - International
License
Algorithmes Efficaces
en Calcul Formel

Alin BOSTAN

Frédéric CHYZAK

Marc GIUSTI

Romain LEBRETON

Grégoire LECERF

Bruno SALVY

Éric SCHOST

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

Version préliminaire du 11 janvier 2017
**Alin Bostan**

INRIA,
Université Paris-Saclay

**Frédéric Chyzak**

INRIA,
Université Paris-Saclay

**Marc Giusti**

CNRS, UMR 7161,
École polytechnique,
Laboratoire d'Informatique

**Romain Lebreton**

Université de Montpellier,
Laboratoire d'informatique, de robotique
et de microélectronique de Montpellier,
CNRS, UMR 5506, CC477

**Grégoire Lecerf**

CNRS, UMR 7161,
École polytechnique,
Laboratoire d'Informatique

**Bruno Salvy**

INRIA,
École Normale Supérieure de Lyon,
Laboratoire de l'Informatique du
Parallélisme,
CNRS, UMR 5668

**Éric Schost**

University of Waterloo,
Cheriton School of Computer Science

Copyright © 2016 Alin Bostan, Frédéric Chyzak, Marc Giusti, Romain Lebreton, Grégoire Lecerf, Bruno Salvy, Éric Schost. Cet ouvrage est couvert par une license CC-BY-SA 4.0, https://creativecommons.org/licenses/by-sa/4.0/.

Publié par les auteurs. Site web : https://hal.archives-ouvertes.fr/AECF/.

Le style LATEX utilisé est une adaptation de THE LEGRAND ORANGE BOOK disponible à l'url http://www.LaTeXTemplates.com/.

Version préliminaire du 11 janvier 2017.
# Table des matières

<table><tr><td>Avant-propos</td><td>11</td></tr><tr><td>Bibliographie</td><td>13</td></tr></table>

<table><tr><td>1</td><td>Calcul formel et complexité</td><td>15</td></tr><tr><td>1.1</td><td>Décider, calculer</td><td>15</td></tr><tr><td>1.2</td><td>Calculer rapidement</td><td>22</td></tr><tr><td></td><td>Conventions</td><td>29</td></tr><tr><td></td><td>Notes</td><td>30</td></tr><tr><td></td><td>Bibliographie</td><td>31</td></tr></table>

## I

### Polynômes et séries

<table><tr><td>2</td><td>Multiplication rapide</td><td>35</td></tr><tr><td>2.1</td><td>Introduction, résultats principaux</td><td>35</td></tr><tr><td>2.2</td><td>Algorithme naïf</td><td>38</td></tr><tr><td>2.3</td><td>Algorithme de Karatsuba</td><td>38</td></tr><tr><td>2.4</td><td>Transformée de Fourier rapide</td><td>40</td></tr><tr><td>2.5</td><td>L'algorithme de Schönhage-Strassen</td><td>46</td></tr><tr><td>2.6</td><td>Algorithmes pour les entiers</td><td>49</td></tr><tr><td>2.7</td><td>Un concept important : les fonctions de multiplication</td><td>50</td></tr><tr><td></td><td>Exercices</td><td>51</td></tr><tr><td></td><td>Notes</td><td>52</td></tr><tr><td></td><td>Bibliographie</td><td>54</td></tr></table>
<table><tr><td>3</td><td>Calculs rapides sur les séries</td><td>57</td></tr><tr><td>3.1</td><td>Séries formelles</td><td>58</td></tr><tr><td>3.2</td><td>La méthode de Newton pour le calcul d'inverses</td><td>61</td></tr><tr><td>3.3</td><td>Itération de Newton formelle et applications</td><td>64</td></tr><tr><td>3.4</td><td>La composition des séries</td><td>72</td></tr><tr><td></td><td>Exercices</td><td>75</td></tr><tr><td></td><td>Notes</td><td>76</td></tr><tr><td></td><td>Bibliographie</td><td>78</td></tr></table>

<table><tr><td>4</td><td>Division euclidienne</td><td>81</td></tr><tr><td>4.1</td><td>Introduction</td><td>81</td></tr><tr><td>4.2</td><td>Division de polynômes</td><td>82</td></tr><tr><td>4.3</td><td>Suites récurrentes linéaires à coefficients constants</td><td>84</td></tr><tr><td>4.4</td><td>Développement de fractions rationnelles</td><td>87</td></tr><tr><td>4.5</td><td>Applications</td><td>89</td></tr><tr><td></td><td>Notes</td><td>91</td></tr><tr><td></td><td>Bibliographie</td><td>92</td></tr></table>

<table><tr><td>5</td><td>Calculs modulaires, évaluation et interpolation</td><td>95</td></tr><tr><td>5.1</td><td>Introduction</td><td>95</td></tr><tr><td>5.2</td><td>Présentation, résultats</td><td>96</td></tr><tr><td>5.3</td><td>Interpolation de Lagrange</td><td>98</td></tr><tr><td>5.4</td><td>Algorithmes rapides</td><td>99</td></tr><tr><td></td><td>Exercices</td><td>106</td></tr><tr><td></td><td>Notes</td><td>106</td></tr><tr><td></td><td>Bibliographie</td><td>108</td></tr></table>

<table><tr><td>6</td><td>Pgcd et résultant</td><td>111</td></tr><tr><td>6.1</td><td>Algorithme d'Euclide</td><td>111</td></tr><tr><td>6.2</td><td>Résultant</td><td>116</td></tr><tr><td>6.3</td><td>Algorithme d'Euclide rapide</td><td>130</td></tr><tr><td></td><td>Exercices</td><td>134</td></tr><tr><td></td><td>Notes</td><td>135</td></tr><tr><td></td><td>Bibliographie</td><td>136</td></tr></table>

<table><tr><td>7</td><td>Approximants de Padé et de Padé-Hermite</td><td>139</td></tr><tr><td>7.1</td><td>Reconstruction rationnelle</td><td>139</td></tr><tr><td>7.2</td><td>Approximants de Padé-Hermite</td><td>145</td></tr><tr><td></td><td>Exercices</td><td>152</td></tr><tr><td></td><td>Notes</td><td>153</td></tr><tr><td></td><td>Bibliographie</td><td>154</td></tr></table>
<table><tr><td>8</td><td>Algèbre linéaire dense : de Gauss à Strassen</td><td>159</td></tr><tr><td>8.1</td><td>Introduction</td><td>159</td></tr><tr><td>8.2</td><td>Multiplication de matrices</td><td>162</td></tr><tr><td>8.3</td><td>Autres problèmes d'algèbre linéaire</td><td>170</td></tr><tr><td></td><td>Exercices</td><td>177</td></tr><tr><td></td><td>Notes</td><td>180</td></tr><tr><td></td><td>Bibliographie</td><td>184</td></tr></table>

<table><tr><td>9</td><td>Algèbre linéaire creuse : algorithme de Wiedemann</td><td>189</td></tr><tr><td>9.1</td><td>Introduction</td><td>189</td></tr><tr><td>9.2</td><td>Polynôme minimal et résolution de systèmes</td><td>190</td></tr><tr><td>9.3</td><td>Calcul du polynôme minimal</td><td>190</td></tr><tr><td>9.4</td><td>Calcul du déterminant</td><td>191</td></tr><tr><td>9.5</td><td>Calcul du rang</td><td>192</td></tr><tr><td></td><td>Exercices</td><td>192</td></tr><tr><td></td><td>Notes</td><td>193</td></tr><tr><td></td><td>Bibliographie</td><td>193</td></tr></table>

<table><tr><td>10</td><td>Algèbre linéaire structurée</td><td>195</td></tr><tr><td>10.1</td><td>Introduction</td><td>195</td></tr><tr><td>10.2</td><td>Le cas quasi-Toeplitz</td><td>201</td></tr><tr><td></td><td>Exercices</td><td>205</td></tr><tr><td></td><td>Notes</td><td>206</td></tr><tr><td></td><td>Bibliographie</td><td>206</td></tr></table>

<table><tr><td>11</td><td>Systèmes linéaires à coefficients polynomiaux</td><td>209</td></tr><tr><td>11.1</td><td>Des séries aux solutions rationnelles</td><td>210</td></tr><tr><td>11.2</td><td>Développement comme une fraction rationnelle</td><td>210</td></tr><tr><td>11.3</td><td>L'algorithme de Storjohann</td><td>211</td></tr><tr><td></td><td>Notes</td><td>214</td></tr><tr><td></td><td>Bibliographie</td><td>215</td></tr></table>

<table><tr><td>12</td><td>Principe de transposition de Tellegen</td><td>217</td></tr><tr><td>12.1</td><td>Introduction</td><td>217</td></tr><tr><td>12.2</td><td>La version en termes de graphes du principe de Tellegen</td><td>219</td></tr><tr><td>12.3</td><td>Principe de Tellegen pour les programmes linéaires</td><td>220</td></tr><tr><td>12.4</td><td>Applications</td><td>223</td></tr><tr><td></td><td>Notes</td><td>230</td></tr><tr><td></td><td>Bibliographie</td><td>233</td></tr></table>
<table><tr><td>13</td><td>Équations différentielles à coefficients séries</td><td>237</td></tr><tr><td>13.1</td><td>Équations différentielles d'ordre 1</td><td>237</td></tr><tr><td>13.2</td><td>Équations d'ordre supérieur et systèmes d'ordre 1</td><td>240</td></tr><tr><td>13.3</td><td>Cas particuliers</td><td>244</td></tr><tr><td>13.4</td><td>Extensions</td><td>246</td></tr><tr><td></td><td>Notes</td><td>247</td></tr><tr><td></td><td>Bibliographie</td><td>247</td></tr></table>

III

Équations différentielles et récurrences linéaires

<table><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>252</td></tr><tr><td>14.2</td><td>Propriétés de clôture</td><td>257</td></tr><tr><td>14.3</td><td>Séries algébriques</td><td>260</td></tr><tr><td>14.4</td><td>Au-delà</td><td>264</td></tr><tr><td></td><td>Exercices</td><td>264</td></tr><tr><td></td><td>Notes</td><td>266</td></tr><tr><td></td><td>Bibliographie</td><td>267</td></tr></table>

<table><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>270</td></tr><tr><td>15.2</td><td>Pas de bébés, pas de géants</td><td>272</td></tr><tr><td>15.3</td><td>Scindage binaire</td><td>275</td></tr><tr><td>15.4</td><td>Récurrences singulières</td><td>278</td></tr><tr><td></td><td>Exercices</td><td>285</td></tr><tr><td></td><td>Notes</td><td>286</td></tr><tr><td></td><td>Bibliographie</td><td>287</td></tr></table>

<table><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>290</td></tr><tr><td>16.2</td><td>Solutions à support fini</td><td>295</td></tr><tr><td>16.3</td><td>Solutions polynomiales</td><td>297</td></tr><tr><td>16.4</td><td>Solutions rationnelles</td><td>303</td></tr><tr><td></td><td>Exercices</td><td>309</td></tr><tr><td></td><td>Notes</td><td>310</td></tr><tr><td></td><td>Bibliographie</td><td>310</td></tr></table>

<table><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>312</td></tr><tr><td>17.2</td><td>Solutions polynomiales</td><td>315</td></tr><tr><td>17.3</td><td>Solutions rationnelles</td><td>315</td></tr></table>
<table><tr><td>17.4</td><td>Séries différentiellement finies analytiques</td><td>317</td></tr><tr><td></td><td>Exercice</td><td>320</td></tr><tr><td></td><td>Notes</td><td>321</td></tr><tr><td></td><td>Bibliographie</td><td>322</td></tr></table>

IV

Factorisation des polynômes

<table><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>327</td></tr><tr><td>18.2</td><td>Factorisation sans carré</td><td>328</td></tr><tr><td>18.3</td><td>Séparabilité et déflation</td><td>331</td></tr><tr><td>18.4</td><td>Factorisation séparable</td><td>334</td></tr><tr><td>18.5</td><td>Réduction à la factorisation des polynômes séparables</td><td>341</td></tr><tr><td></td><td>Notes</td><td>342</td></tr><tr><td></td><td>Bibliographie</td><td>342</td></tr></table>

<table><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>345</td></tr><tr><td>19.2</td><td>Algorithme de Berlekamp</td><td>348</td></tr><tr><td>19.3</td><td>Algorithme de Cantor et Zassenhaus</td><td>353</td></tr><tr><td></td><td>Notes</td><td>356</td></tr><tr><td></td><td>Bibliographie</td><td>357</td></tr></table>

<table><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>359</td></tr><tr><td>20.2</td><td>Le procédé d'orthogonalisation de Gram-Schmidt</td><td>362</td></tr><tr><td>20.3</td><td>L'algorithme LLL</td><td>364</td></tr><tr><td>20.4</td><td>Preuve de l'algorithme LLL</td><td>366</td></tr><tr><td></td><td>Notes</td><td>370</td></tr><tr><td></td><td>Bibliographie</td><td>371</td></tr></table>

<table><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>374</td></tr><tr><td>21.2</td><td>Factorisation sans carré</td><td>379</td></tr><tr><td>21.3</td><td>Algorithme par remontée et recombinaison</td><td>385</td></tr><tr><td>21.4</td><td>Algorithme de van Hoeij</td><td>390</td></tr><tr><td></td><td>Notes</td><td>396</td></tr><tr><td></td><td>Bibliographie</td><td>397</td></tr></table>

<table><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>400</td></tr><tr><td>22.2</td><td>Recombinaison des dérivées logarithmiques</td><td>403</td></tr><tr><td>22.3</td><td>Algorithme de Lecerf</td><td>406</td></tr></table>
<table><tr><td>22.4</td><td>Réduction probabiliste au cas de deux variables</td><td>409</td></tr><tr><td>22.5</td><td>Cas particulier des corps finis</td><td>411</td></tr><tr><td></td><td>Notes</td><td>412</td></tr><tr><td></td><td>Bibliographie</td><td>413</td></tr></table>

V

Systèmes polynomiaux

<table><tr><td>23</td><td>Bases standard</td><td>419</td></tr><tr><td>23.1</td><td>Lien entre algèbre et géométrie</td><td>421</td></tr><tr><td>23.2</td><td>Ordres totaux admissibles sur le monoïde des monômes</td><td>425</td></tr><tr><td>23.3</td><td>Exposants privilégiés et escaliers</td><td>428</td></tr><tr><td>23.4</td><td>Noethérianité et bases standard</td><td>429</td></tr><tr><td>23.5</td><td>Divisions</td><td>431</td></tr><tr><td></td><td>Notes</td><td>434</td></tr><tr><td></td><td>Bibliographie</td><td>435</td></tr></table>

<table><tr><td>24</td><td>Construction de bases standard</td><td>437</td></tr><tr><td>24.1</td><td>Algorithme naïf par algèbre linéaire</td><td>437</td></tr><tr><td>24.2</td><td>Polynôme de syzygie</td><td>438</td></tr><tr><td>24.3</td><td>L'algorithme de construction de Buchberger</td><td>439</td></tr><tr><td>24.4</td><td>Exemple de calcul</td><td>441</td></tr><tr><td>24.5</td><td>Propriétés des bases standard pour quelques ordres</td><td>442</td></tr><tr><td></td><td>Notes</td><td>444</td></tr><tr><td></td><td>Bibliographie</td><td>445</td></tr></table>

<table><tr><td>25</td><td>Nullstellensatz et applications</td><td>447</td></tr><tr><td>25.1</td><td>Nullstellensatz affine</td><td>447</td></tr><tr><td>25.2</td><td>Idéaux de dimension zéro</td><td>450</td></tr><tr><td>25.3</td><td>Projection des variétés affines</td><td>459</td></tr><tr><td>25.4</td><td>Nullstellensatz projectif</td><td>461</td></tr><tr><td>25.5</td><td>Projection des variétés projectives</td><td>462</td></tr><tr><td></td><td>Notes</td><td>464</td></tr><tr><td></td><td>Bibliographie</td><td>464</td></tr></table>

<table><tr><td>26</td><td>Fonction et polynôme de Hilbert, dimension, degré</td><td>467</td></tr><tr><td>26.1</td><td>Fonction et polynôme de Hilbert</td><td>467</td></tr><tr><td>26.2</td><td>Démonstrations et borne supérieure de la régularité</td><td>469</td></tr><tr><td>26.3</td><td>Dimension et géométrie</td><td>470</td></tr><tr><td>26.4</td><td>Position de Noether et degré</td><td>472</td></tr><tr><td>26.5</td><td>Suites régulières</td><td>475</td></tr><tr><td>26.6</td><td>Bases standard pour les suites régulières en position de Noether</td><td>478</td></tr><tr><td></td><td>Notes</td><td>480</td></tr></table>
<table><tr><td></td><td>Bibliographie</td><td>481</td></tr><tr><td>27</td><td>Normalisation de Noether</td><td>483</td></tr><tr><td>27.1</td><td>Introduction</td><td>483</td></tr><tr><td>27.2</td><td>Algorithme de mise en position de Noether</td><td>484</td></tr><tr><td>27.3</td><td>Utilisations de la position de Noether</td><td>489</td></tr><tr><td>27.4</td><td>Paramétrisation en dimension quelconque</td><td>496</td></tr><tr><td></td><td>Notes</td><td>500</td></tr><tr><td></td><td>Bibliographie</td><td>500</td></tr><tr><td>28</td><td>Résolution géométrique</td><td>501</td></tr><tr><td>28.1</td><td>Introduction</td><td>501</td></tr><tr><td>28.2</td><td>Approche incrémentale</td><td>507</td></tr><tr><td>28.3</td><td>Remontée d'une paramétrisation</td><td>515</td></tr><tr><td>28.4</td><td>Algorithme principal</td><td>519</td></tr><tr><td></td><td>Notes</td><td>522</td></tr><tr><td></td><td>Bibliographie</td><td>523</td></tr></table>

VI

Sommation et intégration de suites et fonctions spéciales

<table><tr><td>29</td><td>Sommation hypergéométrique, solutions hypergéométriques</td><td>527</td></tr><tr><td>29.1</td><td>Sommation hypergéométrique indéfinie. Algorithme de Gosper</td><td>528</td></tr><tr><td>29.2</td><td>Sommation hypergéométrique définie. Algorithme de Zeilberger</td><td>534</td></tr><tr><td>29.3</td><td>Solutions hypergéométriques. Algorithme de Petkovšek</td><td>539</td></tr><tr><td></td><td>Notes</td><td>541</td></tr><tr><td></td><td>Bibliographie</td><td>542</td></tr><tr><td>30</td><td>Équations fonctionnelles linéaires et polynômes tordus</td><td>543</td></tr><tr><td>30.1</td><td>Polynômes non commutatifs comme opérateurs linéaires</td><td>543</td></tr><tr><td>30.2</td><td>Morphismes entre anneaux de polynômes tordus</td><td>547</td></tr><tr><td>30.3</td><td>Division euclidienne</td><td>549</td></tr><tr><td>30.4</td><td>Recherche de solutions et factorisation d'opérateurs</td><td>551</td></tr><tr><td>30.5</td><td>Algorithme d'Euclide</td><td>553</td></tr><tr><td>30.6</td><td>Relations de contiguïté</td><td>555</td></tr><tr><td></td><td>Notes</td><td>558</td></tr><tr><td></td><td>Bibliographie</td><td>559</td></tr><tr><td>31</td><td>Algèbres de Ore</td><td>561</td></tr><tr><td>31.1</td><td>Algèbres de Ore rationnelles</td><td>561</td></tr><tr><td>31.2</td><td>Idéal annulateur et module quotient</td><td>562</td></tr><tr><td>31.3</td><td>Bases de Gröbner pour les idéaux à gauche</td><td>564</td></tr><tr><td>31.4</td><td>Module quotient et dimension de l'espace des solutions</td><td>566</td></tr></table>
<table><tr><td>31.5</td><td>Les fonctions ∂-finies et leurs clôtures</td><td>570</td></tr><tr><td></td><td>Notes</td><td>576</td></tr><tr><td></td><td>Bibliographie</td><td>576</td></tr></table>

<table><tr><td>32</td><td>Sommation et intégration symboliques des fonctions spéciales</td><td>579</td></tr><tr><td>32.1</td><td>Expression du télescopage créatif</td><td>580</td></tr><tr><td>32.2</td><td>L'algorithme de sommation ∂-finie définie sur un exemple</td><td>581</td></tr><tr><td>32.3</td><td>Description des algorithmes de sommation et intégration ∂-finies</td><td>584</td></tr><tr><td>32.4</td><td>Bases de Gröbner de modules et découplage de systèmes</td><td>586</td></tr><tr><td></td><td>Notes</td><td>588</td></tr><tr><td></td><td>Bibliographie</td><td>588</td></tr></table>

Annexes

<table><tr><td>33</td><td>Exercices récapitulatifs</td><td>591</td></tr><tr><td></td><td>Notes</td><td>619</td></tr><tr><td></td><td>Bibliographie</td><td>619</td></tr><tr><td></td><td>Bibliographie générale</td><td>621</td></tr><tr><td></td><td>Liste des notations</td><td>653</td></tr><tr><td></td><td>Liste des figures et algorithmes</td><td>655</td></tr><tr><td></td><td>Index des auteurs</td><td>659</td></tr><tr><td></td><td>Index</td><td>667</td></tr></table>
# Avant-propos

Le calcul formel traite des objets mathématiques exacts d'un point de vue informatique. Cet ouvrage « Algorithmes efficaces en calcul formel » explore deux directions : la calculabilité et la complexité. La calculabilité étudie les classes d'objets mathématiques sur lesquelles des réponses peuvent être obtenues algorithmiquement. La complexité donne ensuite des outils pour comparer des algorithmes du point de vue de leur efficacité.

Les trois premières parties passent en revue l'algorithme efficace sur les objets fondamentaux que sont les entiers, les polynômes, les matrices, les séries et les solutions d'équations différentielles ou de récurrences linéaires. Nous y montrons que de nombreuses questions portant sur ces objets admettent une réponse en complexité optimale ou quasi-optimale, en insistant sur les principes généraux de conception d'algorithmes efficaces.

La quatrième partie est dédiée à la factorisation des polynômes. Il s'agit d'un des premiers problèmes d'une importance majeure en calcul formel, pour lequel il existe encore plusieurs questions ouvertes conséquentes. Nous y présentons les algorithmes classiques de factorisation séparable, sans carré et irréductible. Nous traitons tout d'abord le cas des polynômes à une variable dont les coefficients sont dans un corps fini puis des nombres rationnels. Le cas des polynômes à plusieurs variables vient ensuite et nous nous concentrons essentiellement sur le cas de deux variables pour les résultats de complexité.

La cinquième partie aborde les systèmes d'équations polynomiales. Il s'agit de montrer comment répondre à des questions de nature géométrique posées sur les solutions de ces systèmes, tant par des approches à base de récriture (bases standard ou de Gröbner) que plus directement liées à la structure de l'ensemble décrit (résolution géométrique).

La sixième et dernière partie traite de l'intégration et de la sommation définie. De nombreux calculs d'intégrales et de sommes de fonctions ou de suites spéciales de
la combinatoire ou de la physique mathématique peuvent être abordés algorithmiquement. Pour résoudre ces problèmes efficacement, presque tous les algorithmes abordés au cours des cinq premières parties sont utilisés.

Nous n'avons pas cherché à écrire un document de référence pour le calcul formel. Au contraire, nous avons souvent choisi de mettre l'accent sur les idées clés, quitte à présenter sous forme d'exercices les aspects les plus techniques, et à renvoyer le lecteur à des ouvrages ou articles de recherche qui contiennent les solutions. Nous détaillons un grand nombre de résultats très récents qui conduisent souvent le lecteur à la frontière des connaissances. Ainsi, une portion importante des parties III à VI n'était à notre connaissance jamais apparue dans un livre.

Il y a aussi des thèmes classiques du calcul formel que nous n'abordons pas, comme la recherche de primitives (algorithme de Risch) ou de solutions liouvilliennes d'équations différentielles. Nous préférons en effet privilégier un point de vue général où les équations sont vues comme des opérateurs servant à représenter des solutions.

Les prérequis mathématiques nécessaires à la lecture sont modestes et de nombreuses définitions classiques sont rappelées et illustrées au fur et à mesure des besoins. Le texte est par conséquent accessible tant à des étudiants en mathématiques qu'en informatique.

Nous espérons pouvoir être utiles à tous les étudiants et jeunes ou moins jeunes chercheurs et enseignants en mathématiques effectives, ou plus spécifiquement en calcul formel, désireux d'approfondir leurs connaissances ou à la recherche de matériaux pédagogiques.

Cet ouvrage est une synthèse de notes de cours rédigées principalement pour le cours du même nom que nous avons donné pendant plus de dix ans au *Master Pari-sien de Recherche en Informatique* de l'Université Paris Diderot, des Écoles Normales Supérieures de Cachan et de Paris, et de l'École polytechnique. La partie concernant les systèmes polynomiaux provient aussi de notes de cours donnés au DEA *Méthodes algébriques* puis au *Master Algèbre et Géométrie* de l'Université Pierre-et-Marie-Curie, mais aussi au DEA *Informatique Mathématique et Applications* de l'École Normale Supérieure de Paris, l'École polytechnique, l'Université Pierre-et-Marie-Curie, l'Université Paris Diderot, et l'Université Paris-Sud. Plusieurs parties ont aussi fait l'objet de mini-cours plus spécialisés donnés à l'occasion des *Journées Nationales de Calcul Formel* en 2007, 2010, 2011, et 2013.

Des versions préliminaires de parties de ce livre ont bénéficié des commentaires de relecteurs que nous remercions ici : Philippe Dumas, Claude-Pierre Jeannerod, Marc Mezzarobba, Pierre Nicodème, Anne Vaugon, Gilles Villard.

## Ouvrages de référence

Pour les prérequis en informatique, sur les modèles de calcul et la théorie de la complexité, ainsi que sur les algorithmes de base, nous renvoyons le lecteur aux livres de Aho, Hopcroft, Ullman [AHU74], de Knuth [Knu97], de Cormen, Leiserson, Rivest et Stein [Cor+09], et de Stern [Ste94]. Quant aux prérequis en mathématiques, ils concernent essentiellement des points classiques d'algèbre, couverts par exemple par le livre de Lang [Lan02].
Les références générales sur les algorithmes du calcul formel sont deux livres : celui de von zur Gathen et Gerhard [GG03] et celui, plus élémentaire, de Geddes, Czapor et Labahn [GCL92]. La complexité est également utilisée comme fil conducteur dans le livre de Bini et Pan [BP94a] ainsi que dans celui, plus difficile, de Bürgisser, Clausen et Shokrollahi [BCS97]. Mentionnons que les ouvrages, de Geddes, Czapor, et Labahn [GCL92], de Zippel [Zip93], et de Yap [Yap00], couvrent aussi des sujets absents ici comme les manipulations d'expressions symboliques, l'arithmétique des polynômes à plusieurs variables, l'intégration symbolique, ou encore le calcul des racines des polynômes à une variable à coefficients rationnels.

Les aspects calculatoires en algèbre linéaire sont très vastes et il nous a fallu
nous restreindre, en deuxième partie, à couvrir seulement les quelques algorithmes
nécessaires aux parties suivantes. Les livres de Pan [Pan01; Pan84b], et d'Abdeljaoued
et Lombardi [AL04] peuvent être utilisés comme références sur ce sujet.

Une bonne introduction à notre cinquième partie, sur les systèmes polynomiaux,
est le livre de Cox, Little et O'Shea [CLO96]. Parmi les ouvrages sur les bases de
Gröbner, le livre de Becker et Weispfenning [BW93] reste une référence très complète.
Un point de vue plus mathématique est développé dans le livre d'Eisenbud [Eis95].

L'algorithmique des corps finis est assez peu présente ici. Ce sujet intervient
désormais davantage dans les cours de cryptographie et codes correcteurs. Les livres
récents de Shoup [Sho09] et Joux [Jou09] fournissent des contenus très aboutis dans
cette direction.

De premiers éléments pour aborder notre troisième partie, sur les équations
différentielles et les récurrences linéaires, sont donnés dans le livre de Petkovšek,
Wilf et Zeilberger [PWZ96]. Enfin, le livre de Saito, Sturmfels et Takayama [SST00]
couvre, en différentiel, certains aspects de la théorie de l'élimination dans le cadre
non commutatif de notre sixième partie.

Bibliographie

AHU74 AHO, Alfred V., John E. HOPCROFT et Jeffrey D. ULLMAN (1974). *The design and analysis of computer algorithms*. Addison-Wesley Publishing Co.

AL04 ABDELJAOUED, J. et H. LOMBARDI (2004). *Méthodes matricielles : introduction à la complexité algébrique*. Vol. 42. Mathématiques & Applications. Springer-Verlag.

BCS97 BÜRGISSER, Peter, Michael CLAUSEN et M. Amin SHOKROLLAHI (1997). *Algebraic complexity theory*. Vol. 315. Grundlehren der Mathematischen Wissenschaften. Springer-Verlag.

BP94a BINI, Dario et Victor Y. PAN (1994). *Polynomial and matrix computations. Vol. 1*. Progress in Theoretical Computer Science. Fundamental algorithms. Birkhäuser Boston, Inc., Boston, MA,

BW93 BECKER, Thomas et Volker WEISPFENNING (1993). *Gröbner bases. A computational approach to commutative algebra*. Vol. 141. Graduate Texts in Mathematics. Springer-Verlag New York.

CLO96 Cox, David, John LITTLE et Donal O'SHEA (1996). *Ideals, varieties, and algorithms. An introduction to computational algebraic geometry and commutative algebra*. 2e éd. Undergraduate Texts in Mathematics. Springer-Verlag.
Cor+09 CORMEN, Thomas H., Charles E. LEISERSON, Ronald L. RIVEST et Clifford STEIN (2009). *Introduction to algorithms*. Third. MIT Press, Cambridge, MA.

Eis95 EISENBUD, David (1995). *Commutative algebra. With a view toward algebraic geometry*. Vol. 150. Graduate Texts in Mathematics. Springer-Verlag New York.

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

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

Jou09 JOUX, Antoine (2009). *Algorithmic cryptanalysis*. Chapman & Hall/CRC Cryptography and Network Security Series. Chapman et Hall/CRC.

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

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

Pan01 PAN, Victor Y. (2001). *Structured matrices and polynomials*. Unified superfast algorithms. Birkhäuser Boston Inc.

Pan84b PAN, Victor (1984). *How to multiply matrices faster*. Vol. 179. Lecture Notes in Computer Science. Springer-Verlag.

PWZ96 PETKOVŠEK, Marko, Herbert S. WILF et Doron ZEILBERGER (1996). A = B. A. K. Peters.

Sho09 SHOUP, Victor (2009). *A computational introduction to number theory and algebra*. 2e éd. Cambridge University Press.

SST00 SAITO, Mutsumi, Bernd STURMFELS et Nobuki TAKAYAMA (2000). *Gröbner deformations of hypergeometric differential equations*. Springer-Verlag.

Ste94 STERN, Jacques (1994). *Fondements mathématiques de l'informatique*. Éd science international.

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

Zip93 ZIPPEL, Richard (1993). *Effective polynomial computation*. Kluwer Academic Publishers.
# 1. Calcul formel et complexité

## Résumé

Le calcul formel étudie dans quelle mesure un ordinateur peut « faire des mathématiques », et pour celles qu'il peut faire, à quelle vitesse. Ce chapitre introductif présente d'abord les problèmes de calculabilité. Le point de départ est le théorème plutôt négatif de Richardson-Matiyasevich. Ensuite, un survol rapide des thèmes qui seront abordés tout au long de cet ouvrage montre que malgré ce théorème, bien des questions peuvent être traitées algorithmiquement. La seconde partie du chapitre aborde la théorie de la complexité, pose les conventions et les définitions principales, ainsi que les résultats de base sur la méthodologie « diviser pour régner ».

## 1.1 Décider, calculer

### Fondements logiques

D'une certaine manière, le calcul formel est fondé sur une contrainte d'origine logique.

**Théorème 1.1** — Richardson-Matiyasevich. Dans la classe des expressions formées à partir d'une variable X et de la constante 1 par les opérations d'anneau $+, -, \times$ et la composition avec la fonction $sin(\cdot)$ et la fonction valeur absolue $|\cdot|$, le test d'équivalence à 0 est indécidable.

Autrement dit, il n'existe pas d'algorithme permettant pour toute expression de cette classe de déterminer en temps fini si elle vaut 0 ou non. Plus généralement tout test d'égalité peut bien entendu se ramener à tester l'égalité à zéro dès que la soustraction existe. Cette limitation de nature théorique explique la difficulté et
parfois la frustration que rencontrent les utilisateurs débutants des systèmes de calcul formel face à des fonctions de « simplification », qui ne peuvent être qu'heuristiques.

Pour effectuer un calcul, il est pourtant souvent crucial de déterminer si des expressions représentent 0 ou non, en particulier pour évaluer une fonction qui possède des singularités (comme la division). L'approche du calculateur formel expérimenté consiste à se ramener autant que faire se peut à des opérations d'un domaine dans lequel le test à zéro est décidable. Le calcul formel repose ainsi de manière naturelle sur des constructions algébriques qui préservent la décidabilité du test à 0. En particulier, les opérations courantes sur les vecteurs, matrices, polynômes, fractions rationnelles, ne nécessitent pas d'autre test à 0 que celui des coefficients. La notion d'effectivité permet de préciser ce point de vue.

**Définition 1.1** Une structure algébrique (groupe, anneau, corps, espace vectoriel, ...) est dite *effective* si l'on dispose :

- d'une structure de données pour en représenter les éléments;

- d'algorithmes pour en effectuer les opérations et pour y tester l'égalité et autres prédicats.

Par exemple, dans un anneau effectif, outre l'égalité, les opérations requises sont l'addition, la soustraction et la multiplication. Pour une structure ordonnée, des algorithmes sont aussi requis pour effectuer les comparaisons. Cette notion d'effectivité est bien évidemment attachée à un modèle de calcul. Dans cet ouvrage nous considérons principalement les machines à accès direct (voir ci-dessous). Néanmoins il est connu que les fonctions calculables par ces machines sont les mêmes que pour les machines de Turing (mais avec des temps d'exécution bien sûr différents).

## Structures et constructions de base

Les objets les plus fondamentaux sont assez faciles à représenter en machine de manière exacte. Nous considérons tour à tour les plus importants d'entre eux, en commençant par les plus élémentaires. Ils s'assemblent ensuite à l'aide de tableaux ou de listes pour en former de plus complexes.

### Entiers machine

Les entiers fournis par les processeurs sont des entiers modulo une puissance de 2 (le nombre de bits d'un mot machine, typiquement 32 ou 64). Ils sont appelés des *entiers machine*. Les opérations rendues disponibles par le processeur sont l'addition, la soustraction, la multiplication et parfois la division. La norme ANSI du langage C fournit au programmeur la division et le modulo pour ces entiers, c'est-à-dire que le compilateur implante ces opérations si le processeur ne le fait pas.

### Entiers

Pour manipuler des entiers dont la taille dépasse celle d'un mot machine, il est commode de les considérer comme écrits dans une base B assez grande :

$$N = a_0 + a_1 B + \dots + a_k B^k.$$

L'écriture est unique si l'on impose $0 \le a_i < B$. (Le signe est stocké séparément.) Ces nombres peuvent être stockés dans des tableaux d'entiers machine. Les objets obtenus
sont des entiers de taille arbitraire appelés parfois *bignums*.

L'addition et le produit peuvent alors être réduits à des opérations sur des entiers inférieurs à $B^2$, au prix de quelques opérations de propagation de retenue. Le choix de B dépend un peu du processeur. Si le processeur dispose d'une instruction effectuant le produit de deux entiers de taille égale à celle d'un mot machine, renvoyant le résultat dans deux mots machines, alors B pourra être pris aussi grand que le plus grand entier tenant dans un mot machine. Sinon, c'est la racine carrée de ce nombre qui sera utilisée pour B. Quoiqu'il en soit, nous avons donc le résultat suivant :

**Lemme 1.2** L'anneau $\mathbb{Z}$ des entiers relatifs est effectif.

## Entiers modulaires

Les calculs avec des polynômes, des fractions rationnelles ou des matrices à coefficients entiers souffrent souvent d'une maladie propre au calcul formel : la croissance des expressions intermédiaires. Les entiers produits comme coefficients des expressions intervenant lors du calcul sont de taille disproportionnée par rapport à ceux qui figurent dans l'entrée et dans la sortie.

**Exemple 1.1** Voici le déroulement typique du calcul du plus grand diviseur commun (pgcd) de deux polynômes à coefficients entiers par l'algorithme d'Euclide (Chapitre 6), où nous notons $\text{rem}(\cdot, \cdot)$ le reste de la division :

$$P_0 = 7X^5 - 22X^4 + 55X^3 + 94X^2 - 87X + 56,$$

$$P_1 = 62X^4 - 97X^3 + 73X^2 + 4X + 83,$$

$$P_2 = \text{rem}(P_0, P_1) = \frac{113293}{3844}X^3 + \frac{409605}{3844}X^2 - \frac{183855}{1922}X + \frac{272119}{3844},$$

$$P_3 = \text{rem}(P_0, P_1) = \frac{18423282923092}{12835303849}X^2 - \frac{15239170790368}{12835303849}X + \frac{10966361258256}{12835303849},$$

$$P_4 = \text{rem}(P_0, P_1) = -\frac{216132274653792395448637}{44148979404824831944178}X - \frac{631179956389122192280133}{88297958809649663888356},$$

$$P_5 = \text{rem}(P_0, P_1) = \frac{20556791167692068695002336923491296504125}{3639427682941980248860941972667354081}.$$

Les coefficients de ces polynômes intermédiaires font intervenir des entiers qui croissent de manière exponentielle, alors que le résultat recherché est 1.

Grâce à la division euclidienne, on obtient immédiatement.

**Lemme 1.3** Pour tout entier $n \ge 2$, l'anneau $\mathbb{Z}/n\mathbb{Z}$ des entiers modulo $n$ est effectif.

Ces entiers modulaires remédient à ce problème de croissance intermédiaire de deux manières. D'une part, pour un calcul de décision, de dimension, ou de degré, l'exécution de l'algorithme sur la réduction de l'entrée modulo un nombre premier donne un algorithme *probabiliste* répondant à la question. Cette technique peut aussi servir de base à un algorithme *déterministe* lorsque les nombres premiers pour lesquels la réponse est fausse peuvent être maîtrisés. Nous verrons comment exploiter cette
idée pour le calcul du pgcd dans le Chapitre 6.

D'autre part, les entiers modulaires sont utilisés dans les algorithmes reposant sur
le théorème des restes chinois. Ce théorème indique qu'un entier inférieur au produit
de nombres premiers $p_1 \cdots p_k$ peut être reconstruit à partir de ses réductions modulo
$p_1, \dots, p_k$. Lorsqu'une borne sur la taille du résultat est disponible, il suffit d'effectuer
le calcul modulo suffisamment de nombres premiers (choisis en pratique assez grands
pour que leur nombre soit faible et assez petits pour que les opérations tiennent dans
un mot machine), pour ensuite reconstruire le résultat, court-circuitant de la sorte
toute croissance intermédiaire.

## Vecteurs et matrices

Une fois donnée une représentation exacte pour des coefficients, il est facile de
construire des vecteurs ou matrices comme des tableaux, ou plus souvent comme des
tableaux de pointeurs sur les coefficients. Les opérations de produit par un scalaire,
de produit de matrices ou de produit d'une matrice par un vecteur se réduisent aux
opérations d'addition et de multiplication sur les coefficients. Il en va de même pour
la recherche de noyau ou d'inverse de matrices.

**Proposition 1.4** Si $K$ est un corps effectif, l'espace vectoriel $K^n$ l'est aussi, ainsi que l'anneau $M_n(K)$.

## Polynômes

Les polynômes peuvent être représentés de plusieurs manières, et la meilleure
représentation dépend des opérations que l'on souhaite effectuer. Pour un polynôme
en une variable, les choix principaux sont :

- la représentation dense : comme pour les entiers, le polynôme est représenté comme un tableau de (pointeurs sur les) coefficients;

- la représentation creuse : le polynôme est représenté comme une liste de paires (coefficient, exposant) généralement triée par les exposants.

Dans les deux cas, nous avons clairement la propagation de l'effectivité aux anneaux
de polynômes :

**Proposition 1.5** Si $A$ est un anneau effectif, alors $A[X]$ l'est aussi.

L'usage itéré de cette proposition fournit les polynômes à plusieurs variables.

## Fractions rationnelles

Les rationnels peuvent être stockés comme des paires où numérateur et dénomina-
teur sont des entiers de taille arbitraire. Les opérations d'addition et de multiplication
se réduisent aux opérations analogues sur les entiers et le test d'égalité à zéro se
réduit au test d'égalité à 0 sur le numérateur. De même, les fractions rationnelles sont
représentées par des paires de polynômes. Les opérations d'addition, produit, division,
se réduisent aux additions et multiplications sur les coefficients. Plus généralement,
nous obtenons le résultat suivant pour les corps de fractions rationnelles :
**Proposition 1.6** Si $\mathbb{A}$ est un anneau intègre effectif, alors son corps des fractions est effectif.

## Séries tronquées

Les séries tronquées

$$ \sum_{k=0}^{N} a_k X^k + O(X^{N+1}) $$

se représentent pratiquement comme des polynômes. La différence principale ap-
paraît lors du produit : les coefficients des termes d'exposant au moins N + 1 n'ont
pas besoin d'être calculés, ni stockés. La structure considérée alors est l'anneau
quotient $\mathbb{A}[X]/(X^{N+1})$ obtenu en faisant le quotient de l'anneau $\mathbb{A}[X]$ par l'idéal
qu'y engendre $X^{N+1}$. Les éléments de ce quotient sont représentés par le $(N+1)$-
uplet de coefficients $(a_0, \dots, a_N)$. Le test à 0 revient à tester que ces coefficients sont
tous nuls.

**Proposition 1.7** Si $\mathbb{A}$ est un anneau effectif et $N \in \mathbb{N}$, alors $\mathbb{A}[X]/(X^{N+1})$ est un anneau effectif.

Cette structure de données joue un rôle très important non seulement pour des
calculs d'approximations, mais aussi comme une représentation *exacte*. En voici trois
exemples :

1. Une fraction rationnelle dont les numérateurs et dénominateurs ont un degré borné par $d$ peut être reconstruite à partir d'un développement en série à l'ordre $2d+1$. Cette représentation joue ainsi un rôle clé dans la manipulation des nombres algébriques (Chapitre 3), et dans le calcul efficace de la division euclidienne de polynômes, de suites récurrentes linéaires (comme le calcul rapide du 10 000ᵉ nombre de Fibonacci au Chapitre 4) et du polynôme minimal d'une matrice creuse (Chapitre 9).

2. Il est possible de reconstruire une équation différentielle linéaire à coefficients polynomiaux à partir du développement en série d'une solution et de bornes sur l'ordre et le degré des coefficients. De façon analogue, il est possible de reconstruire une récurrence linéaire à coefficients polynomiaux à partir des premières valeurs d'une de ses solutions. L'outil algorithmique pour effectuer ces calculs de *devinette* (*guessing* en anglais) est le calcul rapide d'approximants de Padé-Hermite (Chapitre 7).

3. Un polynôme en deux variables peut être reconstruit à partir du développement en série d'une solution. L'efficacité de la factorisation des polynômes à deux variables, ou encore de la résolution de systèmes polynomiaux par la méthode dite de *résolution géométrique*, abordée au Chapitre 28 repose de manière cruciale sur cette opération, qui doit être effectuée rapidement.

## Équations comme structures de données

Une fois construits les objets de base que sont les polynômes, les séries ou les matrices, il est possible d'aborder des objets mathématiques construits *implicitly*. Ainsi, il est bien connu qu'il n'est pas possible de représenter toutes les solutions
de polynômes de haut degré par radicaux, mais de nombreuses opérations sur ces solutions sont aisées en prenant le polynôme lui-même comme structure de données. Ce point de vue permet d'étendre le domaine d'application du calcul formel pourvu que des algorithmes soient disponibles pour effectuer les opérations souhaitées (typiquement addition, multiplication, multiplication par un scalaire, test d'égalité) par manipulation des équations elles-mêmes.

## Nombres algébriques

C'est ainsi que l'on nomme les solutions de polynômes à une variable. Le résultat est spectaculaire.

**Proposition 1.8** Si $\mathbb{K}$ est un corps effectif, alors sa clôture algébrique $\overline{\mathbb{K}}$ l'est aussi.

Les opérations d'addition et de multiplication peuvent être effectuées à l'aide de résultats (Chapitre 6). Ceux-ci peuvent être calculés efficacement à l'aide de séries (Chapitre 3). La division s'obtient par l'algorithme d'Euclide sur les polynômes (Chapitre 6), et le test d'égalité se déduit du pgcd. Par exemple, il est possible de prouver assez facilement une identité comme

$$ \frac{\sin \frac{2\pi}{7}}{\sin^2 \frac{3\pi}{7}} - \frac{\sin \frac{\pi}{7}}{\sin^2 \frac{2\pi}{7}} + \frac{\sin \frac{3\pi}{7}}{\sin^2 \frac{\pi}{7}} = 2\sqrt{7} \quad (1.1) $$

une fois que l'on reconnaît qu'il s'agit d'une égalité entre nombres algébriques.

## Systèmes polynomiaux

Vu l'importance des systèmes polynomiaux en calcul formel, une grande partie de cet ouvrage leur sera consacrée. Un résultat majeur de cette partie est le suivant.

**Proposition 1.9** Soient $\mathbb{K}$ un corps effectif, et $f_1, \dots, f_p$ des polynômes de l'anneau $\mathbb{K}[X_1, \dots, X_n]$. Alors l'anneau quotient $\mathbb{K}[X_1, \dots, X_n]/(f_1, \dots, f_n)$ est effectif.

Ce quotient est un outil de base pour répondre à de nombreuses questions naturelles sur un système de polynômes, comme l'existence de solutions, la dimension de l'espace des solutions (qui indique s'il s'agit d'une surface, d'une courbe, ou de points isolés), le degré, ou le calcul d'une paramétrisation de l'ensemble des solutions.

Il est également possible d'éliminer une ou des variables entre des polynômes. Cette opération s'interprète géométriquement comme une projection (Chapitre 25). Dans le cas le plus simple, elle permet de calculer un polynôme s'annulant sur les abscisses des intersections de deux courbes (Chapitre 6). Une autre application est le calcul d'une équation d'une courbe donnée sous forme paramétrée.

## Équations différentielles linéaires

Cette structure de données permet de représenter de nombreuses fonctions usuelles (exponentielle, fonctions trigonométriques et trigonométriques hyperboliques, leurs réciproques), ainsi que de nombreuses fonctions spéciales de la physique mathématique (fonctions de Bessel, de Struve, d'Anger, ..., fonctions hypergéométriques
et hypergéométriques généralisées), ainsi bien sûr que de multiples fonctions auxquelles n'est pas attaché un nom classique. Les opérations d'addition et de produit sont effectuées par des variantes noncommutatives du résultant qui se ramènent à de l'algèbre linéaire élémentaire (Chapitre 14). Le test à zéro se réduit à tester l'égalité d'un nombre fini de conditions initiales. Une partie de ces résultats se résume comme suit :

**Proposition 1.10** Si $\mathbb{K}$ est un corps effectif, les séries formelles de $\mathbb{K}[[X]]$ qui sont solutions d'équations différentielles linéaires à coefficients dans $\mathbb{K}[X]$ forment un anneau effectif.

En d'autres termes, des structures de données finies permettent de manipuler ces
objets infinis et d'en tester l'égalité ou la nullité.

Ainsi, des identités élémentaires comme $\sin^2 X + \cos^2 X = 1$ sont non seulement facilement prouvables algorithmiquement, mais elles sont également calculables, c'est-à-dire que le membre droit se calcule à partir du membre gauche.

Les relations étroites entre équations différentielles linéaires et récurrences linéaires — les séries solutions des unes ont pour coefficients les solutions des autres — aboutissent aux mêmes réponses algorithmiques que les questions analogues sur des suites. Par exemple, l'identité de Cassini sur les nombres de Fibonacci

$$
F_{n+2}F_n - F_{n+1}^2 = (-1)^{n+1}, \quad n \ge 0
$$

est exactement du même niveau de difficulté que $\sin^2 X + \cos^2 X = 1$. Le pendant du
résultat précédent est donc la proposition suivante :

**Proposition 1.11** Si $\mathbb{K}$ est un corps effectif, l'ensemble des suites de $\mathbb{K}^\mathbb{N}$ solutions de récurrences linéaires à coefficients dans $\mathbb{K}[n]$ forme un anneau effectif.

Systèmes d'équations différentielles et de récurrences linéaires

Ces systèmes sont aux équations différentielles ou de récurrences ce que les systèmes polynomiaux sont aux polynômes en une variable. Les mêmes opérations sont disponibles. En particulier, l'élimination s'étend dans ce cadre en introduisant des algèbres d'opérateurs adaptés (Chapitre 31). Une application très importante, le télo scopage créatif, permet de calculer automatiquement des sommes et des intégrales définies. Ainsi,

$$
\sum_{k=0}^{n} \left( \sum_{j=0}^{k} \binom{n}{j} \right)^{3} = n 2^{3n-1} + 2^{3n} - 3 n 2^{n-2} \binom{2n}{n},
$$

$$
\sum_{n=0}^{\infty} H_n(x)H_n(y) \frac{u^n}{n!} = \frac{\exp\left(\frac{4u(xy-u(x^2+y^2))}{1-4u^2}\right)}{\sqrt{1-4u^2}},
$$

$$
\frac{1}{2} J_0(x)^2 + J_1(x)^2 + J_2(x)^2 + \dots = \frac{1}{2},
$$
$$ \int_{-1}^{+1} \frac{e^{-px} T_n(x)}{\sqrt{1-x^2}} dx = (-1)^n \pi I_n(p), $$

$$ \int_{0}^{+\infty} x e^{-px^2} J_n(bx) I_n(cx) dx = \frac{1}{2p} \exp\left(\frac{c^2-b^2}{4p}\right) J_n\left(\frac{bc}{2p}\right), $$

$$ \int_{0}^{+\infty} x J_{1}(ax) I_{1}(ax) Y_{0}(x) K_{0}(x) dx = - \frac{\ln(1-a^4)}{2\pi a^2}, $$

$$ \sum_{k=0}^{n} \frac{q^{k^2}}{(q;q)_k (q;q)_{n-k}} = \sum_{k=-n}^{n} \frac{(-1)^k q^{\frac{(5k^2-k)^2}{2}}}{(q;q)_{n-k} (q;q)_{n+k}}, $$

sont des formules qui mettent en jeu diverses fonctions spéciales ou polynômes orthogonaux classiques, et qui peuvent être prouvées automatiquement. Les algorithmes correspondants seront décrits au Chapitre 32. L'énoncé en termes d'anneaux effectifs est un peu lourd et omis ici.

En conclusion, les exemples ci-dessus illustrent bien la manière dont le calcul formel parvient à effectuer de nombreux calculs utiles dans les applications malgré l'indécidabilité révélée par le théorème de Richardson et Matiyasevich.

## 1.2 Calculer rapidement

En pratique, la calculabilité n'indique que la faisabilité. Il faut disposer d'algorithmes efficaces et des bonnes implantations pour pouvoir effectuer des calculs de grande taille. La première partie de cet ouvrage est consacrée aux algorithmes efficaces sur les structures de base du calcul formel. L'efficacité sera mesurée par la théorie de la complexité et nous ferons ressortir des principes récurrents dans la conception d'algorithmes efficaces.

**Exemple 1.2** Pour donner une idée de ce que veut dire rapidement, voici ce qui peut être calculé en *une seconde* avec le système MAPLE sur un ordinateur portable d'aujourd'hui<sup>a</sup>, en notant $\mathbb{K}$ le corps $\mathbb{Z}/p\mathbb{Z}$ à $p$ éléments, $p = 67108879$ étant un nombre premier de 26 bits (dont le carré tient sur un mot machine) :

1. Entiers :

- produit de deux entiers avec 30 000 000 chiffres;

- factorielle de 1 300 000 (environ 7 000 000 chiffres);

- factorisation d'un entier de 42 chiffres (produit de deux nombres premiers de taille la moitié).

2. Polynômes dans $\mathbb{K}[X]$ :

- produit de deux polynômes de degré 650 000;

- pgcd et résultant de deux polynômes de degré 12 500;

- factorisation d'un polynôme de degré 170 (produit comme ci-dessus).

3. Polynômes dans $\mathbb{K}[X, Y]$ :

- résultant de deux polynômes de degré total 20 (sortie de degré 400);
— factorisation d'un polynôme de degré 160 en deux variables.

4. Matrices :

— produit de deux matrices $850 \times 850$ à coefficients dans $K$;

— déterminant d'une matrice $1400 \times 1400$ à coefficients dans $K$;

— polynôme caractéristique d'une matrice $500 \times 500$ à coefficients dans $K$;

— déterminant d'une matrice $200 \times 200$ dont les coefficients sont des entiers 32 bits.

Ces exemples montrent qu'il est relativement aisé de calculer avec des objets de taille colossale. Dans la plupart de ces exemples, les algorithmes naifs mettraient plusieurs années pour le même calcul. On voit aussi qu'en une seconde, les algorithmes sont déjà dans leur régime asymptotique et qu'une analyse asymptotique suffira à obtenir des informations précises. Tout ceci donne envie d'introduire une mesure de complexité des différents algorithmes permettant d'expliquer, voire de prédire, les différences entre les tailles atteintes pour ces questions.

## Mesures de complexité

Pour bien définir la complexité, il faut se donner : un modèle de machine ; les opérations disponibles sur cette machine ; leur coût unitaire. La complexité en espace mesure la mémoire utilisée par l'exécution de l'algorithme, et la complexité en temps, la somme des coûts unitaires des opérations effectuées par l'algorithme. Dans cet ouvrage, nous n'aborderons pas la complexité en espace.

### Machine à accès direct

Le modèle que nous utiliserons est celui de la *machine à accès direct* (MAD), appelée aussi *Random Access Machine* (RAM) en anglais. Dans ce modèle, un programme lit et écrit des entiers sur deux bandes différentes et utilise un nombre arbitraire de registres entiers pour ses calculs intermédiaires. Les opérations élémentaires (l'assemblée de la machine) sont la lecture, l'écriture (sur bande ou en registre), l'addition, la soustraction, le produit, la division et trois instructions de saut : saut inconditionnel, saut si un registre est nul et saut si un registre est positif. Un point technique est que le programme ne fait pas partie des données, il n'est donc pas modifiable. Cet ouvrage ne rentre jamais dans les subtilités de l'utilisation de ce modèle, qui ne sert qu'à fixer précisément les mesures utilisées.

### Complexité arithmétique ou binaire

Nous considérerons deux mesures de complexité :

1. En calcul formel, de nombreux algorithmes peuvent être décrits comme opérant de façon abstraite sur une structure algébrique $\mathbb{A}$ effective. D'une façon concrète, la représentation des éléments de $\mathbb{A}$ importe peu et l'utilisation de chaque opération élémentaire de $\mathbb{A}$ (opération arithmétique ou prédicat) peut être parfaitement identifiée au cours de l'exécution de l'algorithme.

La *complexité arithmétique* d'un tel algorithme est définie comme le nombre total d'opérations arithmétiques et de tests de prédicats effectués dans $\mathbb{A}$. Cette mesure ne prend volontairement pas en compte les opérations de copie, celles

<footnote>a. Ce texte est écrit en 2015. Les ordres de grandeurs sont les mêmes sur différents ordinateurs et pour les différents systèmes de calcul formel principaux, avec des forces et des faiblesses différentes.</footnote>
sur les compteurs de boucle, les indirections, etc. En pratique cette mesure
reflète souvent le coût réel de l'algorithme si le coût des opérations dans $\mathbb{A}$ est
prépondérant et si chaque opération requiert un temps essentiellement constant.
C'est le cas par exemple si $\mathbb{A}$ est un corps fini. Ce n'est en revanche pas le cas
si $\mathbb{A}$ est le corps des nombres rationnels $\mathbb{Q}$ et que la taille de ces nombres croît
de façon significative durant les calculs.

2. Pour étudier le coût des algorithmes opérant sur des entiers, la complexité arithmétique avec $\mathbb{A} = \mathbb{Z}$ n'est pas pertinente. Il convient alors de décomposer les entiers dans une base $B$, qui est en pratique une puissance de 2 fixée (par exemple $2^{32}$). Chaque entier est alors vu comme un vecteur d'éléments de $\{0, \dots, B-1\}$. Nous appellerons complexité binaire la complexité arithmétique associée à $\mathbb{A} = \{0, \dots, B-1\}$, muni des tests d'égalité et comparaisons, ainsi que des opérations arithmétiques modulo $B$. En pratique cette mesure de complexité reflète en général bien les temps observés pour les algorithmes opérant sur des nombres entiers ou rationnels. Néanmoins, précisons que cette mesure de complexité ne correspond pas à celle utilisée classiquement sur une machine de Turing, puisqu'elle néglige les coûts d'accès mémoire. Par exemple le coût binaire pour transposer une matrice de taille $n \times n$ à coefficients dans $\mathbb{A} = \{0, \dots, B-1\}$ est nul avec notre définition : aucune opération arithmétique n'est requise.

La notation $O(\cdot)$

Nous utilisons la notation $O(\cdot)$ pour exprimer une borne sur la complexité des
algorithmes. La signification précise de la notation

$$
f(n) = O(g(n)), \quad n \to \infty
$$

est qu'il existe $K > 0$ et $A > 0$ tels que pour tout $n > A$, $f$ et $g$ sont liés par l'inégalité

$$
|f(n)| \leq K|g(n)|.
$$

Lorsque plusieurs paramètres interviennent dans l'analyse de la complexité, il faut
absolument préciser lequel tend vers l'infini pour que cette notation ait un sens. Si
plusieurs d'entre eux tendent vers l'infini, soit ils sont liés par des inégalités qui seront
précisées, soit la définition ci-dessus s'étend avec une constante $K$ qui ne dépend
d'aucun des paramètres.

Afin de simplifier les bornes de complexité nous utiliserons parfois la notation $\tilde{O}(\cdot)$ de sorte à cacher des facteurs logarithmiques. Plus précisément, nous écrirons

$$
f(n) = \tilde{O}(g(n)), \quad n \to \infty
$$

lorsque qu'il existe un entier $k \ge 0$ tel que

$$
f(n) = O(g(n) \log_2(\max(|g(n)|, 2))).
$$

La notation $O(\cdot)$ intervient aussi dans cet ouvrage pour représenter la troncature
des séries. L'expression

$$
f(X) := g(X) + O(X^N)
$$

signifiera que le polynôme ou la série $g$ est tronqué après son $N$-ième terme et que le
résultat, un polynôme, est stocké dans $f$.
**Exemple 1.3** L'addition de deux entiers de $n$ bits demande donc $O(1)$ opérations arithmétiques et $O(n)$ opérations binaires (quelque soit la base $B$, puisque la différence entre les bases joue sur la constante cachée dans le $O(\cdot)$).

**Exemple 1.4** Le calcul de $n!$ par la méthode naïve requiert $n$ opérations arithmétiques et $O(n^2 \log^2 n)$ opérations binaires. Le choix de la base $B$ n'intervient que dans la constante cachée dans le $O(\cdot)$. Nous verrons au Chapitre 15 qu'il est possible d'abaisser ce coût à seulement $O(n^{1/2} \log n)$ opérations arithmétiques dans $\mathbb{Z}$, et $O(n \log^3 n)$ opérations binaires. Les idées utilisées dans les algorithmes rapides permettant d'atteindre ces complexités fournissent aussi le meilleur algorithme connu de factorisation déterministe d'entiers et des algorithmes très efficaces pour le calcul de millions de décimales de $\pi$, $\log 2$ et de nombreuses autres constantes.

### Taille

Un algorithme et une structure de données sont généralement dotés d'une notion naturelle de taille et il s'agit d'étudier le coût de l'algorithme en fonction de cette taille. Pour simplifier, il est souvent commode de considérer le comportement asymptotique de ce coût lorsque la taille tend vers l'infini. Il est important de comprendre que la complexité d'un problème n'a de sens qu'une fois la structure de données fixée pour l'entrée comme pour la sortie.

Par exemple, pour les polynômes, le choix de la représentation dense conduit à mesurer la complexité par rapport au degré, alors que le choix de la représentation creuse met en avant le nombre de monômes. Pour la factorisation, la complexité est polynomiale en le degré, mais exponentielle en le nombre de monômes, dans le cas le pire.

### Cas le pire, cas moyen

La complexité dans le cas le pire est le maximum des complexités pour toutes les entrées d'une taille donnée. C'est celle que nous étudierons. Il est souvent utile de considérer aussi la complexité en moyenne, lorsque l'on peut mettre une mesure sur l'ensemble des entrées de taille bornée. Pour la plupart des algorithmes que nous étudierons dans la première partie de cet ouvrage, il n'y a pas de différence importante entre les deux. Ce n'est plus le cas en revanche pour la complexité des algorithmes de factorisation, ou pour ceux qui opèrent sur les systèmes polynomiaux.

### Bornes inférieures

La recherche de bornes inférieures de complexité est très difficile. Par exemple, à l'heure actuelle on ne sait pas prouver que la multiplication de matrices est nécessairement plus coûteuse qu'un nombre borné d'additions. Dès qu'il est possible de montrer que tous les bits de l'entrée doivent être pris en compte, la somme de la taille de l'entrée et de la taille de la sortie est une borne inférieure sur la complexité. En effet, dans le modèle MAD, chacune des écritures et des lectures prend une opération.

**Définition 1.2** Si $N$ est la somme de la taille de l'entrée et de la taille de la sortie, un algorithme sera dit *quasi-optimal* lorsque sa complexité est bornée par $\tilde{O}(N)$.
<figure><img src="image_4.png" /><figcaption>FIGURE 1.1 – Appels récursifs d’un algorithme « diviser pour régner ».</figcaption></figure>

L'essentiel de la première partie de cet ouvrage consistera à rechercher des algo-
rithmes quasi-optimaux pour les opérations de base sur les structures de données
fondamentales.

Diviser pour régner

Le principe le plus important de la conception d'algorithmes efficaces est la mé-
thodologie « diviser pour régner ». Il consiste à résoudre un problème en le réduisant
à un certain nombre *m* d'entrées de tailles divisées par *p* (le plus souvent *p* = 2) puis à
recombiner les résultats (voir Figure 1.1). Le coût de la recombinaison et éventuelle-
ment du découpage préliminaire est borné par une fonction *T* de la taille des entrées.
Lorsque les entrées ont une taille inférieure à *p* ou suffisamment petite, notée *s*, un
autre algorithme de coût *κ* indépendant de *n* est invoqué. Le coût total obéit alors
souvent à une récurrence de la forme

$$
C(n) \le \begin{cases} T(n) + mC(\lceil n/p \rceil), & \text{si } n \ge s(\ge p) \\ \kappa, & \text{sinon.} \end{cases} \tag{1.2}
$$

La notation ⌈x⌋ désigne l'entier $k$ tel que $k - 1 < x \le k$ (le « plafond » de $x$).

Qualitativement, le coût total de cette approche dépend fortement de la fonc-
tion T. Lorsque T est relativement élevée, les premiers niveaux de l'arbre contribuent
à l'essentiel du coût et, à une constante près, le coût est donc dominé par la première
étape de récursion. Les algorithmes à base d'itération de Newton du Chapitre 3 sont
de ce type. À l'inverse, pour une fonction T assez faible, le bas de l'arbre domine le
coût qui sera proportionnel au nombre de feuilles <fnref n="1" /> O((n/s)logₚ m). L'algorithme de

<footnote n="1">1. La notation $\log_p x$ représente le logarithme en base $p$ de $x$, c'est-à-dire $\log x/\log p$. L'identité facile à vérifier $a^{\log_p b} = b^{\log_p a}$ sera utile.</footnote>
multiplication de polynômes de Karatsuba du Chapitre 2 et l'algorithme de Strassen pour la multiplication de matrices (Chapitre 8) rentrent dans cette catégorie. Enfin, il est possible que tous les $O(\log_p(n/s))$ niveaux de l'arbre contribuent de manière assez équilibrée, menant à un coût en $O(T(n) \log n)$. La transformée de Fourier ra-pide (Chapitre 2) et le classique algorithme de tri fusion sont dans cette dernière catégorie.

Un cadre commode pour nos applications est résumé dans le lemme suivant, dont une version simplifiée que nous utiliserons dans la plupart des cas est donnée ensuite par le Théorème 1.13. Nous commençons par le cas simple où la taille est une puissance de $p$, où il est possible de donner des inégalités précises, ce qui prépare le passage au cas général et à l'asymptotique.

**Lemme 1.12 — Diviser pour régner.** Soit $C$ une fonction obéissant à l'inéga-lité (1.2) avec $m > 0$, $\kappa > 0$ et $T$ une fonction telle que

$$ T(pn) \ge qT(n), \quad n \in \mathbb{N}, \qquad (1.3) $$

avec $q > 1$. Alors, si $n$ est une puissance positive de $p$,

$$ C(n) \le \begin{cases} \left(1 - \frac{m}{q}\right)^{-1} T(n) + \kappa n^{\log_p m} & \text{si } q > m, \\ T(n) \log_p n + \kappa n^{\log_p q} & \text{si } q = m, \\ n^{\log_p m} \left( \frac{T(n)}{n^{\log_p q} \frac{q}{m-q}} + \kappa \right) & \text{si } q < m. \end{cases} $$

Le premier cas est celui où la complexité est portée par le haut de l'arbre, dans le deuxième elle l'est par tous les niveaux à parts égales, et dans le dernier par les bas niveaux.

*Démonstration.* L'utilisation répétée de l'inégalité (1.2) sur $C$ donne

$$ \begin{align} C(n) &\le T(n) + mC\left(\frac{n}{p}\right), \notag \\ &\le T(n) + mT\left(\frac{n}{p}\right) + \dots + m^{k-1}T\left(\frac{n}{p^{k-1}}\right) + m^k C\left(\frac{n}{p^k}\right), \tag{1.4} \end{align} $$

$$ \le T(n) \left( 1 + \frac{m}{q} + \dots + \left( \frac{m}{q} \right)^{k-1} \right) + m^k \kappa, \qquad (1.5) $$

où la dernière ligne résulte de (1.3) et du choix de $k = \lfloor \log_p \frac{n}{s} \rfloor + 1 \in [\log_p \frac{n}{s}, \log_p n]$. Ce choix entraine la suite d'inégalités suivante :

$$ 1 \le x \Rightarrow x^k \le x \cdot x^{\log_p(n/s)} = n^{\log_p x} \cdot x^{1-\log_p s} \le n^{\log_p x} \qquad (1.6) $$

qui permet de borner le deuxième terme de (1.5) avec $x = m$.

Si $m < q$, la somme entre parenthèses est majorée par la série géométrique. Si $m = q$,
la somme comporte $k$ termes égaux. Si $m > q$, récrire la somme sous la forme

$$ \left(\frac{m}{q}\right)^{k-1}\left(1+\frac{q}{m}+\cdots+\left(\frac{q}{m}\right)^{k-1}\right) $$

montre que le premier terme de (1.5) est majoré par

$$ T(n) \left(\frac{m}{q}\right)^k \frac{q}{m-q}. $$

Ensuite, l'inégalité (1.6) avec $x = m/q$ permet de conclure.

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

En pratique la régularité de T est souvent plus forte que celle donnée par (1.3) et on ne cherche qu'à comparer des estimations de $O(\cdot)$. Le résultat prend alors une forme plus simple, et sans restriction sur $n$ qui n'a plus besoin d'être une puissance de $p$.

**Théorème 1.13** — Diviser pour régner. Soit C une fonction obéissant à l'inégalité (1.2) avec $m > 0$ et $\kappa > 0$, et soit T une fonction croissante telle qu'existent $q$ et $r$ avec $1 < q \le r$ vérifiant

$$ qT(n) \le T(pn) \le rT(n), \quad \text{pour tout } n \text{ assez grand.} \qquad (1.7) $$

Alors, lorsque $n \to \infty$

$$ C(n) = \begin{cases} O(T(n)), & \text{si } q > m, \\ O(T(n) \log n), & \text{si } q = m, \\ O\left(n^{\log_p m} \frac{T(n)}{n^{\log_p q}}\right) & \text{si } q < m. \end{cases} $$

L'hypothèse (1.7) est en particulier vérifiée par les fonctions courantes $n^\alpha \log^\beta n$ pour les valeurs $\alpha > 0$.

*Démonstration.* Soit N la puissance de $p$ telle que $n \le N < pn$. En déroulant l'inégalité (1.2) on obtient une inégalité similaire à (1.4) où la division par $p$ est remplacée par l'opération $x \mapsto [x/p]$. La fonction T étant supposée croissante, toutes ces valeurs de T sont majorées par les valeurs en $N/p^i$, $i = 0, \dots, k-1$ et on obtient donc une majoration (1.5) où $n$ est remplacé par N. Le lemme précédent s'applique alors pour majorer cette expression, majoration qu'il ne reste plus qu'à majorer à son tour.

D'abord, par récurrence, l'inégalité (1.3) entraîne $T(n) \ge q^{\log_p n} T(1) = n^{\log_p q} T(1)$, ce qui permet d'absorber les termes en $\kappa$ dans les $O()$. Ensuite, la croissance de T et l'hypothèse (1.7) donnent pour $n$ assez grand $T(N) \le T(pn) = O(T(n))$, ce qui permet de conclure.

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /><figcaption>FIGURE 1.2 – Borne sur la complexité de l’algorithme de Karatsuba : en noir, la suite solution de $C(n) = 4n + 3C(\lceil n/2 \rceil)$, $C(1) = 1$ ; en rouge, la borne obtenue à l’Exemple 1.5.</figcaption></figure>

**Exemple 1.5** La complexité de l’algorithme de Karatsuba du Chapitre 2 vérifie la récurrence

$$K(n) \le 4n + 3K(\lceil n/2 \rceil), \quad K(1) = 1.$$

Le Lemme 1.12 et le Théorème 1.13 s'appliquent avec $m = 3$, $p = s = 2$, $\kappa = 1$ et T linéaire, si bien que $q = r = 2$, d'où une complexité $K(n) = O(n^{\log_2 3})$, et plus précisément

$$K(n) \le n^{\log_2 3} \left( \frac{4n}{\frac{n}{1} + 1} \right) = 9n^{\log_2 3}$$

pour $n$ une puissance de 2, et donc en majorant par la valeur en la puissance de 2 immédiatement supérieure, le cas général est borné par

$$K(n) \le 9 \cdot 2^{\lceil \log_2 n \rceil \log_2 3} = 9 \cdot 3^{\lceil \log_2 n \rceil}.$$

Cette borne est assez fine, comme le montre la Figure 1.2.

## Conventions

Pour toute la suite de cet ouvrage, et sauf mention contraire, nous adoptons les notations et terminologies suivantes :

- tous les anneaux et corps sont supposés effectifs, commutatifs et unitaires;

- les estimations de complexité sont en nombre d'opérations arithmétiques dans
l'anneau de base;

— le symbole $\mathbf{\blacksquare}$ marque la fin d'une démonstration;

— les logarithmes sont en base 2, lorsque la base n'est pas spécifiée;

— le lemme « diviser pour régner » fait référence à l'énoncé du Lemme 1.12;

— le théorème « diviser pour régner » fait référence à l'énoncé du Théorème 1.13.

## Notes

Le théorème de Richardson–Matiyasevich que nous énonçons a ses origines dans les travaux sur le 10ᵉ problème de Hilbert portant sur l'existence d'un algorithme permettant de trouver les racines entières des polynômes multivariés à coefficients entiers. Une étape importante dans la résolution de ce problème est due à Davis, Putnam et Robinson [DPR61] qui montrent en 1961 l'indécidabilité de ce problème si en plus des polynômes on considère la fonction $x \mapsto 2^x$. Quelques années plus tard, Richardson [Ric68] ramène les problèmes de racines réelles en plusieurs variables à ces questions, et, à l'aide d'un codage astucieux par le sinus, y ramène aussi les problèmes en une variable. Par rapport à notre énoncé, son théorème fait intervenir la fonction $\exp$ et la constante $\ln 2$ (pour utiliser le résultat de Davis, Putnam et Robinson), ainsi que la constante $\pi$ (pour utiliser la nullité du sinus aux multiples entiers de $\pi$). Deux ans plus tard, en 1970, Matiyasevich prouve que la fonction $x \mapsto 2^x$ elle-même peut s'exprimer à l'aide de polynômes, concluant ainsi la preuve de l'indécidabilité du 10ᵉ problème de Hilbert. Le théorème de Richardson peut alors être simplifié, et c'est ce que fait Matiyasevich dans son livre [Mat93], où il montre aussi comment se débarrasser de la constante $\pi$. Il faut noter que par contre, si l'on remplace le sinus par l'exponentielle, alors la simplification devient décidable [Ric69].

Ces théorèmes s'appliquent à des fonctions. Pour des constantes, l'approche la plus récente réduit le test à zéro à une conjecture de théorie des nombres due à Schanuel qui exprime que les seules relations entre exponentielles et logarithmes sont celles qui découlent des formules d'addition et de multiplication. Si l'on accepte cette conjecture de Schanuel, alors un résultat relativement récent de Macintyre et Wilkie en 1996 [MW96] entraîne l'existence d'un algorithme de reconnaissance de 0 pour la classe des constantes obtenues à partir de 1 par addition, soustraction, multiplication et exponentielle. Un tel algorithme, à base d'évaluation numérique et de LLL (voir Chapitre 20) a ensuite été donné en 1997, à nouveau par Richardson [Ric97].

Les différents modèles de complexité (machine MAD, straight-line program, machine de Turing, ...) sont bien présentés par Aho, Hopcroft et Ullman dans leur livre [AHU74], ou par Stern [Ste94]. La jolie identité (1.1) est tirée de [Bec+05]. Les résultats de complexité autour de l'algorithmique du « diviser pour régner » sont nombreux. Le premier remonte sans doute à Bentley, Haken et Saxe [BHS80] qui font apparaître la série géométrique. Une version moderne tenant compte des planchers et des plafonds est présentée par Cormen, Leiserson, Rivest et Stein qui l'appellent *Master Theorem* [Cor+09]. Des variantes sophistiquées sont connues [DS11; Rou01]. Dans cette direction, Yap [Yap11] donne un théorème prêt à l'emploi dans le cas où plusieurs divisions avec des valeurs de $p$ distinctes ont lieu.
Bibliographie

AHU74 AHO, Alfred V., John E. HOPCROFT et Jeffrey D. ULLMAN (1974). *The design and analysis of computer algorithms*. Addison-Wesley Publishing Co.

Bec+05 BECK, Matthias, Bruce C. BERNDT, O-Yeat CHAN et Alexandru ZAHARESCU (2005). « Determinations of analogues of Gauss sums and other trigonometric sums ». In : *International Journal of Number Theory*, vol. 1, n°3, p. 333-356.

BHS80 BENTLEY, Jon Louis, Dorothea HAKEN et James B. SAXE (1980). « A general method for solving divide-and-conquer recurrences ». In : *SIGACT News*, vol. 12, n°3, p. 36-44.

Cor+09 CORMEN, Thomas H., Charles E. LEISERSON, Ronald L. RIVEST et Clifford STEIN (2009). *Introduction to algorithms*. Third. MIT Press, Cambridge, MA.

DPR61 DAVIS, Martin, Hilary PUTNAM et Julia ROBINSON (1961). « The decision problem for exponential diophantine equations ». In : *Annals of Mathematics. Second Series*, vol. 74, p. 425-436.

DS11 DRMOTA, Michael et Wojciech SZPANKOWSKI (2011). « A master theorem for discrete divide and conquer recurrences ». In : *SODA'11 : ACM-SIAM Symposium on Discrete Algorithms*. San Francisco, CA : SIAM, p. 342-361.

Mat93 MATIYASEVICH, Yuri V. (1993). *Hilbert's tenth problem*. Foundations of Computing Series. Translated from the 1993 Russian original by the author, with a foreword by Martin Davis. MIT Press.

MW96 MACINTYRE, Angus et A. J. WILKIE (1996). « On the decidability of the real exponential field ». In : *Kreiseliana*. A. K. Peters, p. 441-467.

Ric68 RICHARDSON, Daniel (1968). « Some undecidable problems involving elementary functions of a real variable ». In : *Journal of Symbolic Logic*, vol. 33, n°4, p. 514-520.

Ric69 — (1969). « Solution of the identity problem for integral exponential functions ». In : *Zeitschrift fur mathematischen Logik und Grundlagen der Mathematik*, vol. 15, p. 333-340.

Ric97 — (1997). « How to recognize zero ». In : *Journal of Symbolic Computation*, vol. 24, n°6, p. 627-645.

Rou01 ROURA, Salvador (2001). « Improved master theorems for divide-and-conquer recurrences ». In : *Journal of the Association for Computing Machinery*, vol. 48, n°2, p. 170-205.

Ste94 STERN, Jacques (1994). *Fondements mathématiques de l'informatique*. Ediscience international.

Yap11 YAP, Chee (2011). « A real elementary approach to the master recurrence and generalizations ». In : *Theory and applications of models of computation*. Vol. 6648. Lecture Notes in Computer Science. Tokyo, Japan : Springer-Verlag, p. 14-26.
<figure><img src="image_8.png" /></figure>
# Polynômes et séries

<table>
  <tbody>
    <tr>
      <td>2</td>
      <td>Multiplication rapide</td>
      <td>35</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>Introduction, résultats principaux</td>
      <td></td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Algorithme naïf</td>
      <td></td>
    </tr>
    <tr>
      <td>2.3</td>
      <td>Algorithme de Karatsuba</td>
      <td></td>
    </tr>
    <tr>
      <td>2.4</td>
      <td>Transformée de Fourier rapide</td>
      <td></td>
    </tr>
    <tr>
      <td>2.5</td>
      <td>L'algorithme de Schönhage-Strassen</td>
      <td></td>
    </tr>
    <tr>
      <td>2.6</td>
      <td>Algorithmes pour les entiers</td>
      <td></td>
    </tr>
    <tr>
      <td>2.7</td>
      <td>Un concept important : les fonctions de multiplication</td>
      <td></td>
    </tr>
    <tr>
      <td>3</td>
      <td>Calculs rapides sur les séries</td>
      <td>57</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Séries formelles</td>
      <td></td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>La méthode de Newton pour le calcul d'inverses</td>
      <td></td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Itération de Newton formelle et applications</td>
      <td></td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>La composition des séries</td>
      <td></td>
    </tr>
    <tr>
      <td>4</td>
      <td>Division euclidienne</td>
      <td>81</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Introduction</td>
      <td></td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Division de polynômes</td>
      <td></td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Suites récurrentes linéaires à coefficients constants</td>
      <td></td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Développement de fractions rationnelles</td>
      <td></td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Applications</td>
      <td></td>
    </tr>
    <tr>
      <td>5</td>
      <td>Calculs modulaires, évaluation et interpolation</td>
      <td>95</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Introduction</td>
      <td></td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Présentation, résultats</td>
      <td></td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Interpolation de Lagrange</td>
      <td></td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>Algorithmes rapides</td>
      <td></td>
    </tr>
    <tr>
      <td>6</td>
      <td>Pgcd et résultant</td>
      <td>111</td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Algorithme d'Euclide</td>
      <td></td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Résultant</td>
      <td></td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>Algorithme d'Euclide rapide</td>
      <td></td>
    </tr>
    <tr>
      <td>7</td>
      <td>Approximants de Padé et de Padé-Hermite</td>
      <td>139</td>
    </tr>
    <tr>
      <td>7.1</td>
      <td>Reconstruction rationnelle</td>
      <td></td>
    </tr>
    <tr>
      <td>7.2</td>
      <td>Approximants de Padé-Hermite</td>
      <td></td>
    </tr>
  </tbody>
</table>
<figure><img src="image_9.png" /></figure>
## 2. Multiplication rapide

### Résumé

Les algorithmes rapides de calcul du produit de polynômes et d'entiers sont au cœur de l'algorithmique efficace en calcul formel. La plupart des gains de complexité dans les chapitres ultérieurs reposent sur l'efficacité de la multiplication. Pour multiplier deux polynômes de degré $n$ à coefficients dans un anneau $\mathbb{A}$, la méthode classique requiert $O(n^2)$ opérations dans $\mathbb{A}$. De même, l'algorithmme scolaire de multiplication de deux entiers à $n$ chiffres nécessite un nombre d'opérations binaires en $O(n^2)$. Nous présentons dans ce chapitre plusieurs algorithmes de multiplication rapide, dont celui de Karatsuba, de complexité $O(n^{1.59})$, ainsi que ceux utilisant la transformée de Fourier rapide, dont la complexité est quasiment linéaire en $n$.

### 2.1 Introduction, résultats principaux

Les problèmes abordés dans ce chapitre concernent la complexité arithmétique de la multiplication des polynômes à une variable et la complexité binaire de la multiplication des entiers. Au vu de l'exemple suivant, il est facile de se convaincre de la similitude des deux questions :

**Polynômes** Soient à multiplier $3X^2 + 2X + 1$ et $6X^2 + 5X + 4$ dans $\mathbb{Z}[X]$,

$$ 
\begin{aligned}
& (3X^2 + 2X + 1) \times (6X^2 + 5X + 4) \\ 
&= (3 \cdot 6)X^4 + (3 \cdot 5 + 2 \cdot 6)X^3 + (3 \cdot 4 + 2 \cdot 5 + 1 \cdot 6)X^2 + (2 \cdot 4 + 1 \cdot 5)X + (1 \cdot 4) \\ 
&= 18X^4 + 27X^3 + 28X^2 + 13X + 4.
\end{aligned}
$$
**Nombres entiers** Soient à multiplier 321 et 654 en base 10,

$$ 
\begin{aligned}
& (3 \cdot 10^2 + 2 \cdot 10 + 1) \times (6 \cdot 10^2 + 5 \cdot 10 + 4) \\
&= (3 \cdot 6)10^4 + (3 \cdot 5 + 2 \cdot 6)10^3 + (3 \cdot 4 + 2 \cdot 5 + 1 \cdot 6)10^2 + (2 \cdot 4 + 1 \cdot 5)10 + (1 \cdot 4) \\
&= 18 \cdot 10^4 + 27 \cdot 10^3 + 28 \cdot 10^2 + 13 \cdot 10 + 4 \\
&= 2 \cdot 10^5 + 9 \cdot 10^3 + 9 \cdot 10^2 + 3 \cdot 10 + 4 = 209934.
\end{aligned}
$$

Dans les deux cas, nous avons retranscrit l'algorithmme naïf, et la suite des calculs est essentiellement la même, si ce n'est que, dans le cas des entiers, il faut en outre gérer les retenues (dernière égalité de l'exemple). On ne sera donc pas surpris que les résultats obtenus dans les deux cas soient très semblables.

### Résultats

Dans toute la suite, $(\mathbb{A}, +, \times)$ désignera un anneau commutatif effectif. Tout d'abord, nous considérons la complexité arithmétique; il s'agit de minimiser le nombre d'opérations $(+, -, \times)$ dans $\mathbb{A}$ pour multiplier des polynômes en degré borné. Les premiers résultats à retenir de ce chapitre sont les suivants.

**Théorème 2.1** La multiplication des polynômes de degré au plus $n$ dans $\mathbb{A}[X]$ requiert :

- $O(n^2)$ opérations dans $\mathbb{A}$ par l'algorithmme naïf;

- $O(n^{1,59})$ opérations dans $\mathbb{A}$ par l'algorithmme de Karatsuba;

- $O(n \log n \log \log n)$ — voire dans certains cas $O(n \log n)$ — opérations dans $\mathbb{A}$, via la transformée de Fourier rapide (FFT).

Ainsi, la multiplication des polynômes peut se faire en un coût arithmétique quasi-linéaire en leur degré.

La multiplication des polynômes est omniprésente : les algorithmes de calcul de pgcd (plus grand commun diviseur), de pgcd étendu, de factorisation en une ou plusieurs variables, de composition des séries formelles, d'évaluation multipoint, d'interpolation, font tous intervenir des produits de polynômes.

L'analogie entre les entiers et les polynômes va très loin; la plupart des réponses apportées dans le cadre de la complexité arithmétique trouvent un équivalent en complexité binaire. Cependant, aucun théorème d'équivalence n'est connu; il se trouve que les mêmes idées algorithmiques s'adaptent plus ou moins facilement dans les deux cadres. Ainsi, on dispose des résultats suivants dans le modèle binaire.

**Théorème 2.2** On peut multiplier des entiers de $n$ chiffres binaires par :

- l'algorithmme naïf en $O(n^2)$ opérations binaires;

- l'algorithmme de Karatsuba en $O(n^{1,59})$ opérations binaires;

- l'algorithmme de Schönhage-Strassen en $O(n \log n \log \log n)$ opérations binaires.

Les preuves de ces résultats de complexité binaire sont plus délicates que celles de leurs analogues polynomiaux, à cause des problèmes de gestion des retenues. Aussi dans la suite ne traitons nous d'abord en détail que les versions polynomiales de
ces résultats (Théorème 2.1), pour passer ensuite brièvement en revue le cas entier (Théorème 2.2).

En pratique

Les constantes cachées dans les O(·) sont déterminantes pour l'efficacité pratique de tels algorithmes. Par exemple, lorsque $\mathcal{A}$ est un corps fini de taille « raisonnable » (typiquement, dont les éléments sont représentés sur quelques mots machine), pour le produit de polynômes dans les meilleures implantations actuelles (MAGMA, NTL, FLINT, MATHEMAGIX) :

— l'algorithme de Karatsuba bat l'algorithme naïf pour des degrés d'environ 20;

— les méthodes à base de FFT en $O(n \log n)$ battent l'algorithme de Karatsuba pour des degrés de l'ordre de 100, mais ne peuvent pas être utilisées pour des degrés arbitrairement grands (vient un moment où l'on manque de racines de l'unité, voir plus loin);

— l'algorithme de type FFT en $O(n \log n \log \log n)$ est utilisé pour des degrés de l'ordre de quelques dizaines ou centaines de milliers.

Certains problèmes, en cryptologie ou en théorie des nombres, nécessitent de manipuler des polynômes de degré de l'ordre de 100 000, tailles pour lesquelles les algorithmes rapides sont indispensables. Plus fréquemment, des degrés de l'ordre de la centaine ou du millier font partie du quotidien du calculateur formel (au moins dans des calculs intermédiaires).

Ces algorithmes ne sont pas concurrents, mais complémentaires. Les bonnes implantations passent automatiquement à l'algorithme le plus efficace en fonction du degré, et ce parfois même au cours d'un même calcul, lorsqu'un algorithme récursif comme celui de Karatsuba atteint une taille pour laquelle il vaut mieux utiliser l'algorithme naïf.

La situation pour les entiers est similaire, même si l'implantation des algorithmes rapides pour les entiers est bien plus délicate en raison des retenues. Dans les meilleures implantations actuelles (MAGMA, GMP) :

— l'algorithme de Karatsuba bat l'algorithme naïf pour des nombres de l'ordre de 30 chiffres décimaux;

— les méthodes à base de FFT (Schönhage–Strassen) gagnent pour des nombres d'environ 3 000 chiffres décimaux.

À nouveau, des entiers de quelques centaines de chiffres sont courants dès que l'on sort des exemples jouets. Au-delà, des problèmes venus de la cryptologie ou de la théorie des nombres demandent de manipuler des nombres de taille colossale (de l'ordre de 100 000 000 chiffres ; il faut 10 Mo pour stocker un tel nombre). Ceci justifie amplement les efforts d'implantation d'algorithmes rapides.

Dans la suite de ce chapitre, on travaille avec des polynômes F et G à coefficients dans un anneau $\mathcal{A}$, ayant un degré au plus $n-1$, formellement

$$F = f_0 + \cdots + f_{n-1} X^{n-1} \quad \text{et} \quad G = g_0 + \cdots + g_{n-1} X^{n-1};$$

le problème est alors de calculer (les coefficients de)

$$H = FG = h_0 + \cdots + h_{2n-2} X^{2n-2}.$$
### 2.2 Algorithme naïf

L'algorithme naïf consiste à développer le produit, c'est-à-dire à écrire

$$H = FG = \sum_{i=0}^{2n-2} h_i X^i \quad \text{avec} \quad h_i = \sum_{j+k=i} f_j g_k.$$

Ainsi, calculer tous les $h_i$ demande $O(n^2)$ opérations dans $\mathbb{A}$. C'est un algorithme de complexité arithmétique quadratique.

**Exercice 2.1** Montrer que, pour multiplier deux polynômes de degrés $m$ et $n$, l'algorithme naïf demande au plus $(m+1) \times (n+1)$ multiplications dans $\mathbb{A}$ et $mn$ additions dans $\mathbb{A}$. ■

**Exercice 2.2** Montrer que, pour multiplier deux entiers à $n$ chiffres chacun, l'algorithme naïf demande $O(n^2)$ opérations binaires. ■

**Exercice 2.3** Estimer la complexité binaire de la méthode naïve lorsque les polynômes ont degré $n$ et des coefficients entiers bornés en valeur absolue par un entier $H$. ■

**Exercice 2.4** Alexander J. Yee et Shigeru Kondo ont calculé $10^{13}$ décimales de $\pi$ sur un PC de bureau en 2011ª. Ce type de calcul repose sur la multiplication d'entiers de même taille que le nombre de décimales cherchées. En supposant que leur machine était capable d'effectuer $10^{12}$ opérations à la seconde, montrer qu'ils n'ont pas utilisé l'algorithme naïf. ■

### 2.3 Algorithme de Karatsuba

Un premier raffinement de l'algorithme naïf repose sur la remarque suivante : il est possible de gagner une multiplication pour le produit des polynômes de degré 1, parmi les 4 du produit par l'algorithme quadratique. Soient en effet à multiplier les polynômes

$$F = f_0 + f_1 X \quad \text{et} \quad G = g_0 + g_1 X.$$

Leur produit $H = FG = h_0 + h_1 X + h_2 X^2$ peut être obtenu par une forme d'interpolation sur les points $0,1,\infty$: en 0 on a $H(0) = h_0 = f_0 g_0$; de même, le coefficient de plus haut degré (qui correspond intuitivement à l'évaluation en l'infini) vaut $h_2 = f_1 g_1$. Enfin, la valeur en 1 vaut $h_0 + h_1 + h_2 = F(1)G(1) = (f_0 + f_1)(g_0 + g_1)$. Ainsi, on obtient $h_1 = (f_0 + f_1)(g_0 + g_1) - h_0 - h_2$ pour seulement une multiplication supplémentaire, donc l'ensemble des coefficients du produit pour 3 multiplications et 4 additions. Quelques additions sont perdues par rapport à l'algorithme naïf, mais le gain d'une multiplication va se transformer en gain dans l'exposant de l'algorithme, par application récursive.

En effet, dans le cas général des degrés quelconques, il suffit de scinder $F$ et $G$ en deux et de procéder de la même manière. Si $F$ et $G$ sont de degré au plus $n-1$, avec

<footnote>a. Voir http://www.numberworld.org/misc_runs/pi-5t/announce_en.html.</footnote>
**Entrée** F, G de degrés au plus $n-1$.

**Sortie** H = FG.

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

2. Décomposer F et G selon l'équation (2.1).

3. Calculer $A_1 = F^{(0)}G^{(0)}$ et $A_2 = F^{(1)}G^{(1)}$ récursivement.

4. Calculer $A_3 = F^{(0)} + F^{(1)}$ et $A_4 = G^{(0)} + G^{(1)}$.

5. Calculer $A_5 = A_3A_4$ récursivement.

6. Calculer $A_6 = A_5 - A_1$ et $A_7 = A_6 - A_2$.

7. Renvoyer $A_1 + A_7X^k + A_2X^{2k}$.

**Algorithme 2.1** – Multiplication de polynômes par l'algorithme de Karatsuba.

$k = ⌊n/2⌋$, on pose

$$F = F^{(0)} + F^{(1)}X^k, \quad G = G^{(0)} + G^{(1)}X^k, \qquad (2.1)$$

pour des polynômes $F^{(0)}$, $F^{(1)}$, $G^{(0)}$, $G^{(1)}$ de degrés au plus $k-1$. Le produit $H = FG$ s'écrit

$$H = F^{(0)}G^{(0)} + (F^{(0)}G^{(1)} + F^{(1)}G^{(0)})X^k + F^{(1)}G^{(1)}X^{2k}.$$

Cette méthode est résumée dans l'Algorithme 2.1. Le théorème qui suit établit la complexité de l'algorithme.

**Théorème 2.3** Si $n$ est une puissance de 2, l'algorithme de Karatsuba calcule le produit de deux polynômes de degrés au plus $n-1$ en au plus $9n^{\log_2 3}$ opérations dans $\mathbb{A}$.

*Démonstration.* Un appel en degré strictement inférieur à $n$ effectue trois appels récurrents en degré strictement inférieur à $n/2$, plus quelques additions, comptées comme suit : l'étape (4) effectue deux additions en taille strictement inférieure à $n/2$; l'étape (6) effectue deux additions en taille strictement inférieure à $n$; quant à l'étape (7), $A_1$ et $A_2X^n$ sont à supports monomiaux distincts, donnant la somme $A_1 + A_2X^n$ sans coût arithmétique, puis la somme totale en une addition en taille strictement inférieure à $n$. Le coût $K(n)$ satisfait donc à la récurrence

$$K(n) \leq 3K(n/2) + 4n,$$

où le terme $4n$ vient compter le nombre d'additions dans $\mathbb{A}$. Le lemme « diviser pour régner » permet alors de conclure avec $p = q = s = 2$, $m = 3$, $T(n) = 4n$, et $\kappa = 1$ (coût de l'algorithme naïf en degré 0).

Le théorème « diviser pour régner » fournit ensuite :

**Corollaire 2.4** On peut multiplier deux polynômes de degré $n$ arbitraire en $O(n^{\log_2 3}) = O(n^{1.59})$ opérations dans $\mathbb{A}$.

Une manière d'estimer plus précisément la constante cachée dans cette borne a été présentée au Chapitre 1.
**Exercice 2.5** Soit $n$ une puissance de 2. Établir un algorithme hybride de multiplication dans $\mathbb{A}[X]$, qui fait appel à l'algorithme de Karatsuba pour $n > 2^d$ et à l'algorithme naïf pour $n \le 2^d$. Montrer que la complexité arithmétique $C(n)$ de cet algorithme vérifie $C(n) \le \gamma(d)n^{\log_2 3 - 8n}$ pour tout $n \ge 2^d$, où $\gamma(d)$ est une fonction qui dépend uniquement de $d$. Trouver la valeur de $d$ qui minimise $\gamma(d)$ et, par comparaison avec le résultat du Théorème 2.3, estimer le gain obtenu par cette optimisation. ■

**Exercice 2.6** Est-il possible de multiplier par un algorithme universel, c'est-à-dire indépendant de l'anneau de base $\mathbb{A}$, deux polynômes de degré au plus 1 en utilisant seulement 2 multiplications dans $\mathbb{A}$ ? ■

**Exercice 2.7 — Algorithme de Toom-Cook.** Soit $\mathbb{A}$ un anneau et soient $A$ et $B$ deux polynômes de degré au plus 3 dans $\mathbb{A}[X]$.

1. Estimer le nombre de multiplications de $\mathbb{A}$ requises par l'algorithme de Karatsuba pour calculer le produit $AB$.

2. On suppose que 2, 3 et 5 sont inversibles dans $\mathbb{A}$ et que la division d'un élément de $\mathbb{A}$ par 2, 3 et 5 est gratuite. Donner un algorithme qui multiplie $A$ et $B$ en utilisant au plus 7 multiplications dans $\mathbb{A}$.

3. On suppose que 2, 3 et 5 sont inversibles dans $\mathbb{A}$. Donner un algorithme de multiplication polynomiale dans $\mathbb{A}[X]$ de complexité arithmétique $O(n^{1,4})$.

Dans la suite de l'exercice, on suppose que l'anneau $\mathbb{A}$ est de caractéristique nulle.

4. Montrer que, pour tout entier $\alpha \ge 2$, il existe un algorithme de multiplication polynomiale dans $\mathbb{A}[X]$ de complexité arithmétique $O(n^{\log_{\alpha}(2\alpha-1)})$.

5. Montrer que pour tout $\varepsilon > 0$, il existe un algorithme de multiplication polynomiale dans $\mathbb{A}[X]$ de complexité arithmétique $O(n^{1+\varepsilon})$, où la constante dans $\O(\cdot)$ dépend de $\varepsilon$, mais pas de $n$.

## 2.4 Transformée de Fourier rapide

Les méthodes à base de transformée de Fourier rapide (appelée aussi FFT pour *Fast Fourier Transform*) sont ce que l'on sait faire de mieux, à l'heure actuelle, pour multiplier les polynômes. Pour simplifier la présentation, on suppose ici que l'on cherche à multiplier des polynômes $F$ et $G$ dans $\mathbb{A}[X]$, de degrés strictement inférieurs à $n/2$ (ou plus généralement tels que $\deg(FG) < n$).

### Idée de l’algorithme

En supposant que l'anneau $\mathbb{A}$ le permette, l'idée générale est synthétisée dans l'Algorithme 2.2. Il s'agit d'évaluer en des points bien choisis, de multiplier les évaluations, et de reconstruire les coefficients du produit à partir de ces valeurs (à condition que cette opération d'interpolation soit possible, voir ci-dessous). Si deux polynômes coïncident sur $1, \omega, \dots, \omega^{n-1}$, nous verrons que leur différence est un multiple de $X^n - 1$, ce qui justifie la correction de l'algorithme.

Lorsque l'algorithme est employé avec l'hypothèse $\deg H < n$, les coefficients de $H$ mod $X^n - 1$ qui sont renvoyés sont bien ceux de $H$. Le coût des étapes de précalcul et de produit point à point est linéaire en $n$, et il reste donc à voir comment effectuer rapidement les opérations d'évaluation et d'interpolation.
**Entrée** F et G deux polynômes, $n$ un entier, et $\omega$ une racine principale $n$-ième de l'unité.

**Sortie** $H = FG \bmod X^n - 1$.

1. Précalcul. Calculer les puissances $\omega^2, \dots, \omega^{n-1}$.

2. Évaluation. Calculer les valeurs :

$$ \mathrm{Ev}(F) = (F(\omega^0), \dots, F(\omega^{n-1})); \quad \mathrm{Ev}(G) = (G(\omega^0), \dots, G(\omega^{n-1})). $$

3. Produit point à point.

$$ (\mathrm{Ev}(F), \mathrm{Ev}(G)) \mapsto \mathrm{Ev}(FG) = (FG(\omega^0), \dots, FG(\omega^{n-1})). $$

4. Interpolation.

$$ \mathrm{Ev}(FG) \mapsto FG. $$

**Algorithme 2.2** – Multiplication de polynômes par transformée de Fourier discrète.

## Racines primitives et principales de l'unité

**Définition 2.1** L'élément $\omega$ de $\mathbb{A}$ est une racine $n$-ième de l'unité si $\omega^n = 1$ ; c'est une racine $n$-ième primitive de l'unité si de plus $\omega^t \neq 1$ pour $t < n$ ; c'est une racine $n$-ième principale de l'unité si de plus $\omega^t - 1$ est non diviseur de zéro dans $\mathbb{A}$ pour $t \in \{1, \dots, n-1\}$ (c'est-à-dire que $\alpha(\omega^t - 1) = 0$ implique $\alpha = 0$).

**Exemple 2.1** Si $\mathbb{A}$ est un corps ou même un anneau intègre, les racines principales et primitives coïncident et la condition revient à dire que $\omega$ engendre le groupe des racines $n$-ièmes de l'unité. Par exemple, dans $\mathbb{C}$, $-1$ n'est pas une racine $4^e$ primitive de l'unité, alors que $i$ l'est. Plus généralement, lorsque $\mathbb{A} = \mathbb{C}$, les racines primitives $n$-ièmes de l'unité sont de la forme $\exp(2qi\pi/n)$, pour $q$ premier avec $n$.

**Exemple 2.2** Dans $\mathbb{A} = \mathbb{Z}/25\mathbb{Z}$, les premières puissances de 6 sont 6, 11, 16, 21, et 1. Comme par exemple $16 - 1 = 15$ est diviseur de 0, 6 est une racine $5^e$ de l'unité qui est primitive mais n'est pas principale. En revanche, les premières puissances de 7 sont 7, 24, 18, 1, ce qui montre que 7 est une racine principale $4^e$ de 1 dans $\mathbb{A}$.

L'intérêt des racines principales de 1 dans des anneaux qui ne sont même pas intègres apparaît en fin de ce chapitre pour l'algorithme de Schönhage et Strassen. Jusque là, il est suffisant de considérer le cas où $\mathbb{A} = \mathbb{C}$ pour comprendre les idées. Les propriétés que nous utiliserons sont résumées dans le lemme suivant.

**Lemme 2.5** Si $\omega$ est racine primitive ou principale $n$-ième de l'unité, alors
1. $\omega^{-1}$ aussi ;
2. si $n = pq$ alors $\omega^p$ est une racine q-ième de l'unité de même nature que $\omega$;

3. pour $\ell \in \{1, \dots, n-1\}$, et $\omega$ une racine principale de l'unité alors

$$ \sum_{j=0}^{n-1} \omega^{\ell j} = 0. $$

*Démonstration.* Tout d'abord $\omega$ est bien inversible : l'identité $\omega^n = 1$ montre que $\omega^{n-1}$ est un inverse de $\omega$. Ensuite, $\omega^{-1}$ est une racine de l'unité : le produit de l'identité précédente par $\omega^{-n}$ donne $1 = \omega^{-n}$. Enfin, elle est principale si $\omega$ l'est : lorsque $\omega^t - 1$ n'est pas diviseur de 0, son produit par l'inversible $\omega^{-t}$ non plus. Le même argument s'applique dans le cas où $\omega$ est primitive.

La deuxième propriété est immédiate. Pour la troisième, le produit de la somme par $1 - \omega^\ell$ télescope les sommants, ce qui donne

$$ (1 - \omega^\ell) \sum_{j=0}^{n-1} \omega^{\ell j} = 1 - \omega^{\ell n} = 1 - (\omega^n)^\ell = 0. $$

Comme $1 - \omega^\ell$ n'est pas diviseur de 0, la somme est bien nulle. ■

R

Dans la définition de racine primitive ou principale, la condition $t$ diviseur strict de $n$ suffit, elle entraîne la propriété pour $t \in \{1, \dots, n-1\}$. En effet, si $t$ n'est pas diviseur de $n$, son pgcd $g$ avec $n$ l'est. Il existe d'après le théorème de Bézout deux entiers $p, q \in \mathbb{N}$ tels que $g = tp + nq$. Alors, l'égalité $\alpha(\omega^t - 1) = 0$ entraîne $0 = \alpha(\omega^t - 1)(1 + \omega^p + \dots + \omega^{t(p-1)}) = \alpha(\omega^{tp} - 1) = \alpha(\omega^{g-nq} - 1) = \alpha(\omega^g - 1)$, et donc $\alpha = 0$, puisque $g$ est un diviseur strict de $n$.

### Transformée de Fourier rapide

L'opération

$$ \text{DFT}: F \in \mathbb{A}[X] \mapsto (F(1), F(\omega), \dots, F(\omega^{n-1})), $$

où $\omega$ est une racine principale $n$-ième de l'unité, s'appelle la *transformée de Fourier discrète*. Son calcul rapide est effectué par un algorithme de type « diviser pour régner ».

Pour appliquer cette idée, supposons que $n$ est pair, $n = 2k$. Alors, $\omega^k = -1$ puisque

$$ (\omega^k - 1)(\omega^k + 1) = \omega^{n-1} = 0 $$

et le premier facteur n'est pas diviseur de 0. Le polynôme $F$ est décomposé par division euclidienne de deux façons :

$$ F = Q_0(X^k - 1) + R_0 \quad \text{et} \quad F = Q_1(X^k + 1) + R_1, $$

avec $\deg R_0 < k$ et $\deg R_1 < k$. Ces décompositions vont nous permettre le calcul de $F$ sur les puissances paires et impaires de $\omega$. En effet, si $\ell$ est pair, $\omega^{k\ell} = 1$ et donc $F(\omega^\ell) = R_0(\omega^\ell)$. De même, si $\ell$ est impair, $F(\omega^\ell) = R_1(\omega^\ell)$. Outre l'application
**Entrée** F = f₀ + ... + fₙ₋₁Xⁿ₋₁; les puissances 1, ω, ..., ωⁿ₋₁ d'une racine n-ième principale de l'unité ω, n étant une puissance de 2.

**Sortie** F(1), ..., F(ωⁿ₋₁).

1. Si n = 1, renvoyer f₀.

2. Sinon, soit k = n/2. Calculer

$$R_0(X) = \sum_{j=0}^{k-1} (f_j + f_{j+k})X^j,$$

$$\bar{R}_1(X) = R_1(\omega X) = \sum_{j=0}^{k-1} (f_j - f_{j+k})\omega^j X^j.$$

3. Calculer récursivement R₀(1), R₀(ω²), ..., R₀((ω²)ᵏ₋₁) et R̄₁(1), R̄₁(ω²), ..., R̄₁((ω²)ᵏ₋₁).

4. Renvoyer

$$R_0(1), \bar{R}_1(1), R_0(\omega^2), \bar{R}_1(\omega^2), \dots, R_0((\omega^2)^{k-1}), \bar{R}_1((\omega^2)^{k-1}).$$

**Algorithme 2.3 – Transformée de Fourier rapide (FFT).**

récursive, le point crucial qui est la source de l'efficacité de l'algorithme de transformée de Fourier rapide (Algorithme 2.3) et qui conduit au choix de racines primitives de l'unité, est que le calcul de R₀ et R₁ est très simple (étape (2)). Lors des appels récursifs, les puissances de ω qui sont utilisées sont des ω²ᵢ, qui sont bien des racines primitives d'après le Lemme 2.5.

**Théorème 2.6** L'Algorithme 2.3 de transformée de Fourier rapide requiert au plus $\frac{3n}{2} \log n$ opérations dans $\mathbb{A}$; les multiplications font toutes intervenir une puissance de ω.

*Démonstration.* Puisque les puissances de ω sont connues, le coût de l'appel en degré n est d'au plus 2 × n/2 additions et soustractions (pour le calcul de R₀ et R₁) et de n/2 multiplications (pour le calcul de $\bar{R}_1$), plus deux appels récursifs en degré n/2. Sa complexité F(n) satisfait donc à la récurrence :

$$F(n) \leq \frac{3n}{2} + 2F\left(\frac{n}{2}\right)$$

et le lemme « diviser pour régner » (avec $p = q = m = s = 2$, T(n) = 3n/2 et κ = 0) permet de conclure. ■

**Exercice 2.8** Montrer que l'Algorithme 2.3 requiert n log n additions dans $\mathbb{A}$, $\frac{1}{2}n \log n$ multiplications d'éléments de $\mathbb{A}$ par des puissances de ω, mais aucune autre multipli-
cation dans $\mathbb{A}$.

R La transformée de Fourier discrète est un morphisme d'algèbres sur $\mathbb{A}$ de $\mathbb{A}[X]/(X^n - 1)$ dans $\mathbb{A}^n$ avec comme multiplication dans $\mathbb{A}^n$ la multiplication coordonnée par coordonnée. Cette observation permet d'économiser des transformées inverses en effectuant plusieurs calculs directement sur les transformées. Une application typique de cette observation est le produit scalaire, ou plus généralement le produit de matrices.

### Interpolation

En termes matriciels, l'opération $F \mapsto \text{Ev}(F)$ est linéaire et sa matrice (pour des polynômes $F$ de degré au plus $n-1$, dans la base monomiale $\{1, X, \dots, X^{n-1}\}$) est la matrice de Vandermonde

$$V_{\omega} = \begin{pmatrix} 1 & 1 & \cdots & 1 \\ 1 & \omega & \cdots & \omega^{n-1} \\ \vdots & & & \vdots \\ 1 & \omega^{n-1} & \cdots & \omega^{(n-1)^2} \end{pmatrix}$$

**Lemme 2.7** Si $\omega \in \mathbb{A}$ est racine $n$-ième principale de l'unité, alors $V_{\omega^{-1}} V_{\omega} = n I_n$.

*Démonstration.* D'après le Lemme 2.5, $\omega$ est inversible. La matrice $V_{\omega^{-1}}$ est donc bien définie. Le coefficient de la $i$-ième ligne et $j$-ième colonne du produit s'écrit

$$\sum_{k=0}^{n-1} \omega^{(i-1)k} \omega^{-(j-1)k} = \sum_{k=0}^{n-1} \omega^{(i-j)k}.$$

D'après le Lemme 2.5, cette dernière somme est nulle pour $0 < i - j < n$. Par symétrie des matrices, elle est aussi nulle si $0 < j - i < n$. Les coefficients en dehors de la diagonale sont donc tous nuls. Sur la diagonale, la somme comporte $n$ termes, tous égaux à 1, d'où le résultat.

Autrement dit, l'interpolation sur les puissances de $\omega$ est calculée efficacement en la ramenant à une FFT sur les puissances de $\omega^{-1}$, qui est bien principale d'après le Lemme 2.5.

### Conclusion

Les résultats de cette section sont résumés dans le théorème suivant.

**Théorème 2.8** Si $2$ est inversible dans $\mathbb{A}$ et $n$ une puissance de $2$, étant donnée une racine principale $n$-ième dans $\mathbb{A}$, le produit de deux polynômes dont la somme des degrés est inférieure à $n$ peut être calculé en $\frac{9}{2}n \log n + O(n)$ opérations dans $\mathbb{A}$. Seuls $n$ des produits sont entre deux éléments de $\mathbb{A}$ qui ne sont pas des puissances de $\omega$.
**Démonstration.** La complexité de l'Algorithme 2.2 est de 3 FFT en degré $n$, soit $\frac{9}{2}n \log n$ opérations, plus $O(n)$ divisions par $n$ pour l'interpolation (ce qui est possible puisque 2 est inversible) et $O(n)$ multiplications pour calculer les puissances de $\omega$. Les multiplications de la FFT font intervenir des puissances de $\omega$, les autres sont celles de l'étape de produit point à point, au nombre de $n$. ■

**Exercice 2.9** Montrer que sous les hypothèses du théorème précédent, on peut multiplier ces polynômes en utilisant $3n \log n + O(n)$ additions dans $\mathbb{A}$, plus $\frac{3}{2}n \log n + O(n)$ multiplications par des puissances de $\omega$, $n$ multiplications arbitraires dans $\mathbb{A}$, et $n$ divisions par $n$. ■

**Exercice 2.10** Soit $n$ dans $\mathbb{N}$, soit $n_0$ la plus petite puissance de 2 supérieure ou égale à $n$, et supposons qu'il existe dans $\mathbb{A}$ une racine $n_0$-ième principale de l'unité. Cette racine étant connue, on peut multiplier les polynômes de degré au plus $\frac{1}{2}n - 1$ en $9n \log n + O(n)$ opérations dans $\mathbb{A}$. ■

**Exercice 2.11** Soit $n$ une puissance de 2, et supposons qu'on dispose d'une racine $n$-ième principale de l'unité $\omega \in \mathbb{A}$. Soit $P$ et $Q$ deux polynômes dans $\mathbb{A}[X]$ de degré au plus $n - 1$. Supposons que les coefficients de $X^0, X^1, \dots, X^{n-1}$ du produit $R = PQ$ sont connus. Montrer que $R$ peut être calculé en $\frac{9}{2}n \log n + O(n)$ opérations dans $\mathbb{A}$. ■

## Mais où sont nos racines ?

Même dans le cas favorable où $\mathbb{A}$ est un corps, il n'y existe pas nécessairement toutes les racines primitives de l'unité nécessaires pour appliquer les techniques des algorithmes précédents. Dans le cas particulier des corps finis, on sait donner une réponse précise à cette question d'existence.

**Proposition 2.9** Le corps fini à $q$ éléments $\mathbb{F}_q$ contient une racine $n$-ième primitive de l'unité si et seulement si $n$ divise $q-1$.

**Exercice 2.12** 1. Prouver la proposition précédente.

2. Montrer que si $n$ divise $q-1$ et si $\alpha$ est un élément primitif de $\mathbb{F}_q$ (c'est-à-dire tel que $\alpha$ engendre le groupe multiplicatif $(\mathbb{F}_q \setminus \{0\}, \times))$ alors $\alpha^{(q-1)/n}$ est une racine $n$-ième primitive de l'unité. ■

Pour pouvoir appliquer le Théorème 2.6 avec $n$ assez grand et $\mathbb{A} = \mathbb{Z}/p\mathbb{Z}$ (c'est-à-dire $q = p$ premier), ce résultat mène à la notion de premiers de Fourier, qui sont les nombres premiers $p$ tels que $p-1$ soit divisible par une grande puissance de 2, donc de la forme $\ell 2^e + 1$, avec $e$ « suffisamment grand » (qui est appelé l'exposant de $p$). Par exemple,

$$4179340454199820289 = 29 \times 2^{57} + 1$$

est un tel nombre premier. Ainsi, dans $\mathbb{Z}/4179340454199820289\mathbb{Z}$, on dispose de racines primitives $2^{57}$-ième de l'unité (21 en est une); on peut donc y multiplier des polynômes de degrés colossaux par cet algorithme en $O(n \log n)$ opérations.
Ces premiers sont donc très utiles pour la FFT sur des corps finis. On peut montrer qu'il y a approximativement $(x/\log x)/2^{e-1}$ premiers de Fourier d'exposant $e$ inférieurs à $x$ (à rapprocher du théorème des nombres premiers qui montre l'existence d'approximativement $x/\log x$ premiers inférieurs à $x$). Il s'ensuit qu'il y a environ 130 corps finis de la forme $k = \mathbb{F}_p$ tel que $p$ a la taille d'un (demi-)mot machine (32 bits) et tel que dans $k$ puisse multiplier par FFT des polynômes de degré de l'ordre d'un million. Les racines de l'unité de $\mathbb{F}_p$ peuvent être calculées à partir des éléments primitifs. Si on choisit au hasard un élément de $\mathbb{F}_p$, lorsque $p \to \infty$, la probabilité qu'il soit primitif tend vers $6/\pi^2$, soit plus de 0,6.

## 2.5 L'algorithme de Schönhage-Strassen

Quand les racines de l'unité font défaut, il reste possible de faire fonctionner les idées à base de transformée de Fourier. Ceci est réalisé par l'algorithme de Schönhage et Strassen. Cet algorithme s'applique quel que soit l'anneau de base, pourvu que $2$ y soit inversible; l'idée est de rajouter les racines de l'unité qui manquent en étendant l'anneau de base de manière judicieuse.

### Racines virtuelles de l'unité

Le point de départ de l'algorithme est résumé dans le résultat suivant.

**Lemme 2.10** Si $2$ est inversible dans $\mathbb{A}$ et $n$ est une puissance de $2$, alors $\omega = X \bmod (X^n + 1)$ est une racine $2n$-ième principale de l'unité dans $\mathbb{A}[X]/(X^n + 1)$.

*Démonstration.* Comme $\omega^n = -1$, $\omega$ est bien une racine $2n$-ième de l'unité. Il reste à prouver qu'elle est principale. D'après la remarque en page 42, il s'agit de montrer que $\omega^t - 1$ n'est pas diviseur de 0, quel que soit le diviseur strict $t$ de $2n$. Puisque $n$ est une puissance de $2$, l'entier $t$ divise $n$, et donc $\omega^t - 1$ divise l'élément $\omega^n - 1 = -2$ qui est supposé inversible, ce qui conclut la preuve. ■

**Exemple 2.3** Même lorsque $\mathbb{A}$ est un corps, ces anneaux ne sont en général pas intègres. Par exemple, dans $\mathbb{Z}/3\mathbb{Z}$, la factorisation $Y^4 + 1 = (Y^2 + Y + 2)(Y^2 + 2Y + 2)$ exhibe deux polynômes diviseurs de 0 dans $\mathbb{Z}/3\mathbb{Z}[Y]/(Y^4 + 1)$.

### L'algorithme

La méthode est résumée dans l'Algorithme 2.4, qui calcule le produit FG modulo $X^n + 1$: ceci permet une utilisation récursive mais ajoute une petite complication pour le produit par FFT, qui calcule un produit modulo $X^n + 1$, d'où les changements de variables et choix de racines de l'unité dans l'étape 2.

**Exemple 2.4** Un exemple jouet permet de mieux comprendre comment se déroule cet algorithme et pourquoi il est correct. La multiplication de polynômes de $\mathbb{Z}/3\mathbb{Z}[X]$ par FFT ne peut pas aller au-delà du degré 2, puisqu'il n'y a pas de racine
**Entrée** F et G de degrés strictement inférieurs à n, où $n = 2^k$, $k > 2$.

**Sortie** FG mod $X^n + 1$.

1. Soient $d = 2^{\lfloor k/2 \rfloor}$ et $\delta = n/d$. Récrire F et G sous la forme

$$ \begin{aligned} \bar{F}(X, Y) &= F_0(X) + F_1(X)Y + \cdots + F_{\delta-1}(X)Y^{\delta-1}, \\ \bar{G}(X, Y) &= G_0(X) + G_1(X)Y + \cdots + G_{\delta-1}(X)Y^{\delta-1}, \end{aligned} $$

de sorte que les $F_i$ et $G_i$ aient des degrés strictement inférieurs à d et que $F(X) = \bar{F}(X, X^d)$ et $G(X) = \bar{G}(X, X^d)$.

2. Calculer $\bar{H} := \bar{F}\bar{G}$ mod $Y^\delta + 1$ dans $\mathbb{B}[Y]$ par FFT, où $\mathbb{B} = A[X]/(X^{2d} + 1)$ et les produits dans $\mathbb{B}$ sont effectués récursivement. Pour cela,

— si $k$ est pair, calculer $\bar{H}(X, X^2Y)$ mod $Y^\delta - 1$ avec $X^4$ comme racine de l'unité, puis reconstruire $\bar{H}$ avec $X^{-2} = -X^{2d-2}$;

— si $k$ est impair, calculer $\bar{H}(X, XY)$ mod $Y^\delta - 1$ avec $X^2$ comme racine de l'unité, puis reconstruire $\bar{H}$ avec $X^{-1} = -X^{2d-1}$.

3. Renvoyer $H(X, X^d)$.

**Algorithme 2.4** – Algorithme de Schönhage-Strassen sur les polynômes.

$n$-ième de l'unité d'ordre plus grand. Pour calculer le produit de

$$ F = 1 + X^7 \quad \text{et} \quad G = 1 + X + X^6 \quad \text{modulo} \quad X^{16} + 1 $$

dans $\mathbb{Z}/3\mathbb{Z}[X]$, l'algorithme introduit $d = 4$, $\delta = 4$ et

$$ \bar{F} = 1 + X^3Y, \quad \bar{G} = 1 + X + X^2Y, $$

et cherche à calculer le produit de ces deux polynômes

$$ \bar{H}(X, Y) = \bar{F}(X, Y)\bar{G}(X, Y) \mod Y^4 + 1, $$

vus comme polynômes en $Y$ à coefficients dans $\mathbb{Z}/3\mathbb{Z}[X]/(X^8 + 1)$. Ce sera suffisant pour connaître le résultat puisqu'en $Y = X^4$, on a $Y^4 + 1 = X^{16} + 1$, polynôme modulo lequel l'algorithme calcule le produit.

Pour se ramener à la situation d'une FFT, il suffit de considérer $\bar{H}(X, X^2Y)$ : le remplacement de $Y$ par $X^2Y$ dans l'identité ci-dessus montre qu'on calcule alors modulo $X^8Y^4 + 1 = -Y^4 + 1$. On commence donc par calculer une DFT de

$$ \bar{F}(X, X^2Y) = 1 + X^5Y, \quad \bar{G}(X, X^2Y) = 1 + X + X^4Y $$
avec $X^4$ comme racine 4e de l'unité. On obtient les évaluations en 1, $X^4$, $X^8 = -1$, $X^{12} = -X^4$, à savoir $(1+X^5, 1-X, 1-X^5, 1+X)$ et $(1+X+X^4, X, 1+X-X^4, 2+X)$ qui se calculent par une FFT très peu chère, puisque les multiplications par des puissances de X ne demandent pas d'opération arithmétique.

Ensuite, il faut multiplier deux à deux ces valeurs dans $\mathbb{Z}/3\mathbb{Z}[X]/(X^8+1)$ pour obtenir les évaluations de $\bar{H}$. L'algorithme procède récursivement parce que c'est encore possible (lorsque $n$ devient trop petit, $2d = n$ ce qui mènerait à une boucle infinie; on a alors recours à une autre méthode, par exemple un produit naïf en ces petits degrés.) Ces calculs récursifs mènent aux valeurs suivantes pour $\bar{H}(X, X^2Y)$ en 1, $X^4$, $X^8 = -1$, $X^{12} = -X^4$ :

$$c_0 = 1 + X^4 + X^5 + X^6, c_1 = X - X^2, c_2 = 1 - X^4 - X^5 - X^6, c_3 = 2 + X^2.$$

L'interpolation par FFT consiste à évaluer le polynôme $c_0 + c_1 Y + c_2 Y^2 + c_3 Y^3$ en les puissances de $X^{-4} = -X^4$, à savoir $1, -X^4, -1, X^4$, ce qui donne

$$1 + X, X^4 + X^5 + X^6, 2X, 0$$

et donc finalement (après multiplication par $1/n = 1/4 = 1$)

$$\bar{H}(X, X^2Y) = 1 + X + X^4(1 + X + X^2)Y + 2XY^2.$$

En utilisant $X^{-2} = -X^6$ pour remplacer $Y$ par $YX^{-2}$, il vient

$$\bar{H}(X, Y) = 1 + X + X^2(1 + X + X^2)Y + X^5Y^2.$$

Le résultat s'obtient enfin en évaluant en $Y = X^4$ :

$$FG = 1 + X + X^6(1 + X + X^2) + X^{13} + X^{12} = 1 + X + X^6 + X^7 + X^8 + X^{13}.$$

En ce qui concerne la complexité, la première étape ne demande pas d'opération arithmétique, et l'évaluation de la dernière étape se fait en complexité linéaire, ainsi que les changements de $Y$ en $X^kY$ avant et après la FFT. Ensuite, d'après le Théorème 2.8, la multiplication $\bar{FG}$ nécessite :

- $O(d \log d)$ opérations $(+, -)$ dans $\mathbb{B}$, chacune en $O(d)$ opérations dans $\mathbb{A}$;

- $O(d \log d)$ multiplications par une puissance de $\omega$ dans $\mathbb{B}$, chacune en $O(d)$ opérations dans $\mathbb{A}$, par simples décalages des indices et éventuels changements de signes;

- $O(d)$ divisions par $d$ dans $\mathbb{B}$, chacune en $O(d)$ opérations dans $\mathbb{A}$;

- au plus $\delta$ produits dans $\mathbb{B}$.

Le coût $C(n)$ de cet algorithme obéit donc à l'inégalité

$$C(n) \leq K d \log d \times d + \delta C(2d)$$

où $K$ est une constante indépendante de $n$ et de $\mathbb{A}$. En divisant par $n = d\delta$, cette inégalité devient

$$\frac{C(n)}{n} \leq K \frac{d}{\delta} \log d + 2 \frac{C(2d)}{2d} \leq K' \log n + 2 \frac{C(2d)}{2d}$$
pour $n$ assez grand, et une constante $K'$ indépendante de $n$ et de $\mathbb{A}$. Ici $2d = 2^{\lfloor k/2 \rfloor+1}$, et pour se ramener au théorème diviser pour régner, il sera plus commode de faire apparaître $\lfloor(k+1)/2\rfloor = \lceil k/2 \rceil$. Pour cela, on considère donc d'abord

$$ \frac{C(2n)}{2n} \le K' \log n + 2 \frac{C(2^{\lceil k/2 \rceil+1})}{2^{\lceil k/2 \rceil+1}} $$

et on pose $c(k) = C(2^{k+1})/2^{k+1}$, ce qui donne par récriture

$$ c(k) \le K''k + 2c(\lceil k/2 \rceil), $$

pour une nouvelle constante $K''$. Le résultat final s'en déduit par application du théorème « diviser pour régner ».

**Théorème 2.11** Soit $\mathbb{A}$ un anneau dans lequel 2 est inversible (d'inverse connu). Alors, deux polynômes de $\mathbb{A}[X]$ de degrés au plus $n$ peuvent être multipliés en $O(n \log n \log \log n)$ opérations $(+, -, \times)$ dans $\mathbb{A}$.

Il est possible d'étendre cette idée au cas où 3 est inversible (FFT triadique), et plus généralement à un anneau quelconque. On obtient alors l'algorithme de complexité $O(n \log n \log \log n)$ mentionné dans l'introduction.

## 2.6 Algorithmes pour les entiers

Les algorithmes ainsi que les résultats présentés ci-dessus s'étendent à la multipli-cation des entiers<fnref n="1" />. Nous allons brièvement présenter cette problématique à travers un exemple.

Soient à multiplier les entiers 2087271 et 1721967, qu'on suppose donnés en base 2,

$$ A = 111111101100101100111 \quad \text{et} \quad B = 110100100011001101111, $$

chacun ayant $D = 21$ chiffres binaires. On peut ramener leur multiplication à un produit de polynômes. Plus exactement, on associe à A et B les polynômes de degré strictement inférieur à D :

$$ P = X^{20} + X^{19} + X^{18} + X^{17} + X^{16} + X^{15} + X^{14} + X^{12} + X^{11} + X^{8} + X^{6} + X^{5} + X^{2} + X + 1 $$

et

$$ Q = X^{20} + X^{19} + X^{17} + X^{14} + X^{10} + X^{9} + X^{6} + X^{5} + X^{3} + X^{2} + X + 1. $$

La stratégie est la suivante : on calcule $R = PQ$ dans $\mathbb{Z}[X]$ puis on évalue $R$ en $X = 2$. Pour multiplier $P$ et $Q$ dans $\mathbb{Z}[X]$, il suffit d'effectuer leur produit dans $\mathbb{A}[X] = (\mathbb{Z}/p\mathbb{Z})[X]$, où $p$ est un nombre premier tel que $2D > p > D$. Par le Théorème 2.11, cette multiplication polynomiale en degré $D$ peut se faire en $O(D \log D \log \log D)$ opérations $(+, -, \times)$ dans $\mathbb{A} = \mathbb{Z}/p\mathbb{Z}$.

<footnote n="1">1. Historiquement, les algorithmes pour la multiplication des entiers ont été introduits *avant* leurs homologues polynomiaux, alors que ces derniers sont souvent bien plus simples à énoncer.</footnote>
Puisque chaque opération de $\mathbb{A}$ nécessite au plus $O(\log^2 D)$ opérations binaires, il s'ensuit que le coût binaire du calcul de $R$ est de $O(D\log^3 D \log \log D)$. Ici $p = 23$ et $R$ (écrit en base 2) vaut $R = X^{40} + 10X^{39} + 10X^{38} + 11X^{37} + 11X^{36} + 11X^{35} + 100X^{34} + 11X^{33} + 11X^{32} + 100X^{31} + 11X^{30} + 100X^{29} + 101X^{28} + 11X^{27} + 101X^{26} + 1000X^{25} + 101X^{24} + 101X^{23} + 1000X^{22} + 1001X^{21} + 1010X^{20} + 1000X^{19} + 111X^{18} + 1000X^{17} + 110X^{16} + 110X^{15} + 110X^{14} + 11X^{13} + 100X^{12} + 110X^{11} + 100X^{10} + 11X^{9} + 100X^{8} + 100X^{7} + 100X^{6} + 11X^{5} + 10X^{4} + 11X^{3} + 11X^{2} + 10X + 1$. Enfin, l'évaluation de $R$ en 2 revient à gérer les retenues et cela peut se faire en un coût négligeable. Ici $AB = R(2) = 110100010011010111101101110110110110101001$, ou encore, en écriture décimale $AB = 3594211782057$.

Une légère amélioration est possible si l'on choisit pour $p$ un premier de Fourier supérieur à 2D. Dans notre exemple, on peut prendre $p = 193 = 3 \cdot 2^6 + 1$. En effet, il existe dans $\mathbb{A} = \mathbb{Z}/193\mathbb{Z}$ une racine primitive $\omega = 125$, d'ordre $2^6 = 64$, donc supérieur à $2D = 40$. Ce choix mène à un algorithme de complexité binaire $O(D\log^3 D)$.

Une autre approche est de calculer PQ dans $\mathbb{Z}[X]$ par DFT en calculant avec des nombres complexes représentés par des approximations numériques flottantes. Une estimation des erreurs numériques montre qu'il suffit de calculer en précision fixe $O(\log^2 D)$. La complexité binaire est donc toujours $O(D\log^3 D)$ via cette approche.

On peut faire un peu mieux, en remplaçant la base 2 par une base B telle que $\text{BlogB}$ soit de l'ordre de $\log D$ et en appliquant l'un des raisonnements précédents; on peut aboutir ainsi à un algorithme de complexité binaire $O(D\log^2 D)$. Qui plus est, en appelant récursivement l'algorithme pour multiplier les petits entiers, on peut descendre cette complexité à $O(D \log D \log \log D \log \log \log D \dots)$.

Une meilleure solution est cependant détenue par la version entière, à base de FFT dans des anneaux de type $\mathbb{Z}/(2^{2^k} + 1)\mathbb{Z}$, de l'algorithme de Schönhage-Strassen, dont nous nous contentons de mentionner l'existence.

**Théorème 2.12** On peut multiplier deux entiers de $D$ chiffres binaires en utilisant $O(D \log D \log \log D)$ opérations binaires.

En développant largement ces techniques de FFT, la complexité du produit entier peut être réduite à $O(n \log n \cdot 2^{O(\log^* n)})$, où $\log^* n$ est le nombre d'itérations du logarithme pour ramener $n$ en dessous de 1. Par exemple (pour la base 2), $\log^* 256 = 4$, car $\log 256 = 8$, $\log 8 = 3$, $1 < \log 3 < 2$ et enfin $\log \log 3 < 1$.

## 2.7 Un concept important : les fonctions de multiplication

Un bon nombre des résultats de complexité donnés dans la suite de cet ouvrage reposent sur la notion de *fonction de multiplication*. Une fonction $M: \mathbb{N} \rightarrow \mathbb{N}$ sera dite fonction de multiplication polynomiale (pour un anneau $\mathbb{A}$) si :

- on peut multiplier les polynômes de $\mathbb{A}[X]$ de degrs au plus $n$ en au plus $M(n)$ opérations dans $\mathbb{A}$;

- $M$ vérifie l'inégalité $M(n+n') \geq M(n)+M(n')$.

Ainsi, on sait d'après ce qui précède que des fonctions de la forme $n^{\log_2 3}$ ou $n \log n$ sont (à des constantes près) des fonctions de multiplication (respectivement pour tous les anneaux ou ceux qui permettent la transformée de Fourier). L'intérêt de cette
notion est qu'elle permet d'énoncer des résultats de complexité indépendants du choix
de l'algorithme utilisé pour multiplier les polynômes (même si parfois cela mène à
des estimations légèrement pessimistes, parce que les algorithmes n'exploitent pas la
remarque page 44).

La seconde condition est utilisée pour écarter l'hypothèse d'une fonction qui croît trop lentement (si M est constante, elle ne vérifie pas cette condition). Concrètement, elle permettra par exemple d'établir que dans des algorithmes du type « inversion de série formelle par l'opérateur de Newton » (Chapitre 3), le coût total est essentiellement celui de la dernière étape.

De la même manière, on est amené à introduire la notion de fonction de multiplication pour les entiers. Une fonction $M_Z : N \to N$ est une fonction de multiplication pour les entiers si :

— on peut multiplier des entiers de *n* chiffres en $M_Z(n)$ opérations binaires;

— $M_Z$ vérifie l'inégalité $M_Z(n + n') \ge M_Z(n) + M_Z(n')$.

Ce concept est très proche de son analogue polynomial et les contextes d'utilisation
sont souvent les mêmes : on utilise la seconde condition pour contrôler les coûts de
multiplication lors de l'analyse d'algorithmes récursifs.

Exercices

**Exercice 2.13** Soit $\mathbb{A}$ un anneau, soit $a \in \mathbb{A}$ et soit $P$ un polynôme de $\mathbb{A}[X]$, de degré au plus $n$. On se propose de calculer le polynôme $Q(X) = P(X + a)$.

1. Donner un algorithme direct pour le calcul de Q et estimer sa complexité en termes de n;

2. Montrer qu'il est possible de calculer $Q(X)$ en utilisant $O(M(n) \log n)$ opérations $(+, -, \times)$ dans $\mathbb{A}$. (Un algorithme en $O(M(n))$ opérations sera donné au Chapitre 3 par somme composée, et une autre solution en $O(M(n))$ est donnée en Exercice 33.2 en Annexe.)

**Exercice 2.14** Soit $\mathbb{A}$ un anneau, soient $\kappa, n \in \mathbb{N}$ et soit $P$ un polynôme de $\mathbb{A}[X]$, de degré au plus $n$.

1. Donner un algorithme direct pour le calcul de $P(X)^\kappa$ et estimer sa complexité en fonction de $\kappa$ et de $n$;

2. Montrer qu'il est possible de calculer $P(X)^\kappa$ en n'utilisant que $O(M(n\kappa))$ opérations $(+, -, \times)$ dans $\mathbb{A}$.

**Exercice 2.15** — **Produit court.** Pour deux polynômes $A, B \in \mathbb{Q}[X]$ de degré strictement inférieur à $n$, le polynôme $AB$ mod $X^n$ est appelé *produit court* de $A$ et $B$. Le but de cet exercice est de montrer qu'il est possible de calculer plus efficacement le produit court que le produit entier AB. Le gain est d'un facteur constant.

1. Prouver que tel est le cas avec la multiplication polynomiale naïve.

On suppose que l'on dispose d'un algorithme de multiplication polynomiale qui, en degré $n$, utilise $M_\alpha(n) = K n^\alpha$ opérations dans $\mathbb{Q}$, avec $\alpha > 1$ et $K > 0$.
2. Donner un algorithme de type « diviser pour régner » pour le calcul du produit court de A, B.

3. Montrer que la complexité de cet algorithme est $\frac{1}{2^{\alpha-2}} M_{\alpha}(n) + O(n \log n)$. En déduire l'impact sur la multiplication naïve, et celle de Karatsuba.

4. Modifier la conception de l'algorithme précédent, en découpant en deux sous-problèmes de tailles $\beta n$ et $(1-\beta)n$, avec $1/2 \le \beta \le 1$ laissé en paramètre. Prouver que le nouvel algorithme utilise

$$ \frac{\beta^\alpha}{1 - 2(1 - \beta)^\alpha} \cdot M_\alpha(n) + O(n \log n) $$

opérations dans $\mathbb{Q}$.

5. Montrer que le minimum de la fonction $\beta \mapsto \frac{\beta^\alpha}{1-2(1-\beta)^\alpha}$ est atteint en $\beta_{\min} = 1 - \left(\frac{1}{2}\right)^{\frac{1}{\alpha-1}}$ et en déduire qu'on peut calculer le produit court en $C(\alpha)M_\alpha(n) + O(n \log n)$ opérations dans $\mathbb{Q}$, où

$$ C(\alpha) = \frac{\left(2^{\frac{1}{\alpha-1}} - 1\right)^{\alpha}}{2^{\frac{\alpha}{\alpha-1}} - 2} $$

6. Conclure sur le gain obtenu (à l'aide d'un calcul numérique).

**Exercice 2.16** Soient $P_1, \dots, P_t$ des polynômes de $\mathbb{A}[X]$, de degrés $d_1, \dots, d_t$. Montrer que le produit $P_1 \cdots P_t$ peut s'effectuer en $O(M(n) \log t)$ opérations dans $\mathbb{A}$, où $n = \sum_i d_i$.

### Notes

Le mathématicien russe A. N. Kolmogorov avait conjecturé au début des années 1960 qu'il serait impossible de multiplier deux entiers en un coût binaire sous-quadratique. En 1962 cette conjecture fut infirmée par Karatsuba et Ofman [KO63]. Une généralisation de leur algorithme a été proposée peu de temps après par Toom et Cook [Coo66; Too63]. L'algorithme de Toom-Cook a une complexité binaire de $O(n \cdot 32^{\sqrt{\log_2 n}})$ pour multiplier deux entiers de taille binaire $n$; ce résultat peut être importé dans le monde polynomial (voir l'Exercice 2.7 pour une version plus faible). L'algorithme de FFT a une longue histoire, qui remonte à Gauss [Coo90; CT93; HJB85]. Il s'agit d'un progrès algorithmique très célèbre : Dongarra et Sullivan [DS00] le placent parmi les dix algorithmes qui ont le plus marqué le développement des sciences de l'ingénieur du 20e siècle. L'article fondateur de Cooley et Tukey [CT65] est l'un des plus cités en informatique<fnref n="2" />. La variante de la DFT décrite dans ce chapitre, appelée *decimation-in-frequency* en anglais, est due à Gentleman et Sande [GS66]; il s'agit d'une version duale (au sens du principe de transposition de Tellegen évoqué au Chapitre 12) de l'algorithme *decimation-in-time* [CT65]. De nombreux livres [Bri74; Nus81; Van92] et articles [Bera; DV90; FJ05] permettent de se familiariser avec la

<footnote n="2">2. Plus de 2000 citations, d'après la base bibliographique *Science Citation Index (SCI®)*.</footnote>
myriade de techniques de type FFT; en particulier, Frigo et Johnson [FJ05] décrivent l'une des implantations les plus efficaces de la transformée de Fourier complexe (appelée FFTW, pour *Fastest Fourier Transform in the West*.)

Une avancée importante a été la découverte de Schönhage et Strassen [SS71] du résultat équivalent pour la multiplication des nombres entiers. Pendant longtemps, on a cru que cet algorithme ne pourrait présenter qu'un intérêt purement théorique.

À ce jour, la plupart des logiciels généralistes de calcul formel disposent d'une multiPLICATION rapide d'entiers qui inclut la FFT. MAPLE, MATHEMATICA et SAGE utilisent la bibliothèque GMP. Il faut noter que cette bibliothèque est le résultat d'un travail de nombreuses années, qui comporte une partie importante de code assembleur consacré à la multiplication sur chacun des processeurs produits dans une période récente.

L'analogue polynomial de l'algorithme de Schönhage-Strassen traité en Section 2.5 a été suggéré par Nussbaumer [Nus80]. Le cas particulier de la caractéristique 2 est traité par Schönhage [Sch77]. Plus récemment, Cantor et Kaltofen [CK91] ont donné un algorithme qui multiplie des polynômes de degré $n$ sur une algèbre $\mathbb{A}$ (non nécessairement commutative, ni associative) en $O(n \log n \log \log n)$ opérations dans $\mathbb{A}$.

En ce qui concerne la terminologie, nous suivons Demazure [Dem97] pour la définition des racines principales de l'unité. D'autres auteurs [Für09] utilisent la troisième propriété du Lemme 2.5. D'autres enfin [GG03] les appellent simplement primitives.

La borne de complexité binaire $O(n \log n \cdot 2^{O(\log n^*)})$ mentionnée en fin de Section 2.6, pour multiplier des entiers de $n$ chiffres, a été obtenue récemment par Fürer [De+08; Für07; Für09]. Harvey, van der Hoeven et Lecerf ont ensuite montré que l'algorithme de Fürer pouvait être un peu modifié pour atteindre un coût binaire en $O(n \log n \cdot 16^{\log_2 n})$ [HHL16]. Puis ils ont conçu un autre algorithme ayant un coût $O(n \log n \cdot 8^{\log_2 n})$. À ce jour, on ne connaît pas de borne analogue pour la complexité arithmétique de la multiplication polynomiale en degré $n$ en toute généralité. Néanmoins si les coefficients sont dans le corps fini $\mathbb{F}_q$ alors la complexité binaire du produit en degré $n$ est en $O(n \log q \log(n \log q) \cdot 8^{\log_2(n \log q)})$ [HHL14].

Un problème ouvert est de trouver un algorithme de multiplication polynomiale en complexité $O(n)$, ou de prouver qu'un tel algorithme n'existe pas. Morgenstern [Mor73] a donné une première réponse partielle, en montrant que dans un modèle simplifié où $\mathbb{A} = \mathbb{C}$ et où les seules opérations admises sont les additions et les multiplications par des nombres complexes de module inférieur à 1, au moins $\frac{1}{2}n \log n$ opérations sont nécessaires pour calculer une DFT en taille $n$. Plus récemment, Bürgisser et Lotz [BL04] ont étendu ce type de borne à la multiplication dans $\mathbb{C}[X]$. Concernant les entiers, Cook et Aanderaa [CA69] ont prouvé une borne inférieure de la forme $cn \log n/(\log \log n)^2$ pour la multiplication en taille binaire $n$ sur une machine de Turing. Cette borne a été améliorée à $O(n \log n)$ pour un modèle restreint (multiplication en ligne) [PFM74].

Malgré la simplicité de leur idée de base, les algorithmes de Karatsuba et de Toom-Cook soulèvent encore des questions délicates; plusieurs travaux récents traitent diverses généralisations (degrés déséquilibrés ou non-puissances, caractéristique non nulle de l'anneau des scalaires) [BZ07; CH07; Mon05; WP06]. Par exemple, connaître le nombre minimum de multiplications $sm(n)$ (resp. $sm_p(n)$) suffisantes pour multi-
plier des polynômes de degré donné $n$ sur un anneau arbitraire (resp. sur un anneau de caractéristique $p$) est un problème important et difficile. Pour les premières valeurs de $n$, les meilleures bornes connues actuellement, dues à Montgomery [Mon05], sont $sm(1) \le 3$, $sm(2) \le 6$, $sm(3) \le 9$, $sm(4) \le 13$, $sm(5) \le 17$ et $sm(6) \le 22$. L'optimálité de ces bornes est un problème difficile. Pour $sm_p(n)$, on sait par exemple que $sm_2(1) = 3$, $sm_2(2) = 6$, $sm_2(3) = 9$ [Kam85] et que $sm_7(4) = 10$, $sm_5(4) \le 11$ [CKO09].

Une autre question importante est de savoir si les $n$ coefficients du produit court FG mod $X^n$ de deux polynômes F et G de degrés au plus $n$ peuvent être calculés plus efficacement que l'ensemble des $2n$ coefficients du produit FG. L'Exercice 2.15, inspiré d'un résultat de Mulders [Mul00], apporte une réponse positive dans le cas où l'algorithmme de multiplication est celui de Karatsuba. Par contre, on ne connaît aucun élément de réponse dans le cas où la multiplication polynomiale repose sur la DFT. Une question liée, le calcul du produit médian [HQZ04], sera abordée au Chapitre 12.

Une faiblesse de l'algorithmme DFT est son comportement quasi-linéaire par morceaux, dû aux sauts de complexité entre deux puissances successives de 2. Récemment, van der Hoeven a donné un nouvel algorithmme, appelé TFT (de Truncated Fourier Transform en anglais), qui coïncide avec la DFT si $n$ est une puissance de 2 et a une courbe de complexité « plus lisse » que la DFT [Hoe04].

Le produit de $t$ polynômes de degrés arbitraires $d_1, \dots, d_t$ peut s'effectuer [Str83] en $O(M(n)(1+H(d_1, \dots, d_t)))$ opérations, où $H(d_1, \dots, d_t) = -\sum_i \frac{d_i}{n} \cdot \log \frac{d_i}{n}$ est l'entropie du vecteur $(d_1/n, \dots, d_t/n)$. La multiplication des polynômes à plusieurs variables peut se ramener à une variable grâce à une technique connue sous le nom de substitution de Kronecker [Moe76] (voir aussi l'Exercice 5.9).

### Bibliographie

Bera BERNSTEIN, Daniel J. (2001). *Multidigit multiplication for mathematicians*. URL: http://cr.yp.to/papers.html (visible en 2001).

BL04 BÜRGISser, Peter et Martin Lotz (2004). « Lower bounds on the bounded coefficient complexity of bilinear maps ». In: *Journal of the Association for Computing Machinery*, vol. 51, n°3, p. 464–482.

Bri74 BRIGHAM, E. Oran (1974). *The fast Fourier transform*. Prentice-Hall.

BZ07 BODRATO, Marco et Alberto ZANONI (2007). « Integer and polynomial multiplication : towards optimal Toom–Cook matrices ». In: *ISSAC’07 : International Symposium on Symbolic and Algebraic Computation*. Éd. par Dongming WANG. ACM Press, p. 17–24.

CA69 COOK, S. A. et S. O. AANDERAA (1969). « On the minimum computation time of functions ». In: *Transactions of the American Mathematical Society*, vol. 142, p. 291–314.

CH07 CHUNG, Jaewook et M. Anwar HASAN (2007). « Asymmetric squaring formulae ». In: *ARITH’07 : IEEE Symposium on Computer Arithmetic*. IEEE, p. 113–122.

CK91 CANTOR, D. G. et E. KALTOFEN (1991). « On fast multiplication of polynomials over arbitrary algebras ». In: *Acta Informatica*, vol. 28, n°7, p. 693–701.
CKO09 CENK, M., C. K. KOC et F. OZBUDAK (2009). « Polynomial multiplication over finite fields using field extensions and interpolation ». In : *ARITH'09: IEEE Symposium on Computer Arithmetic*. IEEE Computer Society, p. 84–91.

Coo66 COOK, S. (1966). « On the minimum computation time of functions ». Thèse de doctorat. Harvard University.

Coo90 COOLEY, James W. (1990). « How the FFT gained acceptance ». In : *A history of scientific computing*. ACM Press Hist. Ser. P : ACM, p. 133–140.

CT65 COOLEY, James W. et John W. TUKEY (1965). « An algorithm for the machine calculation of complex Fourier series ». In : *Mathematics of Computation*, vol. 19, p. 297–301.

CT93 — (1993). « On the origin and publication of the FFT paper ». In : *Current Contents*, vol. 33, n°51–52, p. 8–9.

De+08 DE, A., P. P. KURUR, C. SAHA et R. SAPTHARISHI (2008). « Fast integer multiplication using modular arithmetic ». In : *STOC'08: ACM Symposium on Theory of Computing*. Victoria, British Columbia, Canada : ACM, p. 499–506.

Dem97 DEMAZURE, Michel (1997). *Cours d'algèbre*. Nouvelle Bibliothèque Mathématique n°1. Cassini, Paris.

DS00 DONGARRA, J. et F. SULLIVAN (2000). « Top ten algorithms ». In : *Computing in Science & Engineering*, vol. 2, n°1, p. 22–23.

DV90 DUHAMEL, P. et M. VETTERLI (1990). « Fast Fourier transforms : a tutorial review and a state of the art ». In : *Signal Processing. An Interdisciplinary Journal*, vol. 19, n°4, p. 259–299.

FJ05 FRIGO, M. et S. G. JOHNSON (2005). « The design and implementation of FFTW3 ». In : *Proceedings of the IEEE*, vol. 93, n°2, p. 216–231.

Für07 FÜRER, Martin (2007). « Faster integer multiplication ». In : *STOC'07: ACM Symposium on Theory of Computing*. New York : ACM, p. 57–66.

Für09 — (2009). « Faster integer multiplication ». In : *SIAM Journal on Computing*, vol. 39, n°3, p. 979–1005.

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

GS66 GENTLEMAN, W. M. et G. SANDE (1966). « Fast Fourier transforms : for fun and profit ». In : *AFIPS'66*. San Francisco, California : ACM, p. 563–578.

HHL14 HARVEY, David, Joris van der HOEVEN et Grégoire LECERF (2014). « Faster polynomial multiplication over finite fields ». http://arxiv.org/abs/1407.3361.

HHL16 — (2016). « Even faster integer multiplication ». In : *Journal of Complexity*, vol. 36, p. 1–30.

HJB85 HEIDEMAN, Michael T., Don H. JOHNSON et C. Sidney BURRUS (1985). « Gauss and the history of the fast Fourier transform ». In : *Archive for History of Exact Sciences*, vol. 34, n°3, p. 265–277.

Hoe04 HOEVEN, Joris van der (2004). « The truncated Fourier transform and applications ». In : *ISSAC'04: International Symposium on Symbolic and Algebraic Computation*. New York : ACM Press, p. 290–296.
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.

Kam85 KAMINSKI, Michael (1985). « A lower bound for polynomial multiplication ». In : *Theoretical Computer Science*, vol. 40, n°2-3, p. 319-322.

KO63 KARATSUBA, A. et Y. OFMAN (1963). « Multiplication of multidigit numbers on automata ». In : *Soviet Physics Doklady*, vol. 7, p. 595-596.

Moe76 MOENCK, Robert T. (1976). « Another polynomial homomorphism ». In : *Acta Informatica*, vol. 6, n°2, p. 153-169.

Mon05 MONTGOMERY, Peter L. (2005). « Five, six, and seven-term Karatsuba-like formulae ». In : *IEEE Transactions on Computers*, vol. 54, n°3, p. 362-369.

Mor73 MORGENSTERN, Jacques (1973). « Note on a lower bound of the linear complexity of the fast Fourier transform ». In : *Journal of the Association for Computing Machinery*, vol. 20, p. 305-306.

Mul00 MULDERS, Thom (2000). « On short multiplications and divisions ». In : *Applicable Algebra in Engineering, Communication and Computing*, vol. 11, n°1, p. 69-88.

Nus80 NUSSBAUMER, Henri J. (1980). « Fast polynomial transform algorithms for digital convolution ». In : *IEEE Transactions on Acoustics, Speech, and Signal Processing*, vol. 28, n°2, p. 205-215.

Nus81 — (1981). *Fast Fourier transform and convolution algorithms*. Vol. 2. Springer Series in Information Sciences. Springer-Verlag.

PFM74 PATERSON, M. S., M. J. FISCHER et A. R. MEYER (1974). « An improved overlap argument for on-line multiplication ». In : *Complexity of computation*. AMS, p. 97-111.

Sch77 SCHÖNHAGE, A. (1976/77). « Schnelle Multiplikation von Polynomen über Körpern der Charakteristik 2 ». In : *Acta Informatica*, vol. 7, n°4, p. 395-398.

SS71 SCHÖNHAGE, A. et V. STRASSEN (1971). « Schnelle Multiplikation großer Zahlen ». In : *Computing*, vol. 7, p. 281-292.

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

Too63 TOOM, A. L. (1963). « The complexity of a scheme of functional elements simulating the multiplication of integers ». In : *Doklady Akademii Nauk SSSR*, vol. 150, p. 496-498.

Van92 VAN LOAN, Charles (1992). *Computational frameworks for the fast Fourier transform*. Vol. 10. Frontiers in Applied Mathematics. SIAM.

WP06 WEIMERSKIRCH, André et Christof PAAR (2006). « Generalizations of the Karatsuba algorithm for efficient implementations ». Cryptology ePrint Archive. URL : https://eprint.iacr.org/2006/224.pdf.
## 3. Calculs rapides sur les séries

### Résumé

La multiplication rapide de polynômes vue au chapitre précédent permet des calculs efficaces comme l'inversion, l'exponentiation ou la prise de logarithme d'une série. Ces opérations sont effectuées grâce à une version formelle de la méthode de Newton. La composition de séries est en général plus coûteuse que le produit, mais peut aussi tirer parti d'une multiplication rapide.

Ce chapitre porte sur le calcul des premiers termes de développements en séries. Pour fixer les notations, N sera utilisé pour représenter le nombre de termes à calculer (en comptant les coefficients nuls), et « série » sera employé pour « série formelle tronquée à l'ordre N ». Ainsi, « calculer une série » voudra toujours dire en calculer les N premiers termes. Comme dans le Chapitre 2, l'expression M(N) dénote une borne supérieure sur le nombre d'opérations arithmétiques nécessaires pour multiplier deux polynômes de degré inférieur à N. La fonction M étant une fonction de multiplication polynomiale, rappelons qu'elle est supposée vérifier l'inégalité $M(N + N') \geq M(N) + M(N')$. L'efficacité des algorithmes sera mesurée par leur complexité arithmétique.

L'opérateur de Newton est très largement utilisé en calcul numérique, pour trouver des solutions approchées d'équations. C'est un processus itératif, qui consiste à chaque étape à remplacer la fonction dont on cherche un zéro par son linéarisé au voisinage de la dernière approximation trouvée. Ainsi, si $\Phi$ est une fonction $C^1$ de $\mathbb{R}$ dans $\mathbb{R}$, pour résoudre $\Phi(y) = 0$ dans $\mathbb{R}$, on choisit $y_0 \in \mathbb{R}$ et on itère

$$y_{k+1} = \mathcal{N}(y_k) = y_k - \frac{\Phi(y_k)}{\Phi'(y_k)}. \qquad (3.1)$$

Cette itération est représentée en Figure 3.1. La solution serait exacte si $\Phi$ était une fonction linéaire; dans le cas général, si $y_0$ est bien choisi (par exemple assez
<figure><img src="image_10.png" /><figcaption>FIGURE 3.1 – Opérateur de Newton de $\mathbb{R}$ dans $\mathbb{R}$.</figcaption></figure>

proche d'une racine simple), l'itération est bien définie ($\Phi'(y_k)$ n'est jamais nul) et la suite $y_k$ converge vers une limite $y_\infty$. En outre, à l'itération $k$ la distance à la limite est de l'ordre du carré de cette distance à l'itération $k-1$. On parle dans ce cas de convergence *quadratique*. Dans ces conditions, pour $k$ suffisamment grand, le nombre de décimales correctes est approximativement *doublé* à chaque itération. De nombreuses généralisations de cette méthode existent, d'abord à des systèmes (la division par l'inverse est alors remplacée par une multiplication à gauche par l'inverse de la matrice jacobienne), puis à des espaces de Banach (où la matrice jacobienne est remplacée par la différentielle de l'opérateur).

Dans le cadre formel, la fonction $\Phi$ est remplacée par une application de l'algèbre des séries formelles (définie en Section 3.1) dans elle-même. Cet algèbre n'est pas un espace normé (donc pas un espace de Banach), mais dans bien des cas, l'itération de Newton converge (pour la métrique de la valuation, voir plus loin), et le nombre de coefficients corrects est doublé à chaque itération, ce qui mène à une bonne complexité. La composition des séries reste une opération plus coûteuse, à laquelle est consacrée la fin de ce chapitre.

## 3.1 Séries formelles

Dans cette section et la suivante les séries ont leurs coefficients dans un anneau qui n'est pas nécessairement commutatif, sauf lorsque c'est explicitement indiqué. À partir de la Section 3.3, l'utilisation de la formule de Taylor contraint à se restreindre au cas commutatif.

### Définition

Si $\mathbb{A}$ est un anneau (non nécessairement commutatif), on note $\mathbb{A}[[X]]$ l'ensemble des séries formelles sur $\mathbb{A}$. Ses éléments sont des suites $(f_i)_{i \in \mathbb{N}}$ de $\mathbb{A}$, notées

$$F(X) = \sum_{i \ge 0} f_i X^i.$$
Le coefficient $f_i$ est appelé le *i*-ième coefficient de $F(X)$, le coefficient $f_0$ est appelé terme constant de $F(X)$, et parfois noté $F(0)$. L'indice du premier coefficient non nul est appelé la *valuation* de $F$ et noté val $F$. Par convention, val $0 = \infty$.

Les opérations de $\mathbb{A}[[X]]$ sont l'addition des suites et une multiplication (appelée parfois produit de Cauchy) qui généralise la multiplication des polynômes :

$$ \sum_{i \ge 0} f_i X^i \times \sum_{i \ge 0} g_i X^i = \sum_{i \ge 0} h_i X^i, \quad \text{avec} \quad h_i = \sum_{j+k=i} f_j g_k, $$

la dernière somme étant finie.

**Exemple 3.1** La série formelle

$$ 1 = 1 \cdot X^0 + 0 \cdot X + 0 \cdot X^2 + \dots, $$

où 1 est l'unité de $\mathbb{A}$, est élément neutre pour la multiplication de $\mathbb{A}[[X]]$. La formule donnant les coefficients du produit se réduit alors à un terme.

**Exemple 3.2** Si $F = 1+X$ et $G = 1-X+X^2-X^3+\dots$, alors le produit vaut $H = FG = 1$ : le coefficient de $X^n$ dans $H$ vaut $1-1=0$ pour $n>0$ et 1 sinon. La série $G$ est donc l'inverse de $F$ pour la multiplication, que l'on peut noter $(1+X)^{-1}$.

### Métrique

Lorsque les coefficients sont des nombres complexes, la question de la convergence des séries entières représentées par ces séries formelles ne se posera pas pour les algorithmes considérés dans ce chapitre. En revanche, quel que soit l'anneau de coefficients $\mathbb{A}$, il est possible de définir une distance entre deux séries F et G par $d(F,G) = 2^{-\text{val}(F-G)}$.

**Exercice 3.1** Vérifier que la fonction $d$ définie ci-dessus est bien une distance. ■

Muni de cette distance, l'ensemble des séries formelles forme un espace métrique complet (les suites de Cauchy convergent). En effet, dire qu'une suite $(S_n)$ de séries est de Cauchy signifie qu'étant donné $k \in \mathbb{N}$, il existe un entier $K$ tel que pour tous $m, n \ge K$, on ait $d(S_m, S_n) < 2^{-k}$, autrement dit après l'indice $K$, les $k$ premiers termes des $S_n$ sont fixés, ce sont ceux de la série limite.

### Composition

Si F et G sont deux séries formelles, avec terme constant $G(0) = 0$, on définit la composition comme

$$ F(G(X)) = f_0 + f_1 G(X) + f_2 G(X)^2 + \dots $$

Les points de suspension signifient que l'on considère la limite des sommes $H_n$ obtenues en arrêtant la somme au terme $f_n G(X)^n$. Cette limite existe puisque pour $m \ge n$, la distance obéit à $d(H_m, H_n) \le 2^{-n}$, ce qui fait de $(H_n)_n$ une suite de Cauchy.
Inverse

Si F est de la forme $F = 1 + XG$ avec $G \in A[[X]]$, alors F est inversible et son inverse est donné par

$$1 - XG + X^2G^2 - \dots,$$

composition de $(1 + X)^{-1}$ par GX. La preuve découle de cette composition :

$$(1 + XG)(1 - XG + \dots) = H(XG)$$

où $H = (1 + X)(1 + X)^{-1} = 1$.

Si $a = F(0)$ est inversible, la série F se récrit $F = a(1 + XG)$ avec $G = a^{-1}(F - a)/X$, donc F est inversible, d'inverse

$$(1 + XG)^{-1}a^{-1} = a^{-1} - Ga^{-1}X + G^2a^{-1}X^2 + \dots$$

Ces ingrédients munissent l'ensemble des séries formelles d'une structure d'anneau.

**Proposition 3.1** L'ensemble $A[[X]]$ des séries formelles à coefficients dans $A$ est un anneau, qui est commutatif si $A$ l'est. Ses éléments inversibles sont les séries de terme constant inversible.

*Démonstration.* L'addition est commutative comme celle de $A$. L'associativité et la distributivité du produit sont obtenues comme pour les polynômes. L'unité pour le produit est la série 1. La première partie de la proposition est donc prouvée.

Si $F = \sum f_i X^i \in A[[X]]$ est inversible, et $G = \sum g_i X^i$ est son inverse, alors l'extraction du coefficient de $X^0$ dans l'identité $FG = 1$ donne $f_0 g_0 = 1$ ce qui prouve qu'une série inversible a un terme constant inversible. La réciproque a été présentée ci-dessus. ■

Séries à plusieurs variables

Comme $A[[X]]$ est un anneau, il peut être utilisé comme anneau de base pour définir des séries formelles en une autre variable Y, ce qui définit de la même manière l'anneau des séries formelles à deux variables $A[[X]][[Y]]$, noté aussi $A[[X, Y]]$.

Dérivation

La dérivée d'une série est définie formellement coefficient par coefficient via l'identité

$$\left( \sum_{i \ge 0} f_i X^i \right)' = \sum_{i \ge 0} (i+1) f_{i+1} X^i.$$

Les relations habituelles $(F + G)' = F' + G'$ et $(FG)' = F'G + FG'$ sont prouvées comme pour les polynômes. Dans le cas de séries en plusieurs variables, on utilise la notation des dérivées partielles $\partial F / \partial X$ pour désigner la dérivée par rapport à la variable X.

Troncatures

Algorithmiquement, on ne manipule pas de série à précision « infinie », mais seulement des troncatures, c'est-à-dire un certain nombre des premiers termes. Les séries tronquées deviennent alors de simples polynômes, pour lesquels on dispose en particulier d'algorithmes rapides de multiplication. Étant donnés les N premiers
termes d'une série F, l'objectif de ce chapitre est de donner des algorithmes permettant
de calculer efficacement les N premiers termes d'autres séries définies à partir de F.

Pour étendre la notation utilisée avec les polynômes, étant donnée la série

$$
S = \sum_{i \ge 0} a_i X^i,
$$

on notera S mod XN le polynôme

$$
S \bmod X^N := \sum_{0 \le i < N} a_i X^i.
$$

On notera $S = f + O(X^N)$ si les séries ou polynômes $S$ et $f$ coïncident jusqu'au terme de degré $N-1$, et même plus généralement $S = O(T^k)$ si $S = O(X^{k \text{val}(T)})$.

Formule de Taylor

Dans le cas où l'anneau des coefficients $\mathbb{A}$ est commutatif, le lien entre la composition et la dérivation est donné par la formule de Taylor. Si $F, G, H$ sont trois séries formelles, avec $G(0) = H(0) = 0$, et les entiers $2, 3, \dots, k-1$ sont inversibles dans $\mathbb{A}$, alors

$$
F(G + H) = F(G) + F'(G)H + F''(G)\frac{H^2}{2!} + \dots + O(H^k).
$$

La formule est classique pour les polynômes, et les coefficients de X^N dans les deux
membres de la formule sont les mêmes que ceux de l'identité entre polynômes obtenue
en considérant les troncatures de F, G et H modulo X^{N+1}, ce qui permet de conclure.

Pour cette identité, la commutativité de $\mathbb{A}$ est cruciale comme le montre l'exemple
de $F = X^2$ : on a alors $(G+H)^2 = G^2 + GH + HG + H^2$ et la partie linéaire $GH + HG$ est
en général différente de $2GH = F'(G)H$.

Intégration

Il s'agit de l'opération inverse de la dérivation. On suppose que les entiers sont
inversibles dans $\mathbb{A}$ ($\mathbb{A}$ est une $\mathbb{Q}$-algèbre) et on définit alors

$$
\int \sum_{i \ge 0} f_i X^i = \sum_{i \ge 0} f_i \frac{X^{i+1}}{i+1}.
$$

Exemple 3.3 La série log(1 + X) vaut

$$
\log(1 + X) = \int (1 + X)^{-1} = X - \frac{1}{2}X^2 + \frac{1}{3}X^3 + \dots
$$

3.2 La méthode de Newton pour le calcul d'inverses

Partant de la fonction $\Phi(y) = 1/y - f$, dont l'unique racine est l'inverse $1/f$, l'itération (3.1) devient

$$
\mathcal{N}(g_k) = g_k + g_k^2 (1/g_k - f) = g_k + g_k (1 - g_k f).
$$

Cette itération permet le calcul d'inverses dans plusieurs contextes.
**Entrée** Un entier $N > 0$, F mod $X^N$ une série tronquée.

**Sortie** $F^{-1}$ mod $X^N$.

Si $N = 1$, alors renvoyer $f_0^{-1}$, où $f_0 = F(0)$.

Sinon :

1. Calculer récursivement l'inverse $G$ de $F$ mod $X^{\lceil N/2 \rceil}$.

2. Renvoyer $G + (1 - GF)G$ mod $X^N$.

**Algorithme 3.2** – Inverse de série par itération de Newton.

## Convergence quadratique pour l'inverse d'une série formelle

**Lemme 3.2** Soient $\mathbb{A}$ un anneau non nécessairement commutatif, $F \in \mathbb{A}[[X]]$ une série formelle de terme constant inversible et $G$ une série telle que $G - F^{-1} = O(X^n)$ ($n \ge 1$). Alors la série

$$ \mathcal{N}(G) = G + (1 - GF)G \quad (3.2) $$

vérifie $\mathcal{N}(G) - F^{-1} = O(X^{2n})$.

*Démonstration.* Par hypothèse, on peut définir $H \in \mathbb{A}[[X]]$ par $1 - GF = X^nH$. Il suffit alors de récrire $F = G^{-1}(1 - X^nH)$ et d'inverser :

$$ F^{-1} = (1 + X^n H + O(X^{2n}))G = G + X^n HG + O(X^{2n})G = \mathcal{N}(G) + O(X^{2n}). \quad \blacksquare $$

## Algorithme

**Lemme 3.3** L'Algorithme 3.2 d'inversion est correct.

*Démonstration.* La preuve est une récurrence sur les entiers. Pour $N = 1$ la propriété est claire. Si la propriété est vraie pour $1 \le N < k$, alors elle l'est pour $1 \le \lceil N/2 \rceil < k$ donc le résultat de l'étape (1) vaut $F \mod X^{\lceil N/2 \rceil}$. Alors, d'après le Lemme 3.2, $G + (1 - GF)G = F^{-1} \mod X^{2\lceil N/2 \rceil}$ et la troncature modulo $X^N$ est donc bien correcte. ■

## Complexité

Chaque itération coûte essentiellement deux multiplications, et l'ensemble est donc effectué en $O(M(N))$ opérations d'après le lemme « diviser pour régner ». Une estimation plus précise est donnée par la proposition suivante.

**Proposition 3.4** Soient $\mathbb{A}$ un anneau non nécessairement commutatif, F dans $\mathbb{A}[[X]]$ avec $F(0)$ inversible et $N \ge 1$. On peut calculer l'inverse de F modulo $X^N$ en $3M(N) + O(N)$ opérations dans $\mathbb{A}$ lorsque $N$ est une puissance de 2.

*Démonstration.* L'étape (2) de l'algorithme demande deux multiplications : comme $G$ est de degré inférieur à $\lceil N/2 \rceil$, celle de $G$ par $F$ coûte $2M(\lceil N/2 \rceil)$ et produit un résultat de la forme $X^{\lceil N/2 \rceil}R_N$ avec $R_N$ de degré inférieur à $\lceil N/2 \rceil$. Sa multiplication par $G$
coûte alors $M(\lceil N/2 \rceil)$ opérations. S'ajoutent à ces multiplications $\lambda \cdot N$ opérations supplémentaires (additions, ...), $\lambda$ étant une constante que l'on ne va pas déterminer. Notons $C(N)$ le coût du calcul modulo $N$; on a alors

$$C(N) \leq C(\lceil N/2 \rceil) + 3M(\lceil N/2 \rceil) + \lambda N.$$

L'application du lemme « diviser pour régner » donne alors

$$C(N) \leq 3M(N) + 2\lambda N$$

pour $N$ puissance de 2 en utilisant l'inégalité $M(d) \leq \frac{1}{2}M(2d)$.

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

Comme d'habitude, lorsque $N$ n'est pas une puissance de 2, la complexité est bornée par celle de la puissance de 2 immédiatement supérieure, ce qui rajoute au plus un facteur 2.

### Optimisations

Le coût de la Proposition 3.4 peut être encore amélioré. Le produit GF est de la forme $1 + X^{\lceil N/2 \rceil}R_N$. Il s'agit donc de calculer le produit d'un polynôme de degré $N/2$ par un polynôme de degré $N$, alors que l'on connaît à l'avance les $N/2$ premiers coefficients du produit. Dans cette situation, il est possible d'abaisser le coût du produit de $M(N)$ à $M(N/2)$, en utilisant un produit médian et des techniques non triviales de « transposition » de code abordées au Chapitre 12.

### Division de séries

Le quotient $H/F$ où $H$ et $F$ sont des séries et $F$ est inversible, peut être obtenu comme le produit $H \times F^{-1}$. Ce calcul peut être amélioré d'un facteur constant. Il suffit pour cela de remplacer la dernière itération $G + (1 - GF)G$ par $Y + (H - YF)G$ où $Y := HG \mod X^{\lceil N/2 \rceil}$. Cette astuce remplace un produit de la forme $(N/2, N)$ (pour GF), un produit de la forme $(N/2, N/2)$ (pour $(1 - GF)G$ et un produit final $(N, N)$ par deux produits $(N/2, N/2)$ (pour $Y$ et $(H - YF)G$) et un produit $(N/2, N)$ pour calculer $YF$. L'économie est donc la différence $M(N) - M(N/2)$.

### Application à la division euclidienne

L'Algorithmme 3.2 est l'élément clé de la division euclidienne rapide, traitée au Chapitre 4.

### Application au logarithme

Si $F \in \mathbb{A}[[X]]$ est telle que $F(0) = 0$, on définit la série $\log(1 + F)$ par composition avec $\log(1 + X)$. Pour calculer cette série lorsque $\mathbb{A}$ est commutatif, il suffit d'utiliser l'identité

$$\log(1 + F) = \int \frac{F'}{1 + F}.$$

**Exercice 3.2** Prouver cette identité à partir des définitions de $\log$, $\int$ et de la dérivation données plus haut. ■
Le calcul demande une division de séries, une dérivation et une intégration, mais ces deux dernières opérations sont de complexité linéaire. Le bilan est donc une complexité en O(M(N)).

Inverse de matrices

L'anneau A n'étant pas supposé commutatif, il est possible de traiter directement des matrices de séries comme des séries de matrices : le Lemme 3.2 et l'Algorithmme 3.2 s'appliquent. La complexité de cet algorithme s'exprime à l'aide de la fonction MM(k, N) complexité du produit de matrices k×k de polynômes de degré au plus N. Le Corollaire 5.11 du Chapitre 5 montre que MM(k, N) = O(k<sup>θ</sup>N + k<sup>2</sup>M(N)).

**Proposition 3.5** Soient A un anneau non nécessairement commutatif, A une matrice k × k à coefficients dans A[[X]], telle que A(0) soit inversible. La matrice A<sup>-1</sup> mod X<sup>N</sup> peut être calculée en O(MM(k, N)) opérations dans A.

Démonstration. Comme pour l'inverse dans l'Algorithmme 3.2, le coût est dominé à
une constante près par celui du calcul de G + (I - GF)G mod XN : la complexité vérifie
l'inégalité

$$
C(N) \leq C(\lceil N/2 \rceil) + 2MM(k, N)
$$

et le résultat se déduit de l'hypothèse MM(k, N) ≤ 1/2 MM(k, 2N).

3.3 Itération de Newton formelle et applications

Dans toute la suite de ce chapitre, l'anneau A est supposé commutatif.

Un résultat général

**Théorème 3.6** — Itération de Newton sur les séries. Soit $\Phi \in A[[X, Y]]$ une série à deux variables en $X$ et $Y$, telle que $\Phi(0, 0) = 0$ et $\frac{\partial\Phi}{\partial Y}(0, 0)$ est inversible dans $A$. Il existe alors une unique série $S \in A[[X]]$, telle que $\Phi(X, S) = 0$ et $S(0) = 0$. Si $F$ est une série telle que $S - F = O(X^n) (n \ge 1)$, alors

$$
\mathcal{N}(F) = F - \frac{\Phi(X, F)}{\frac{\partial \Phi}{\partial Y}(X, F)}
$$

vérifie $S - \mathcal{N}(F) = O(X^{2n})$.

Ce résultat est un résultat local en (0,0). Il exprime que si le point (0,0) est solution,
il s'étend en une courbe solution. Par translation, d'autres situations où la solution
n'est pas en (0,0) s'y ramènent. La première partie est une version du théorème des
fonctions implicites pour les séries formelles.

Démonstration. D'abord il faut observer que l'itération est bien définie : la composition de Φ et de dΦ/dY avec F sont possibles parce que F(0) = S(0) = 0. Pour la même raison,
**Entrée** Un entier $N > 0$, $\Phi \bmod (X^N, Y^N)$ une série tronquée avec $\Phi(0, 0) = 0$ et $(\partial\Phi/\partial Y)(0, 0)$ inversible.

**Sortie** S mod $X^N$, T mod $X^{\lceil N/2 \rceil}$, telles que
$S(0) = 0$, $\Phi(X, S) = 0$, et $T = (\partial\Phi/\partial Y(X, S))^{-1}$.

1. Si $N = 1$, alors renvoyer 0, $((\partial\Phi/\partial Y)(0, 0))^{-1}$.

2. Sinon :

a. Calculer récursivement $F := S \bmod X^{\lceil N/2 \rceil}$ et $G := T \bmod X^{\lceil N/2 \rceil}$.

b. Calculer $G := G + (1 - G \frac{\partial \Phi}{\partial Y}(X, F))G \bmod X^{\lceil N/2 \rceil}$.

c. Calculer $F := F - G\Phi(X, F) \bmod X^N$.

d. Renvoyer $F, G$.

**Algorithme 3.3 – Résolution de $\Phi(X, Y) = 0$ par itération de Newton.**

$\frac{\partial \Phi}{\partial Y}(X, F)$ est inversible puisque son terme constant $\frac{\partial \Phi}{\partial Y}(0, 0)$ est inversible. De ce fait, on déduit aussi $\text{val}(\mathcal{N}(F) - F) = \text{val}(\Phi(X, F))$.

Ensuite, la preuve se déroule en deux temps : l'itération est d'abord utilisée pour montrer l'existence d'une série solution $S$, puis on montre l'unicité en même temps que la vitesse de convergence vers $S$.

La valuation de $\Phi(X, F)$ est doublée par l'itération : la formule de Taylor donne

$$ \Phi(X, \mathcal{N}(F)) = \Phi(X, F) + \frac{\partial \Phi}{\partial Y}(X, F)(\mathcal{N}(F) - F) + O((\mathcal{N}(F) - F)^2) = O(\Phi(X, F)^2). \quad (3.3) $$

La suite définie par $S_0 = 0$ et $S_{k+1} = \mathcal{N}(S)$ vérifie donc

$$ \text{val}(S_{k+2} - S_{k+1}) = \text{val}(\Phi(X, S_{k+1})) \geq 2 \text{val}(S_{k+1} - S_k). $$

Cette suite est donc de Cauchy et sa limite $S$ existe et vérifie $\Phi(X, S) = \lim(\Phi(X, S_k)) = 0$. La vitesse de convergence vers $S$ vient encore de la formule de Taylor. En effet, l'égalité

$$ \Phi(X, F) = \Phi(X, S) + \frac{\partial \Phi}{\partial Y}(X, F)(S - F) + O((S - F)^2) \quad (3.4) $$

donne $S - F = O(\Phi(X, F))$. On en déduit $S - \mathcal{N}(F) = O(\Phi(X, \mathcal{N}(F))) = O(X^{2n})$ au vu des hypothèses et de (3.3). Cette égalité donne aussi l'unicité de la solution : si $F$ est une autre solution avec $F(0) = 0$, elle doit vérifier $S - F = O((S - F)^2)$, ce qui entraîne $\text{val}(S - F) = \infty$ et donc $F = S$. ■

## Algorithme

**Lemme 3.7** L'Algorithme 3.3 de résolution est correct.

*Démonstration.* Pour $N = 1$, les valeurs renvoyées sont les termes constants de $S$ et $T$. Si la propriété est vraie pour $1 \le k < N$, alors l'étape (1.a) renvoie $T \bmod X^{\lceil N/2 \rceil}$ et
S mod X[N/2]. D'après le Lemme 3.2, l'étape (1.b) calcule T mod X[N/2]. L'étape (1.c) calcule alors S mod XN : les formules (3.3) et (3.4) montrent que seuls les [N/2] premiers coefficients de ∂Φ/∂Y(X, S) sont utiles pour doubler la précision. ■

Applications

Un exemple d'application est fourni par l'inverse de la section précédente avec Φ(X, Y) = f − 1/(f (0) − 1 + Y) (lorsque A est commutatif). La complexité de l'Algo- rithme 3.3 dépend en particulier de celle du calcul de Φ et de ∂Φ/∂Y . Les applications les plus directes de la méthode de Newton sont résumées dans les théorèmes suiv- vants.

Théorème 3.8 Dans le cas où la série Φ du Théorème 3.6 est un polynôme en Y de degré O(1), les N premiers termes de sa série solution S, telle que S(0) = 0, sont obtenus en O(M(N)) opérations dans A.

Démonstration. L'argument est le même que pour la complexité de l'inverse. Pour Φ polynôme de degré O(1), les évaluations de Φ et de Φ' coûtent O(M(N)), le reste de l'argument est inchangé. ■

Dans le cas où Φ est un polynôme à la fois en X et en Y, un algorithme de com- plexité O(N) (donc linéaire au lieu de quasi-linéaire, et ce, quelle que soit la multipli- cation utilisée) sera présenté au Chapitre 14.

Théorème 3.9 Soit f une série de A[[X]]. On peut calculer en O(M(N)) opérations dans A les N premiers termes de :

1. f − $^{1}$, lorsque f (0) est inversible;

2. log f et f α lorsque f (0) = 1 et α ∈ K ;

3. exp(f ) lorsque f (0) = 0.

Pour les points (2) et (3), nous supposons également que 2, 3, ..., N − 1 sont des éléments inversibles de A.

Démonstration. La preuve pour l'inverse et le logarithme provient des Sections 3.2 et 3.2. L'exponentielle est traitée dans la section qui vient, et la puissance se déduit des précédentes par

$f^{\alpha} = \exp(\alpha \log f).$

Une application remarquable de cette identité évidente apparaît dans la Section 3.3 sur l'inverse compositionnel. ■

Exponentielle

Si F ∈ A[[X]] est telle que F(0) = 0, on définit la série exp(F) par

$$
\exp(F) = 1 + F + \frac{F^2}{2!} + \frac{F^3}{3!} + \dots
$$

lorsque 2, 3, ... sont inversibles dans A. Pour calculer cette série, l'idée est d'appliquer une méthode de Newton avec Φ(Y) = F − log Y, Y(0) = 1.
**Entrée** Un polynôme $P \in \mathbb{K}[X]$ de degré $d$, un entier $k$.

**Sortie** Les sommes de Newton $p_0, \dots, p_k$ de P.

1. Former les polynômes $N = T^{d-1}P'(1/T)$ et $D = T^d P(1/T)$.

2. Calculer $N/D \bmod T^{k+1} = p_0 + p_1 T + \cdots + p_k T^k$.

3. Renvoyer les coefficients.

**Algorithme 3.4** – Calcul des sommes de Newton à partir des coefficients.

**Entrée** Les sommes de Newton $p_0, \dots, p_d$ du polynôme unitaire $P \in \mathbb{K}[X]$.

**Sortie** Les coefficients de P.

1. Former la série $S_1 := p_1 + p_2 T + \cdots + p_d T^{d-1}$.

2. Calculer $S_2 := \exp(-\int S_1) \bmod T^{d+1}$.

3. Renvoyer $X^{p_0} S_2(1/X)$.

**Algorithme 3.5** – Calcul des coefficients à partir des sommes de Newton.

**Exercice 3.3** Montrer que les conditions du Théorème 3.6 sont vérifiées. ■

Il s'ensuit une convergence quadratique vers $\exp(F)$ pour l'itération

$$E_{k+1} = \mathcal{N}(E_k) = E_k + E_k(F - \log E_k).$$

Chaque itération demande donc le calcul d'une multiplication et d'un logarithme, d'où la complexité en $O(M(N))$. Une autre itération légèrement plus efficace est présentée en Exercice 3.8.

**Exercice 3.4** Donner une itération de Newton qui calcule directement la racine carrée, sans passer par l'exponentielle et le logarithme. Plus difficile, estimer le gain par rapport à l'algorithme général de calcul de puissance. ■

**Exercice 3.5** Donner des bornes sur les constantes dans les $O(\cdot)$ pour le logarithme, l'exponentielle et la puissance. ■

## Sommes de Newton

Une conséquence intéressante de l'efficacité du calcul de l'exponentielle est l'utilisation efficace des sommes de Newton comme structure de données.

Si $P \in \mathbb{K}[X]$ est un polynôme de degré $d$, avec $\mathbb{K}$ un corps, alors $P$ a $d$ racines $\alpha_1, \dots, \alpha_d$ dans la clôture algébrique $\bar{\mathbb{K}}$ de $\mathbb{K}$. Rappelons qu'un polynôme est dit unitaire lorsque son coefficient de plus haut degré est 1. Les sommes de Newton de P sont les sommes $p_i = \alpha_i^i + \cdots + \alpha_d^i$. Ce sont des éléments de $\mathbb{K}$, puisqu'ils sont fonctions symétriques des racines de P. (Tout polynôme symétrique en les racines d'une équation algébrique s'exprime comme un polynôme en les coefficients de l'équation.) Leur utilisation efficace repose sur la proposition suivante.
**Proposition 3.10** Les sommes de Newton $p_1, \dots, p_d$ d'un polynôme $P \in \mathbb{K}[X]$ de degré $d$ peuvent être calculées par l'Algorithme 3.4 en $O(M(d))$ opérations dans $\mathbb{K}$. Lorsque la caractéristique de $\mathbb{K}$ est 0 ou strictement supérieure à $d$, la conversion inverse est également possible en $O(M(d))$ opérations par l'Algorithme 3.5.

*Démonstration.* Si $P = c \prod_{\alpha} (X - \alpha)$, la décomposition en éléments simples

$$S = \frac{XP'}{P} = \sum_{P(\alpha)=0} \frac{X}{X-\alpha} = \sum_{\substack{P(\alpha)=0 \\ i \ge 0}} \alpha^i X^{-i},$$

où $\alpha$ décrit les zéros de P comptés avec multiplicité, montre que les coefficients du développement en puissances de $X^{-1}$ de $XP'/P$ sont les $p_i$. Pour se ramener à des calculs de séries formelles, la première étape introduit une nouvelle variable par $X = 1/T$. Le calcul de la dérivée à une complexité linéaire en $d$ ; les conversions en puissances de T ne demandent pas d'opérations arithmétiques. Les $k = d$ premiers termes de la division sont obtenus $O(M(d))$ opérations, comme montré en Section 3.2.

L'opération inverse, à savoir le calcul de P à partir des $d$ premiers termes de la série $S = XP'/P$ est effectué en $O(M(d))$ opérations par la formule $P = \exp \int S/X$. À nouveau, il est plus simple de travailler en T pour rester dans le domaine des séries formelles. Le signe '-' dans le calcul de $S_2$ provient du signe de l'intégration des puissances négatives de X. La contrainte sur la caractéristique permet de définir l'exponentielle tronquée par son développement en série. ■

## Application à la somme et au produit composés

On suppose dans cette section que $\mathbb{K}$ est un corps de caractéristique nulle ou strictement supérieure à $d$.

**Théorème 3.11** Soient P et Q deux polynômes unitaires de $\mathbb{K}[X]$, de degrés $d_p$ et $d_q$ et soit $d = d_p d_q$. Alors, les polynômes

$$P \oplus Q = \prod_{P(\alpha)=0, Q(\beta)=0} (X - (\alpha + \beta)), \quad P \otimes Q = \prod_{P(\alpha)=0, Q(\beta)=0} (X - \alpha\beta)$$

peuvent être calculés en $O(M(d))$ opérations dans $\mathbb{K}$ par les Algorithmes 3.6 et 3.7.

Comme ces polynômes sont de degré $d$, il s'agit là encore d'une complexité quasi-optimale. Dans cet énoncé, les polynômes sont définis à l'aide des racines de P et Q dans une clôture algébrique, mais leurs coefficients sont dans $\mathbb{K}$. Il est également possible de les définir sur un anneau quelconque comme des *résultants* en deux variables. Les résultats et leur calcul sont présentés au Chapitre 6. En général, avec deux variables, il n'existe pas encore d'algorithmes quasi-optimal pour le résultant et le résultat ci-dessus est l'un des rares cas particuliers qui se traite efficacement.

Les polynômes constituent une structure de données bien commode pour coder des nombres algébriques et calculer avec : les polynômes $P \oplus Q$ et $P \otimes Q$ fournissent l'addition et la multiplication avec cette structure de données.
**Entrée** Des polynômes unitaires P et Q de $K[X]$, de degrés $d_p$ et $d_q$.

**Sortie** Leur produit composé $P \otimes Q$ de degré $d := d_p d_q$.

1. Calculer les sommes de Newton $(p_0, \dots, p_{d^2})$ et $(q_0, \dots, q_{d^2})$ de P et Q.

2. Renvoyer le polynôme dont les sommes de Newton sont $(p_0 q_0, \dots, p_{d^2} q_{d^2})$.

### Algorithme 3.6 – Produit composé.

**Entrée** les polynômes unitaires P et Q de $K[X]$, de degrés $d_p$ et $d_q$.

**Sortie** leur somme composée $P \oplus Q$ de degré $d := d_p d_q$.

1. Calculer les sommes de Newton $(p_0, \dots, p_{d^2})$ et $(q_0, \dots, q_{d^2})$ de P et Q.

2. Former les polynômes

$$S_p := p_0 + p_1 \frac{T}{1!} + \dots + p_{d^2} \frac{T^{d^2}}{(d^2)!}, \quad S_q := q_0 + q_1 \frac{T}{1!} + \dots + q_{d^2} \frac{T^{d^2}}{(d^2)!}$$

3. Calculer $S := S_p \times S_q \mod T^{d^2+1} = s_0 + s_1 T + s_2 T^2 + \dots$

4. Renvoyer le polynôme dont les sommes de Newton sont $(s_0, 1! s_1, 2! s_2, \dots, (d^2)! s_{d^2})$.

### Algorithme 3.7 – Somme composée.

*Démonstration.* Le produit de sommes de Newton est la somme de Newton du produit :

$$\sum_{\alpha} \alpha^i \sum_{\beta} \beta^i = \sum_{\alpha, \beta} (\alpha \beta)^i.$$

Il suffit donc de multiplier terme à terme les séries génératrices des sommes de Newton de P et Q pour obtenir la série génératrice des sommes de Newton de $P \otimes Q$. Si l'on a calculé $d$ termes de ces séries, le résultant, dont le degré est $d$, est alors reconstruit par une exponentielle en $O(M(d))$ opérations.

Diviser le $i$-ième coefficient de $XP'/P$ par $i!$, pour tout $i \ge 0$, produit la série

$$\sum_{\substack{P(\alpha)=0 \\ i \ge 0}} \frac{\alpha^i X^{-i}}{i!} = \sum_{P(\alpha)=0} \exp(\alpha/X).$$

En effectuant la même opération sur Q et en multipliant les séries, on obtient donc

$$\sum_{P(\alpha)=0} \exp(\alpha/X) \sum_{Q(\beta)=0} \exp(\beta/X) = \sum_{\substack{P(\alpha)=0 \\ Q(\beta)=0}} \exp((\alpha + \beta)/X).$$
Il suffit alors de remultiplier le i-ième coefficient par i! pour obtenir la série génératrice des sommes de Newton du polynôme P ⊗ Q, qui est reconstruit par une exponentielle. L'ensemble des opérations a une complexité bornée par O(M(d)). ■

## Systèmes

Le Théorème 3.6 s'étend à des systèmes d'équations.

**Théorème 3.12** Soient $\Phi = (\Phi_1, \dots, \Phi_k)$ des séries en $k+1$ indéterminées X et Y = (Y_1, \dots, Y_k), telles que $\Phi(0) = 0$ et la matrice jacobienne $(\frac{\partial \Phi}{\partial Y})$ est inversible dans $\mathbb{A}$ en 0. Alors, le système

$$ \Phi_1(X, Y_1, \dots, Y_k) = \dots = \Phi_k(X, Y_1, \dots, Y_k) = 0 $$

admet une solution $S = (S_1, \dots, S_k) \in \mathbb{A}[[X]]^k$ telle que $S(0) = 0$. Si $F = (F_1, \dots, F_k)$ est tel que $S - F = O(X^n)$ ($n \ge 1$), alors

$$ \mathcal{N}(F) = \begin{pmatrix} F_1 \\ \vdots \\ F_k \end{pmatrix} - \left( \frac{\partial \Phi}{\partial Y}(X, F_1, \dots, F_k) \right)^{-1} \begin{pmatrix} \Phi_1(X, F_1, \dots, F_k) \\ \vdots \\ \Phi_k(X, F_1, \dots, F_k) \end{pmatrix} $$

vérifie $S - \mathcal{N}(F) = O(X^{2n})$.

La preuve est la même que celle du Théorème 3.6, la formule de Taylor pour les fonctions de plusieurs variables s'exprimant à l'aide de la matrice jacobienne.

**Lemme 3.13** L'Algorithmme 3.3 fonctionne aussi, sans changement, pour un système $\Phi$ satisfaisant aux hypothèses du Théorème 3.12.

*Démonstration.* La preuve est la même, en observant que l'étape (1.b) est écrite de telle sorte que le Lemme 3.2 s'applique : cette étape calcule l'inverse de la matrice jacobienne à la précision requise et la dernière étape effectue un produit matrice-vecteur pour en déduire la correction de la nouvelle itération. ■

Le corollaire suivant joue un rôle important dans l'algorithm de résolution géométrique du Chapitre 28.

**Corollaire 3.14 — Systèmes polynomiaux.** Sous les mêmes hypothèses, si $\Phi_i \in \mathbb{A}[[X]]$ $[Y_1, \dots, Y_k]$ pour $1 \le i \le k$, alors les solutions séries à précision N peuvent être calculées en $O(MM(k, N))$ opérations dans $\mathbb{A}$.

*Démonstration.* Dans ce calcul, le nombre d'opérations nécessaires à l'évaluation des polynômes ou de leurs dérivées en une série à précision N est borné par une constante fois MM(N), constante qui ne dépend que des polynômes $\Phi_i$. Ensuite, l'étape (b.) demande des produits de matrices $k \times k$ de polynômes de degré $\lfloor N/2 \rfloor$, et l'étape (c.) demande un produit matrice-vecteur en précision N. La complexité de l'application
récursive est donc dominée par celle de l'étape (2), le résultat est donc le même que pour l'inverse de matrices de séries.

Inverse compositionnel

Étant donnée une série F avec $F(0) = 0$ et $F'(0)$ inversible, il s'agit ici de calculer une série $F^{(-1)}$ telle que $F(F^{(-1)}(X)) = F^{(-1)}(F(X)) = X$.

Les N premiers termes

Les conditions d'application de la méthode de Newton (Théorème 3.6) sont vérifiées pour la fonction

$$ \Phi(X, Y) = F(Y) - X. $$

L'opérateur de Newton correspondant est donc

$$ \mathcal{N}(G) = G - \frac{F(G) - X}{F'(G)}. $$

La dérivation de $F(G) = X$ donne $F'(G)G' = 1$ et cet opérateur se simplifie en

$$ \mathcal{N}(G) = G - G' \times (F(G) - X). $$

Le coût est alors dominé soit par celui du produit, soit plus généralement par celui de la composition par $F$ (voir Section 3.4).

**Exercice 3.6** La k-ième racine positive de l'équation $x = \tan x$ admet un développement asymptotique de la forme

$$ r_k = (2k + 1)\frac{\pi}{2} + \frac{a_1}{k} + \frac{a_2}{k^2} + \dots $$

Pour tout $i$, le coefficient $a_i$ s'écrit $a_i = f_i(1/\pi)$, où $f_i$ est un polynôme impair de degré $2i + 1$. Borner le nombre d'opérations dans $\mathbb{Q}$ nécessaires pour calculer $f_1, \dots, f_N$, pour $N$ grand.

Exemple

La série génératrice exponentielle $T \in \mathbb{Q}[[X]]$ des arbres généraux étiquetés (appelés aussi arbres de Cayley) vérifie l'équation

$$ T = X \exp(T), $$

c'est-à-dire que le coefficient de $X^n$ dans la série $T$ est le nombre d'arbres de taille $n$ divisé par $n!$. Ces coefficients sont obtenus en $O(M(N))$ opérations dans $\mathbb{Q}$ en calculant d'abord ceux de $Y \exp(-Y)$ par les méthodes de la section précédente, puis en inversant ce développement par la méthode présentée ci-dessus.

Pour déterminer le comportement asymptotique des coefficients ainsi calculés lorsque $N$ tend vers l'infini (autrement dit, l'asymptotique du nombre d'arbres de taille $N$), une méthode utilisée en analyse d'algorithmes passe par le calcul du développement de $T$ au voisinage du point singulier $\exp(-1) = 1/e$, où $T$ vaut 1. En posant $u = \sqrt{2\sqrt{1 - eX}}$, on voit qu'il s'agit alors de calculer le développement de $y$ solution de

$$ \sqrt{2 - 2y \exp(1 - y)} = u, $$
au voisinage de $y = 1$ et $u = 0$, équation à laquelle la méthode présentée dans cette section s'applique pour donner le résultat en $O(M(N))$ opérations.

### Le N-ième terme

Bien qu'en général le coût de l'inverse compositionnel soit dominé par celui de la composition et donc en $O(\sqrt{N \log N} M(N))$ (voir la section suivante), il est possible d'obtenir le N-ième terme sans calculer les précédents pour $O(M(N))$ opérations. L'idée repose sur la formule d'inversion de Lagrange :

$$ [X^N]F^{(-1)}(X) = \frac{1}{N}[X^{N-1}]\frac{1}{(F(X)/X)^N} \quad (3.5) $$

où la notation $[X^k]F$ représente le coefficient de $X^k$ dans la série F. Ainsi, le calcul par cette formule ne requiert qu'une puissance et donc peut être effectué en $O(M(N))$ opérations.

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

Cette formule donne immédiatement une jolie expression pour les coefficients de T solution de $T = X \exp(T)$ de l'exemple précédent : le n-ième coefficient vaut $n^{n-1}/n!$. Cette expression permet de calculer les N premiers coefficients en $O(N)$ opérations, ce qui est légèrement mieux que le $O(M(N))$ du cas général. L'asymptotique mentionnée plus haut redonne dans cet exemple la formule de Stirling.

## 3.4 La composition des séries

Si $F(X)$ et $G(X)$ sont des séries, avec $G(0) = 0$, il s'agit de calculer efficacement la série $F(G(X))$. À la différence des algorithmes vus dans les sections précédentes, il s'avère que les algorithmes connus n'atteignent pas une complexité quasi-optimale (c'est-à-dire linéaire en $N$, à des facteurs logarithmiques près).

### Méthode naïve

La composition peut être calculée par la méthode de Horner. Si $F(X) = \sum_{i \ge 0} f_i X^i + O(X^N)$, le calcul utilise alors

$$ F(G(X)) = f_0 + G \times (f_1 + G \times (f_2 + \dots)) + O(G^N) $$

et $O(G^N) = O(X^N)$.

**Exercice 3.7** Montrer que cette formule amène à un coût $O(NM(N))$.

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

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

Une technique dite « pas de bébés, pas de géants » réduit ce coût à $O(\sqrt{N}(M(N) + MM(\sqrt{N}))) = O(\sqrt{N}(M(N) + N^{\theta/2}))$, où $\theta$ est un exposant réalisable pour la complexité du produit de matrices, présenté au Chapitre 8. La série F est d'abord écrite

$$ F(X) = F_0(X) + X^k F_1(X) + \dots + X^{k[\lceil N/k \rceil]} F_{\lceil N/k \rceil}(X) + O(X^N), $$
où les polynômes $F_i$, en nombre $1 + \lceil N/k \rceil$, ont un degré au plus $k-1$. Ensuite, on calcule les puissances $G^2, \dots, G^k$ de G en $kM(N)$ opérations. Soient $f_{i,j}$ les coefficients des $F_i$ et $g_{i,j}$ ceux des $G^i$. Le développement

$$F_i(G) = \sum_{j=0}^{N-1} \left( \sum_{\ell=0}^{k-1} f_{i,\ell} g_{\ell,j} \right) X^j + O(X^N)$$

montre que les coefficients des $F_i(G)$ sont les coefficients du produit de la matrice ($f_{i,j}$) de taille ($\lceil N/k \rceil + 1) \times k$ par la matrice ($g_{i,j}$) de taille $k \times N$. En découpant chacune de ces matrices en blocs de taille $k \times k$ — la première devient un vecteur colonne d'environ $N/k^2$ blocs, la seconde un vecteur ligne d'environ $N/k$ blocs —, ce produit est effectué en au plus $O(N^2k^{\theta-3})$ opérations. Ensuite, il ne reste plus qu'à effectuer $\lceil N/k \rceil$ produits à précision N suivis d'autant d'additions pour un coût total de $O(NM(N)/k)$. Le coût total est minimisé par le choix $k = \lceil \sqrt{N} \rceil$ qui mène à la complexité annoncée.

### L'algorithme de Brent-Kung

Bien qu'il n'y ait pas d'algorithme permettant d'abaisser la complexité de la composition au même ordre que $M(N)$, il est possible de faire sensiblement mieux que la méthode précédente, en évitant la multiplication matricielle, grâce à l'Algorithme 3.8 sophistiqué, dû à Brent et Kung.

**Théorème 3.15** Si $2, 3, \dots, \lceil \sqrt{N \log N} \rceil$ sont inversibles dans $\mathbb{A}$ et si $F(X)$ et $G(X)$ sont deux séries de $\mathbb{A}[[X]]$ avec $G(0) = 0$ et $G'(0)$ inversible, alors l'Algorithme 3.8 calcule la série $F(G(X)) \mod X^N$ en $O(\sqrt{N \log NM(N)})$ opérations arithmétiques, ou en $O(\sqrt{NM(N)})$ opérations arithmétiques si $\log M(N)/\log N \to \gamma > 1$.

L'idée de départ de l'algorithme permettant d'aboutir à cette complexité est d'utiliser la formule de développement de Taylor, sous la forme

$$F(G_1 + X^m G_2) = F(G_1) + F'(G_1)X^m G_2 + F''(G_1) \frac{X^{2m} G_2^2}{2!} + \dots, \quad (3.6)$$

où $G_1$ est un polynôme de degré inférieur à $m$; le choix de $m$ est ajusté plus loin pour minimiser la complexité. Le premier gain de complexité par rapport à la méthode naïve provient de l'observation suivante.

**Lemme 3.16** Étant données les séries $G$ et $F \circ G$, avec $G'(0)$ inversible, la série $F' \circ G$ peut être calculée en $O(M(N))$ opérations arithmétiques.

*Démonstration.* La dérivation de séries a une complexité linéaire, et la division est en $O(M(N))$. Le résultat provient alors de la formule de dérivation d'une composée (déjà utilisée en page 71) :

$$F' \circ G = \frac{(F \circ G)'}{G'}.$$

L'hypothèse sur $G'(0)$ permet de garantir l'inversion du dénominateur.

<figure><img src="image_14.png" /></figure>
**Entrée** Deux séries tronquées F(X) mod XN et G(X) mod XN
avec G(0) = 0 et G'(0) inversible.

**Sortie** Leur composition F(G(X)) mod XN.

1. Soit $m = \lfloor \sqrt{N} \rfloor$ si la multiplication utilisée emploie la FFT,
$m = \lfloor \sqrt{N}/\log \lfloor \sqrt{N} \rfloor \rfloor$ sinon.

2. Découper G en $G_1 = G \mod X^m$ et $G_2 = G - G_1$, poser $H := 1$.

3. Calculer $U = 1/G'_1$ mod XN.

4. Calculer $R = S_0 := F(G_1) \mod X^N$ par la méthode récursive suivante :

a. Si $\text{deg} F = 1$ renvoyer $F(0) + (F - F(0))G_1 \mod X^N$,

b. sinon

i. soit $k = \lceil \text{deg}(F)/2 \rceil$,

ii. découper F en A := F mod X^k et B := (F - A)/X^k,

iii. calculer récursivement $G_1^{\lceil k/2 \rceil}$, A(G₁), et B(G₁) modulo XN,

iv. renvoyer A(G₁) + G₁^k B(G₁) mod XN.

5. Pour $i = 1, \dots, \lceil N/m \rceil$ faire

a. $H := HG_2 \mod X^N$,

b. $S_i := S'_{i-1} U/i \mod X^N$,

c. $R := R + S_i H \mod X^N$.

6. Renvoyer R.

Algorithme 3.8 – Algorithme de composition de Brent-Kung.

L'utilisation répétée de cette idée dans la formule (3.6) mène à une complexité

$$C(F(G_1)) + \frac{N}{m}O(M(N))$$

pour l'ensemble de la composition, où $C(F(G_1))$ représente la complexité du calcul de la première composition avec un polynôme, donnée par le lemme suivant.

**Lemme 3.17** Soient F et G deux polynômes de degrés k et m, avec $G(0) = 0$. Alors, le calcul des N premiers coefficients de la série F ◦ G peut être effectué en $O(km \log NM(N)/N)$ opérations arithmétiques, ou en $O(kmM(N)/N)$ opérations si $\log M(N)/\log N \to \gamma > 1$.

*Démonstration du théorème à l'aide de ce lemme.* Ce lemme fournit l'estimation

$$C(F(G_1)) = O(mM(N) \log N).$$

La complexité totale est donc bornée par

$$O(mM(N) \log N) + \frac{N}{m}O(M(N)) = O\left(M(N)\left(m \log N + \frac{N}{m}\right)\right).$$

Cette dernière somme est minimisée par le choix de $m = \sqrt{N/\log N}$ qui donne le résultat du théorème. Le cas sans le facteur $\log N$ est laissé en exercice. ■
*Démonstration du lemme.* Comme d'habitude, on commence par traiter le cas où le degré *k* est une puissance de 2. Ensuite, quitte à considérer que les coefficients de F pour les degrés allant de *k* + 1 jusqu'à la puissance de 2 immédiatement supérieure sont nuls, la perte de complexité est d'au plus un facteur 2, absorbé dans la constante du O(·).

L'idée est d'effectuer ce calcul par « diviser pour régner » en récrivant le poly-
nôme F sous la forme

$$
F = F_1 + X^{k/2}F_2,
$$

où $F_1$ et $F_2$ sont deux polynômes de degré au plus $k/2$. Dans la composition

$$
F(G) = F_{1}(G) + G^{k/2}F_{2}(G),
$$

les deux polynômes du second sommant ont degré au plus km/2. La complexité C_k^m
découlant de l'utilisation de cette formule vérifie donc

$$
C_k^m \le 2C_{k/2}^m + 2M(\min(N, km/2)) + O(N).
$$

C'est dans cette prise de minimum que réside toute la subtilité : tant que *k* est grand,
les calculs sont tronqués à l'ordre N, et dès que *k* devient suffisamment petit, on
exploite l'arithmétique polynomiale. Il vient donc

$$
\begin{align*}
C_k^m \le 2M(N) + 4M(N) + \dots + 2^{\ell-1}M(N) + 2^\ell \left( M\left(\frac{km}{2^\ell}\right) + 2M\left(\frac{km}{2^{\ell+1}}\right) + \dots \right), \\
\le 2^\ell M(N) + 2^\ell M\left(\frac{km}{2^\ell}\right) \log\left(\frac{km}{2^\ell}\right),
\end{align*}
$$

où la valeur de $\ell$, qui conclut la preuve du lemme, est déterminé par

$$
\frac{km}{2^l} < N \le \frac{km}{2^{l-1}}.
$$

Exercices

Exercice 3.8 — Exponentielle de série. Soit $s$ une série de terme constant nul et $n = 2^p$. On étudie l'itération définie par $y_0 = z_0 = 1$ et pour $k > 0$

$$
\begin{align*}
z_k &= z_{k-1} + z_{k-1} (1 - y_{k-1} z_{k-1}) \mod X^{2^{k-1}}, \\
y_k &= y_{k-1} - y_{k-1} \int z_k (y'_{k-1} - s' y_{k-1}) \mod X^{2^k}.
\end{align*}
$$

1. Montrer que $y_p - \exp(s) = O(X^n)$;

2. Estimer la complexité du calcul des *n* premiers coefficients de la série exp(*s*) par cette itération;

3. (Plus difficile) Interpréter cette itération comme l'itération de Newton associée à l'opérateur $Y \mapsto Y' - s'Y$.
**Exercice 3.9 — Composition avec l'exponentielle.** Soit $f(X) \in \mathbb{Q}[[X]]$ une série à coefficients rationnels. Le but de cet exercice est de montrer que les N premiers coefficients de la série $S(X) = f(e^X - 1)$ peuvent être calculés en $O(M(N) \log N)$ opérations arithmétiques dans $\mathbb{Q}$.

1. Montrer que $S(X) - A(e^X - 1) = O(X^N)$, où $A(X)$ est l'unique polynôme de degré inférieur à N tel que $f(X) = A(X) + O(X^N)$.

2. On pose $B(X) = A(X-1)$; donner un algorithme pour calculer les coefficients de B à partir de ceux de A en $O(M(N) \log N)$ opérations dans $\mathbb{Q}$.

3. La transformée de Laplace formelle $\mathcal{L}$ est définie sur $\mathbb{Q}[[X]]$ comme l'application $\mathbb{Q}$-linéaire telle que $\mathcal{L}(X^k) = k! X^k, k \in \mathbb{N}$. Si $B(X) = b_0 + b_1 X + \dots + b_{N-1} X^{N-1}$, montrer que

$$ \mathcal{L}(B(e^X)) = \sum_{i=0}^{N-1} \frac{b_i}{1 - iX} $$

4. Donner un algorithme pour calculer les N premiers coefficients de $B(e^X)$ à partir de ceux de $B(X)$ en $O(M(N) \log N)$ opérations dans $\mathbb{Q}$.

5. Donner finalement un algorithme calculant les N premiers coefficients de S à partir de ceux de f en $O(M(N) \log N)$ opérations dans $\mathbb{Q}$.

■

**Exercice 3.10 — Composition de séries avec arcsinus.** Soit $\mathbb{K}$ un corps de caractéristique nulle, et soit $F \in \mathbb{K}[X]$ une série formelle de terme constant nul. Écrire un algorithme à base d'itérations de Newton permettant de calculer les N premiers termes de la série composée $(\arcsin \circ F)$, à partir des N premiers termes de la série F, en $O(M(N))$ opérations arithmétiques dans $\mathbb{K}$.

■

**Exercice 3.11 — Composition de séries formelles en petite caractéristique.** Soit $\mathbb{K} = \mathbb{F}_p$ le corps fini à $p$ éléments, où $p$ est un nombre premier. On rappelle que tout élément $a$ de $\mathbb{K}$ vérifie l'égalité $a^p = a$. Le but de cet exercice est de concevoir un algorithme pour la composition des séries dans $\mathbb{K}[[X]]$, de complexité asymptotique quasi-optimale en la précision. Soient $f$ et $g$ deux séries dans $\mathbb{K}[[X]]$ tronquées à précision N, avec $g(0) = 0$, et soit $h$ leur composée $h = f \circ g \bmod X^N$. On suppose $N \ge p$.

1. Montrer que $g(X)^p \bmod X^N$ s'écrit $s(X^p)$ pour une certaine série $s \in \mathbb{K}[[X]]$ connue à précision $\lfloor N/p \rfloor$.

2. Montrer que le calcul de $h$ se ramène à $p$ compositions de séries à précision $\lfloor N/p \rfloor$.

3. En déduire un algorithme récursif pour le calcul de $h$, et estimer sa complexité.

■

## Notes

La méthode de Newton remonte à 1669, et se trouve bien décrite dans son traité des fluxions [New40]. Cette méthode a ensuite été raffinée et étendue de nombreuses façons [Yam00; Ypm95].

L'utilisation de l'itération de Newton en informatique est très ancienne. L'inversion
de matrice de la Section 3.2 remonte à Schulz [Sch33] pour des matrices réelles. L'importance de la méthode de Newton pour le calcul rapide avec des séries formelles a été soulignée par Lipson [Lip76] : « In a symbolic mathematics system setting, we feel that Newton's method offers a unified and easily implemented algorithm for solving a wide variety of power series manipulation problems that can be cast into a root-finding mold. Our experience attests to the practicality of Newton's method. We feel that it is a great algebraic algorithm. »

Du point de vue de la complexité, Sieveking [Sie72] a donné l'Algorithme 3.2 d'inversion rapide de séries. Kung [Kun74] a observé qu'il s'agissait d'une itération de Newton et a amélioré la constante dans la complexité. Ensuite, Brent [Bre76] a montré comment l'exponentielle et le logarithme s'obtiennent également en O(M(N)) opérations. Les algorithmes de la Section 3.3 sont dus à Schönhage [Sch82]; l'application aux sommes et produits composés en Section 3.3 est beaucoup plus récente [Bos+06b]. La méthode de Newton permet également la résolution de systèmes linéaires à coefficients polynomiaux (Chapitre 11) et la résolution d'équations ou de systèmes différents (Chapitre 13). L'idée d'incorporer le dividende dans la dernière itération de Newton, afin de gagner un facteur constant pour la division de séries (Section 3.2), est due à Karp et Markstein [KM97].

Les algorithmes non triviaux de composition décrits en Section 3.4 sont dus à Brent et Kung [BK78]. Une version optimisée est décrite dans l'article de Johansson [Joh15], qui utilise la formule d'inversion de Lagrange (3.5). L'idée de l'algorithme « pas de bébés, pas de géants » en Section 3.4 remonte à Paterson et Stockmeyer [PS73]. La composition $F \circ G$ mod $X^N$ peut s'effectuer en coût quasi-optimal $O(M(N) \log N)$ si $G$ est un polynôme [BK78], ou plus généralement lorsque $G$ est une série algébrique [Hoe02c]. Des compositions avec d'autres fonctions particulières peuvent aussi être calculées en $O(M(N) \log N)$, ou parfois $O(M(N))$ [BSS08]. Savoir si cela est encore possible pour des séries quelconques est un grand problème ouvert. La formule d'inversion évoquée en page 72 a été découverte par Lagrange en 1768 [Lag68].

Dans le cas où la caractéristique $p$ de $A$ est finie, Bernstein [Ber98] a montré qu'il est possible de descendre la composition à une complexité quasi-optimale en la précision $N$ des séries, pour le prix d'une dépendance linéaire en $p$. Cet algorithme est intéressant lorsque la caractéristique $p$ est petite. L'exercice 3.11 présente l'idée de l'algorithme dans un cas particulier. Une avancée récente importante est due à Kedlaya et Umans [KU08], qui ont donné le premier algorithme de complexité binaire quasi-linéaire à la fois en $N$ et en $\log p$, pour le calcul de $F \circ G \mod X^N$, lorsque $F$ et $G$ sont des séries à coefficients dans le corps fini $\mathbb{F}_p$.

Ritzmann [Rit86] a montré que la composition des séries entières peut s'effectuer en complexité binaire quasi-optimale : si $F, G \in \mathbb{Z}[X]$ ont des coefficients bornés par $\ell$ en valeur absolue, alors il est possible de calculer $F \circ G \mod X^N$ en $O(M(\mathbb{Z}[N^2 \log N \log \ell]))$ opérations binaires. Il montre aussi que le problème de la composition de séries à coefficients dans un corps $K$ est d'une difficulté équivalente à celui du calcul du $N$-ième coefficient de $F \circ G$. Le problème de l'inverse compositionnel leur est également équivalent [BK78].

En présence d'une structure supplémentaire, il est possible de calculer la composition ou l'inverse compositionnel en complexité arithmétique optimale par rapport à la précision. C'est le cas par exemple pour l'inverse compositionnel d'un polynôme,
ou plus généralement d'une série algébrique. C'est également le cas pour la composition $F \circ G$ d'une série différentiellement finie F par une série algébrique G. Dans tous ces cas de figure, la série calculée est différentiellement finie, et les techniques des chapitres sur ces séries (Corollaire 14.10 et Théorème 15.1) s'appliquent : les N premiers termes peuvent se calculer en $O(N)$ opérations arithmétiques, et le N-ième terme en un nombre d'opérations arithmétiques quasi-linéaire en $\sqrt{N}$. Il n'est actuellement pas connu si une complexité linéaire en N peut être obtenue pour d'autres types de structures, comme par exemple pour le calcul de l'inverse compositionnel $F^{(-1)}(X) \mod X^N$ d'une série différentiellement finie F, ou de la composition $F \circ G \mod X^N$ de deux séries différentiellement finies F et G.

En pratique, les constantes cachées dans les $O(\cdot)$ jouent un grand rôle, et nécessitent généralement l'implantation à la fois des algorithmes asymptotiquement meilleurs et des algorithmes plus classiques, souvent plus efficaces en petite taille. Dans le modèle de multiplication polynomiale par FFT, les meilleurs constantes devant $M(N)$ sont dues à Harvey [Har10; Har11] : $\frac{13}{9}$ pour l'inverse, $\frac{4}{3}$ pour la racine carrée, $\frac{13}{6}$ pour l'exponentielle. Ces problèmes ont aussi été étudiés par van der Hoeven [Hoe10]. Des articles de synthèse dus à Bernstein [Berb; Ber08] décrivent diverses techniques utilisées pour éliminer certains calculs redondants dans la méthode de Newton.

L'Exercice 3.8 est inspiré par [HZ], et l'Exercice 3.9 est tiré de [BSS08].

D'autres modèles, appelés *paresseux* et *détendus*, permettent de manipuler les séries de façon plus naturelle sans se soucier de la précision, qui est gérée de façon automatique. Chaque série est vue comme le vecteur des coefficients déjà calculés et une fonction permettant de calculer un coefficient de degré N en fonction des précédents [BHL11; Hoe02c; Hoe03; Hoe07].

Bibliographie

Berb BERNSTEIN, Daniel J. (2001). *Removing redundancy in high-precision Newton iteration*. URL: http://cr.yp.to/fastnewton.html (visible en 2001).

Ber08 — (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.

Ber98 — (1998). « Composing power series over a finite ring in essentially linear time ». In : *Journal of Symbolic Computation*, vol. 26, n°3, p. 339–341.

BHL11 BERTHOMIEU, Jérémy, Joris van der HOEVEN et Grégoire LECERF (2011). « Relaxed algorithms for *p*-adic numbers ». In : *Journal de Théorie des Nombres de Bordeaux*, vol. 23, n°3, p. 541–577.

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+06b BOSTAN, Alin, Philippe FLAJOLET, Bruno SALVY et Éric SCHOST (2006). « Fast computation of special resultants ». In : *Journal of Symbolic Computation*, vol. 41, n°1, p. 1–29.
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.

BSS08 BOSTAN, Alin, Bruno SALVY et Éric SCHOST (2008). « Power series composition and change of basis ». In : *ISSAC'08 : International Symposium on Symbolic and Algebraic Computation*. Éd. par David J. JEFFREY. ACM Press, p. 269–276.

Har10 HARVEY, David (2010). « Faster exponentials of power series ». url : https://arxiv.org/abs/0911.3110.

Har11 — (2011). « Faster algorithms for the square root and reciprocal of power series ». In : *Mathematics of Computation*, vol. 80, p. 387–394.

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.

Hoe03 — (2003). « Relaxed multiplication using the middle product ». In : *ISSAC'03 : International Symposium on Symbolic and Algebraic Computation*. Éd. par Manuel BRONSTEIN. Philadelphia, USA : ACM Press, p. 143–147.

Hoe07 — (2007). « New algorithms for relaxed multiplication ». In : *Journal of Symbolic Computation*, vol. 42, n°8, p. 792–802.

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

HZ HANROT, Guillaume et Paul ZIMMERMANN (2004). *Newton iteration revisited*. url : http://www.loria.fr/~zimmerma/papers (visible en 2004).

Joh15 JOHANSSON, Fredrik (2015). « A fast algorithm for reversion of power series ». In : *Mathematics of Computation*, vol. 84, n°291, p. 475–484.

KM97 KARP, A. H. et P. MARKSTEIN (1997). « High-precision division and square root ». In : *ACM Transactions on Mathematical Software*, vol. 23, n°4, p. 561–589.

KU08 KEDLAYA, Kiran S. et Christopher UMANAS (2008). « Fast modular composition in any characteristic ». In : *FOCS'08 : IEEE Conference on Foundations of Computer Science*. Washington, DC, USA : IEEE Computer Society, p. 146–155.

Kun74 KUNG, H. T. (1974). « On computing reciprocals of power series ». In : *Numerische Mathematik*, vol. 22, p. 341–348.

Lag68 LAGRANGE, Joseph Louis, comte de (1768). « Nouvelle méthode pour résoudre les équations littérales par le moyen des séries ». In : *Mémoires de l'Académie Royale des Sciences et Belles-Lettres de Berlin*, vol. 24, p. 251–326.

Lip76 LIPSON, J. D. (1976). « Newton's method : a great algebraic algorithm ». In : *SYMSAC'76 : ACM Symposium on Symbolic and Algebraic Computation*. ACM Press, p. 260–270.

New40 NEWTON, Isaac (1740). *La méthode des fluxions, et les suites infinies*. Traduction française par G. Buffon du texte de 1671 de Newton. de Bure. url : http://www.gallica.fr.
PS73 PATERSON, M. S. et L. J. STOCKMEYER (1973). « On the number of nonscalar multiplications necessary to evaluate polynomials ». In : *SIAM Journal on Computing*, vol. 2, n°1, p. 60–66.

Rit86 RITZMANN, P. (1986). « A fast numerical algorithm for the composition of power series with complex coefficients ». In : *Theoretical Computer Science*, vol. 44, p. 1–16.

Sch33 SCHULZ, G. (1933). « Iterative Berechnung der reziproken Matrix ». In : *Zeitschrift für angewandte Mathematik und Physik*, vol. 13, p. 57–59.

Sch82 SCHÖNHAGE, A. (1982). *The fundamental theorem of algebra in terms of computational complexity*. Preliminary Report. 73 pages. Tübingen, Germany : Mathematisches Institut der Universität.

Sie72 SIEVEKING, M. (1972). « An algorithm for division of powerseries ». In : *Computing*, vol. 10, p. 153–156.

Yam00 YAMAMOTO, Tetsuro (2000). « Historical developments in convergence analysis for Newton's and Newton-like methods ». In : *Journal of Computational and Applied Mathematics*, vol. 124, n°1-2, p. 1–23.

Ypm95 YPMA, Tjalling J. (1995). « Historical development of the Newton–Raphson method ». In : *SIAM Review*, vol. 37, n°4, p. 531–551.
# 4. Division euclidienne, fractions rationnelles et récurrences linéaires à coefficients constants

## Résumé

La multiplication et l'inversion rapide de séries permettent le calcul efficace de la division euclidienne de polynômes et du développement en série des fractions rationnelles. Ils mènent en particulier au calcul rapide d'un ou de plusieurs termes d'une suite récurrente linéaire à coefficients constants.

## 4.1 Introduction

Les suites récurrentes linéaires à coefficients constants sont extrêmement présentes dans la nature. On estime ainsi qu'environ 12% des suites de la littérature<fnref n="1" /> sont de ce type. L'exemple le plus célèbre est celui de la suite de Fibonacci, définie par $F_0 = 0$, $F_1 = 1$ et

$$F_{n+2} = F_{n+1} + F_n, \quad n \ge 0. \qquad (4.1)$$

L'utilisation directe de cette récurrence permet de calculer les N premiers éléments de la suite de manière optimale du point de vue de la complexité. En revanche, pour le calcul d'un seul terme d'indice élevé, le passage par le calcul de tous ses prédécesseurs est inutilement coûteux. Une première idée de ce chapitre consiste à utiliser une version matricielle des récurrences pour aboutir à une complexité quasi-optimale pour cette question.

La récurrence (4.1) n'est que d'ordre 2. Pour des suites définies par des récurrences d'ordre $d$ plus grand, il faut aussi chercher à limiter le coût des calculs vis-à-vis de $d$. Simplement dérouler la récurrence n'est plus alors la meilleure façon de calculer les N premiers termes de la suite. Une autre idée importante de ce chapitre est d'exploiter le lien entre les récurrences linéaires à coefficients constants et les fractions

<footnote n="1">1. Recensées très largement par l'encyclopédie en ligne des suites d'entiers, http://oeis.org.</footnote>
rationnelles, qui mène à une complexité quasi-optimale vis-à-vis à la fois de N et d pour cette question, ainsi que pour celle du calcul des N premiers termes, et aussi pour le développement en série des fractions rationnelles.

Les applications algorithmiques de ces idées sont nombreuses. Certaines d'entre elles sont présentées dans des exercices de ce chapitre, comme le calcul rapide d'un polynôme sur 1, 2, ..., N, une partie du calcul du polynôme minimal d'une matrice, et des tests probabilistes de primalité d'entiers.

## 4.2 Division de polynômes

Étant donnés deux polynômes F et G de $A[X]$, avec G unitaire, à coefficients dans un anneau $A$, il existe d'uniques polynômes Q et R, quotient et reste de la division euclidienne de F par G, tels que :

$$F = QG + R \quad \text{avec} \quad \deg R < \deg G.$$

Calculer rapidement le quotient Q et le reste R est d'importance vitale dans toute la suite de cet ouvrage. Outre l'application aux suites récurrentes qui sera détaillée plus loin, les algorithmes de division sont utilisés au Chapitre 5 pour l'évaluation multipoint et l'interpolation et au Chapitre 6 pour le calcul de pgcd. À leur tour, ces deux algorithmes sont centraux dans nombre d'applications.

### Méthode naïve

L'algorithme naïf pour calculer Q et R consiste à poser la division. Pour cela, les termes dominants sont d'abord isolés :

$$F = aX^m + T_F, \quad G = X^n + T_G, \quad \text{avec } \deg T_F < m, \deg T_G < n.$$

Si $m < n$, il n'y a rien à faire. Sinon, la boucle élémentaire de la division de F par G consiste à « tuer » le terme de plus haut degré de F en effectuant la soustraction

$$F - aX^{\delta}G = (aX^{m} + T_{F}) - (aX^{\delta+n} + aX^{\delta}T_{G}) = T_{F} - aX^{\delta}T_{G},$$

où $\delta = m - n$. Le polynôme obtenu est congru à F modulo G, mais de degré strictement inférieur à $m$. Il n'y a plus qu'à itérer ce procédé jusqu'à obtenir un polynôme de degré strictement inférieur à $n$, et à garder trace des quotients successifs.

La complexité de cet algorithme est facile à estimer. La boucle élémentaire présentée ci-dessus s'effectue en $O(n)$ opérations de type $(+, -, \times)$; dans le pire des cas, l'algorithme effectue $(m-n+1)$ passages dans cette boucle, de sorte que la complexité de la division de F par G est de $O(n(m-n))$ opérations. Le pire des cas est atteint pour $m = 2n$. C'est donc un algorithme *quadratique* (mais un bon algorithme quadratique : la constante dans le $O(\cdot)$ est en fait petite, c'est 2). Le même genre d'estimation s'obtient pour la division euclidienne classique des entiers.

### Algorithme rapide

Un bien meilleur résultat peut être obtenu à base d'itération de Newton. L'idée principale de l'algorithme part de la récriture de $F = QG + R$ en

$$\frac{F}{G} = Q + \frac{R}{G}.$$