Multiplication rapide par FFT

Depuis l’avènement des machines à calculer, c’est à dire la fin des années 40, le manque d’algorithme de multiplication efficace oblige à segmenter chaque nombre en tranches, par exemple $20 10 B = B210 +B110 + B0$ comme à l’école. La multiplication de nombres de taille $n$ nécessite alors un temps (ou nombre d’opérations) proportionnel à $2 n$ sans compter l’utilisation mémoire proportionnelle à $n$ . Afin de conserver le rythme de progression du record de calcul de décimales de $p$ au cours des années 50 (voir pages Histoire ou Décimales), des améliorations théoriques et algorithmiques semblent nécessaires. C’est alors à cette époque que les choses s’accélèrent.

En 1965, Cooley and Tukey introduisent sous sa forme moderne une méthode de réduction de la complexité du calcul des séries de Fourier connue aujourd’hui sous le nom de Transformée de Fourier Rapide [1]. Il semble que Gauss ait déjà deviné l’astuce de la factorisation critique dès 1805. Quel génie, celui-là... Quoi qu’il en soit, Schönhage et Strassen en tirent un algorithme de multiplication de grands entiers en compexité $O (n.log(n).log(log(n)))$ ce qui est considérablement mieux que $( 2) O n$ [2] !

Cette page explique le fonctionnement de cet algorithme de multiplication rapide de deux grands entiers par FFT, qui a tant fait pour le calcul des décimales de $p$ ! Pour les ordis, c’est bien... Mais pour nous, finalement, euh... ce n’est pas si naturel !

2 L’algorithme de multiplication rapide de deux grands entiers par FFT

Soient X et Y deux grands entiers de taille $n$ en base $B = 10$ usuellement (ou une puissance de $10$ ).

2.1 Décomposition polynomiale

On écrit $X$ et $Y$ sous forme polynomiale en les décomposant (de manière unique) dans la base $B$

où $x j$ et $y j$ sont donc respectivement les $j$ -ème décimales de $X$ et $Y$ ( $0 < x ,y < B - 1 j j$ ). Cette opération est de complexité à peu près proportionnelle à $n$ . On cherche maintenant à effectuer “rapidement” la multiplication $R(z) = P(z)× Q(z)$ puis à évaluer à la fin $R(B)$ . $R(z)$ est polynôme de degré $< 2n$ et peut être retrouvé par interpolation à partir de son évaluation en $2n$ points, autrement dit par l’évaluation de $P$ et $Q$ en ces mêmes $2n$ points. Le choix des points, c’est la finesse de l’algorithme ! On utilise les racines primitives $2n$ -ième de l’unité, i.e. les

car cette évaluation (qui n’est alors rien d’autre que le calcul d’une série de Fourier) peut être réalisée en complexité $O(n.logn)$ si $n$ est une puissance de $2$ et avec un peu d’astuce. C’est la Transformée de Fourier Rapide ou FFT en anglais.

2.2 Transformée de Fourier rapide

En conservant les notations précédentes, évaluons le polynôme $P (wk)$ et $Q (wk)$ en se plaçant sous l’hypothèse que $n = 2d$ . Posons

Maintenant, il faut se rappeler que puisque $w$ est une racine $2n$ -ième de l’unité, alors pour $k (- {1,...,2n}$

En conséquence, évaluer $P$ aux racines $2n$ -ièmes de l’unité revient à évaluer $P1$ et $P2$ chacun aux $n$ points $( ) ( ) ( ) w2 1, w2 2,..., w2 n$ et à reporter les résultats dans $P$ à l’aide de l’équation 6. Si $F (2n)$ est le nombre d’opérations élémentaires (addition, multiplication) requises pour évaluer $P$ en $w1,w2,...,w2n$ alors en vertu du principe précédent et de 6, on a

où le terme $2 × 2n$ provient des dernières addition et multiplication requises pour obtenir $P$ à partir de $P1$ et $P2$ dans 6.

Comme $2 w$ est une racine primitive $n$ -ième de l’unité, on peut réappliquer le même procédé pour chaque polynôme $P1$ et $P2$ . Par suite, comme $d n = 2$ est une puissance de $2$ , le processus (et donc 6) s’itère encore $d$ fois. On obtient finalement un nombre d’opérations

d’où la complexité en $O(n.log n)$ de la FFT. De même pour la FFT inverse, qui utilise $w- k$ au lieu de $wk$ .

2.3 Interpolation

Nous avons maintenant calculé nos $( k) P w$ et $( k) Q w$ pour $k = 0,...,2n- 1$ . On forme les $2n$ produits $( k) ( k) ( k) R w = P w Q w$ qui reviennent à des multiplications élémentaires puisque $( k) P w$ est forcément de module inférieur à $sum n -1 j=0 xj < n(B - 1)$ d’après 1, de même pour $Q$ . Comme l’on cherche à retrouver $R(z) = sum 2nj=-0 1rjzj$ , il faut l’interpoler à partir des $2n$ valeurs $( ) ak = R wk$ . Autrement dit, on cherche à résoudre le système

Le calcul de $(r0,r1,...,r2n-1)$ n’est donc rien d’autre que le conjugué de la transformée de Fourier de $(a0,a1,...,a2n-1)$ , autrement dit sa transformée de Fourier inverse, dont on a vu qu’elle est en $O(n.log n)$ comme la FFT. On obtient alors $R(z)$ d’où l’on tire $R(B) = XY$ . Cette dernière opération est à peu près en complexité $O(n)$ car les $r0,r1,...,r2n-1$ sont déjà proches des décimales de $XY$ , aux retenues près (dès que $rj > B$ ). Toutes les retenues sont en fait calculées simultanément, puis reportées, puis on recalcule les éventuelles nouvelles retenues créées, etc... (opération appelée vectorisation).

La combinaison de tous ces algorithmes et quelques menus raffinements selon les programmeurs permettent d’atteindre à peu près la borne optimale théorique de Schönhage et Strassen pour la FFT en $O (n.log(n).log(log(n)))$ .

References

[1] J. W. Cooley, J. W. Tukey, “An algorithm for the machine calculation of complex Fourier series”, Math. Comp. 1965, vol 19, pp. 297-301.

[2] A. Schönhage, V. Strassen, “Schennelle Multiplikation grosser Zahlen”, Computing, 1971, vol. 7, pp. 281-292.