www.pi314.net


Histoire
Mathématiciens
Toutes les formules
Approx. numériques
Programmes
Algos perso/divers
Décimales
Poèmes
Articles/vidéos
Délires !
 Pi-Day
Images/Fonds
Musique
Liens sur Pi
Bibliographie



Boris Gourévitch
L'univers de Pi - V2.57
modif. 13/04/2013



Pi-Day dans
Pi Day Countdown
Google
Accueil Historique/Actu (Pi, site, moi) Edito Livre d'or Pages en .pdf Je me présente Quelques photos Remerciements Page des nets d'or Sites qui m'indexent Derniers changements Contact

Cette page en français This page in English


pictpict

J.W. Cooley et J.W. Tukey

Multiplication rapide par FFT



1 Retour à l’école

Si l’on vous demande de multiplier 32*14 dans votre tête, comment allez-vous faire ? Sans doute, tout en vous rappelant douloureusement des règles enseignées a l’école primaire, utiliser 32=3*10+2, 14=10+4, puis multiplier membre à membre, gérer les retenues et tout ça. Ouf, ça parait le plus simple... Mais est-ce vraiment le plus simple ? Question bizarre, n’est-ce pas ?

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

pict

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

e2i2knp= wk,    w = e22ipn
(3)

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

pict

et      2
y = z  . On obtient alors

P(z) = P1(y)+ z.P2(y).
(6)

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

       (      )
(wk)2 = w 22n+k 2 = (wn+k)2 .
(7)

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

F(2n) = 2F (n)+ 2 2n, F(1) = 0
(8)

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

                          (        )
                           logn
F(2n) = 2  2n (d+ 1) = 4n log2 + 1
(9)

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

W  (r0,r1,...,r2n-1) = (a0,a1,...,a2n-1)
(10)

    (                              )
       1    1      12    ...    12n-1
       1    w2    w4        w 4n-2
W =    1   w      w         w
       ...                      ...
       1  w2n-1  w4n-2      w2n2- n
(11)

ce qui est équivalent à

                    -1
(r0,r1,...,r2n-1) = W   (a0,a1,...,a2n-1)
(12)

  -1
W  a la bonne idée d’être simple

          ( 1     1       1    ...    1    )
            1    w-1     w-2        w-2n+1
  -1   -1-  1    w-2     w-4        w-4n+2
W    = 2n   .                          .     .
            ..                          ..2
            1  w- 2n+1  w- 4n+2      w -2n+n
(13)

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.


Retour à la page d'accueil