Arithmétique multi-précision

grands entiers ;
réels représentés avec une grande précision ;
arithmétique modulaire avec un grand modulo, éléments de
grands corps finis ;
polynômes de grand degré (à coefficients dans Z, Q, un corps
fini, etc.).
Entiers
Pn et polynômes : cas très semblables, car manipuler
i
i=0 i β (représentation en base β, typiquement une grande
a
puissance de 2) ou
Xn
ai X i
i=0

n’est pas très différent, l’arithmétique polynomiale étant souvent

plus simple (pas de retenues).
Taille : de quelques centaines à plusieurs milliards de bits

analyse numérique, lorsque les formats VF usuels ne suffisent

pas ; vérification et mise au point d’algorithmes (ex.
approximations polynomiales, valeurs tabulées, calcul de points
difficiles à arrondir) : centaines de bits
Taille : de quelques centaines à plusieurs milliards de bits

analyse numérique, lorsque les formats VF usuels ne suffisent

pas ; vérification et mise au point d’algorithmes (ex.
approximations polynomiales, valeurs tabulées, calcul de points
difficiles à arrondir) : centaines de bits
cryptographie : RSA, Diffie-Hellman. Besoin d’entiers de
grande taille pour résister aux attaques : millers de bits
Taille : de quelques centaines à plusieurs milliards de bits

analyse numérique, lorsque les formats VF usuels ne suffisent

pas ; vérification et mise au point d’algorithmes (ex.
approximations polynomiales, valeurs tabulées, calcul de points
difficiles à arrondir) : centaines de bits
cryptographie : RSA, Diffie-Hellman. Besoin d’entiers de
grande taille pour résister aux attaques : millers de bits
“mathématiques expérimentales” : recherches de propriétés, de
contre-exemples, etc. : millions de bits et au-delà.
Á quoi ça sert ?

D’abord, c’est beau. . . mais parfois ça sert :

Physique :
Frolov : calculs à précision > 100 chiffres décimaux en 2003
pour le calcul d’interactions électromagnétiques ;
Laskar (Observatoire Paris) : stabilité à long terme (> 108
années) du système solaire.
Á quoi ça sert ?

D’abord, c’est beau. . . mais parfois ça sert :

Physique :
Frolov : calculs à précision > 100 chiffres décimaux en 2003
pour le calcul d’interactions électromagnétiques ;
Laskar (Observatoire Paris) : stabilité à long terme (> 108
années) du système solaire.
Théorie des nombres :
certains problèmes nécessitent la manipulation de polynômes
de degré ≈ 100000 ;
plus grand nombre premier explicitement connu

224036583 − 1

(environ 7 millions de chiffres décimaux), prouvé en utilisant le

test de Lucas-Lehmer ;
Á quoi ça sert ?

découverte “expérimentale” de la formule de

Bailey-Borwein-Plouffe
∞  
X 1 4 2 1 1
π= − − − ,
16k 8k + 1 8k + 4 8k + 5 8k + 6
k=0

qui permet d’avoir “directement” le jème chiffre binaire de π.

Á quoi ça sert ?

découverte “expérimentale” de la formule de

Bailey-Borwein-Plouffe
∞  
X 1 4 2 1 1
π= − − − ,
16k 8k + 1 8k + 4 8k + 5 8k + 6
k=0

qui permet d’avoir “directement” le jème chiffre binaire de π.

contre-exemples à des conjectures, ou renforcement de la
confiance qu’on leur accorde.
Á quoi ça sert ?

découverte “expérimentale” de la formule de

Bailey-Borwein-Plouffe
∞  
X 1 4 2 1 1
π= − − − ,
16k 8k + 1 8k + 4 8k + 5 8k + 6
k=0

qui permet d’avoir “directement” le jème chiffre binaire de π.

contre-exemples à des conjectures, ou renforcement de la
confiance qu’on leur accorde. P
Exemple (théorème de Robin). Notons σ(n) = d|n d .
L’hypothèse de Riemann est vraie si et seulement si

σ(n)
∀n > 5040, < eγ ,
n ln ln n
où γ est la constante d’Euler.
Multiplication entière et réelle : algorithmes simples

La multiplication est l’opération la plus importante en

multi-précision : la division et la racine carrées sont construites à
partir de la multiplication (et sans coûter beaucoup plus cher, grâce
à la méthode de Newton).
base β : entier assez grand (les “chiffres” sont des entiers
représentables sur le processeur utilisé, appelés “petits
entiers”) ;
d’abord : multiplication et division d’un grand entier par un
petit entier.
Multiplication par un “petit” entier, ou “scalaire”
Produit a × b, où b est multiprécision, et a est “atomique”

b = bn−1 β n−1 + bn−2 β n−2 + · · · + b0

Multiplication par un “petit” entier, ou “scalaire”
Produit a × b, où b est multiprécision, et a est “atomique”

b = bn−1 β n−1 + bn−2 β n−2 + · · · + b0

et on suppose 0 ≤ a, bi ≤ β − 1. L’algorithme est élémentaire :

c0 ← 0
for i = 0 to n − 1 do
d ← a · bi + ci
pi ← d mod β
ci+1 ← d div β
end for
pn ← cn
Multiplication par un “petit” entier, ou “scalaire”
Produit a × b, où b est multiprécision, et a est “atomique”

b = bn−1 β n−1 + bn−2 β n−2 + · · · + b0

et on suppose 0 ≤ a, bi ≤ β − 1. L’algorithme est élémentaire :

c0 ← 0
for i = 0 to n − 1 do
d ← a · bi + ci
pi ← d mod β
ci+1 ← d div β
end for
pn ← cn
On montre par récurrence l’invariant :

(ci pi−1 pi−2 · · · p0 )β = a × (bi−1 bi−2 · · · b0 )β

et que les pi et les ci sont ≤ β − 1.

Division par un “petit” entier
Divisons b par a, où a et b sont écrits comme précédemment. On
obtient :
qn−1 ← bn−1 div a
rn−1 ← bn−1 mod a
for i = n − 2 to 0 do
ti ← βri+1 + bi
qi ← ti div a
ri ← ti mod a
end for
r ← r0
Division par un “petit” entier
Divisons b par a, où a et b sont écrits comme précédemment. On
obtient :
qn−1 ← bn−1 div a
rn−1 ← bn−1 mod a
for i = n − 2 to 0 do
ti ← βri+1 + bi
qi ← ti div a
ri ← ti mod a
end for
r ← r0
et on montre par récurrence sur p que

a × β p−1 qn−1 + β p−2 qn−2 + · · · + qn−p + rn−p

 

= β p−1 bn−1 + β p−2 bn−2 + · · · + bn−p

et que 0 ≤ rn−p ≤ a.
Multiplication naïve

multiplier deux nombres multi-précision de même taille n, a et

b
se généralise sans problème à des tailles différentes
la méthode naïve consiste à multiplier a par chacun des
chiffres bi à l’aide de la “multiplication scalaire” et à accumuler
les résultats obtenus, avec le décalage adéquat
Temps O(n2 ) (que Kolmogorov pensait optimal).
Multiplication de Karatsuba
A et B sur 2k “chiffres” → on les décompose en A = A1 β k + A0 et
B = B1 β k + B0 , où A1 , A0 , B1 et B0 s’écrivent sur k chiffres.
Multiplication de Karatsuba
A et B sur 2k “chiffres” → on les décompose en A = A1 β k + A0 et
B = B1 β k + B0 , où A1 , A0 , B1 et B0 s’écrivent sur k chiffres.
On définit C = (A1 − A0 )(B1 − B0 ), on obtient facilement :

AB = β 2k A1 B1 + β k (A1 B1 + A0 B0 − C ) + A0 B0 ,
Multiplication de Karatsuba
A et B sur 2k “chiffres” → on les décompose en A = A1 β k + A0 et
B = B1 β k + B0 , où A1 , A0 , B1 et B0 s’écrivent sur k chiffres.
On définit C = (A1 − A0 )(B1 − B0 ), on obtient facilement :

AB = β 2k A1 B1 + β k (A1 B1 + A0 B0 − C ) + A0 B0 ,

→ calcule AB en n’effectuant que 3 multiplications de

nombres de k chiffres : A1 B1 , A0 B0 et (A1 − A0 )(B1 − B0 ).
Multiplication de Karatsuba
A et B sur 2k “chiffres” → on les décompose en A = A1 β k + A0 et
B = B1 β k + B0 , où A1 , A0 , B1 et B0 s’écrivent sur k chiffres.
On définit C = (A1 − A0 )(B1 − B0 ), on obtient facilement :

AB = β 2k A1 B1 + β k (A1 B1 + A0 B0 − C ) + A0 B0 ,

→ calcule AB en n’effectuant que 3 multiplications de

nombres de k chiffres : A1 B1 , A0 B0 et (A1 − A0 )(B1 − B0 ).
Chaque fois que l’on double la taille n des nombres manipulés, on
triple le temps de calcul Kar (n). On trouve facilement
 
Kar (n) = O nlog(3)/ log(2) .
Multiplication de Karatsuba
A et B sur 2k “chiffres” → on les décompose en A = A1 β k + A0 et
B = B1 β k + B0 , où A1 , A0 , B1 et B0 s’écrivent sur k chiffres.
On définit C = (A1 − A0 )(B1 − B0 ), on obtient facilement :

AB = β 2k A1 B1 + β k (A1 B1 + A0 B0 − C ) + A0 B0 ,

→ calcule AB en n’effectuant que 3 multiplications de

nombres de k chiffres : A1 B1 , A0 B0 et (A1 − A0 )(B1 − B0 ).
Chaque fois que l’on double la taille n des nombres manipulés, on
triple le temps de calcul Kar (n). On trouve facilement
 
Kar (n) = O nlog(3)/ log(2) .
√
log(3)/ log(2) ≈ 1.585 → on gagne presque un facteur n par
rapport à la méthode naïve.
Multiplication de Karatsuba

En pratique devient intéressant dès que la taille des nombres

manipulés dépasse quelques dizaines de mots machine.
S’adapte trivialement pour effectuer le produit de 2 nombres
complexes en ne faisant que 3 multiplications réelles.
Digression : produit matriciel de Strassen

La méthode de Karatsuba a très certainement inspiré Strassen

pour la mise au point de son algorithme de produit de matrices
n × n en temps nlog(7)/ log(2) .
Digression : produit matriciel de Strassen

La méthode de Karatsuba a très certainement inspiré Strassen

pour la mise au point de son algorithme de produit de matrices
n × n en temps nlog(7)/ log(2) .
Supposons que l’on veuille multiplier deux matrices d’ordre 2n,
A et B, décomposées en blocs de taille n comme suit :
   
A1,1 A1,2 B1,1 B1,2
A= , et B = ,
A2,1 A2,2 B2,1 B2,2

pour en calculer le produit

 
C1,1 C1,2
C= .
C2,1 C2,2
si on note :


 M1 = (A1,1 + A2,2 )(B1,1 + B2,2 )



 M2 = (A2,1 + A2,2 )B1,1
 M3 = A1,1 (B1,2 − B2,2 )


M4 = A2,2 (B2,1 − B1,1 ) ,
M5 = (A1,1 + A1,2 )B2,2




 M6 = (A2,1 − A1,1 )(B1,1 + B1,2 )



M7 (A1,2 − A2,2 )(B2,1 + B2,2 )

=
si on note :


 M1 = (A1,1 + A2,2 )(B1,1 + B2,2 )



 M2 = (A2,1 + A2,2 )B1,1
 M3 = A1,1 (B1,2 − B2,2 )


M4 = A2,2 (B2,1 − B1,1 ) ,
M5 = (A1,1 + A1,2 )B2,2




 M6 = (A2,1 − A1,1 )(B1,1 + B1,2 )



M7 (A1,2 − A2,2 )(B2,1 + B2,2 )

=

alors 

 C1,1 = M1 + M4 − M5 + M7
C1,2 = M3 + M5

.
 C2,1
 = M2 + M4
C2,2 = M1 − M2 + M3 + M6

Multiplication de Toom-Cook

Andrei Toom, 1963 — Stephen Cook, 1966 ;

Pour simplifier, on supposera être en base 2, mais se généralise
sans peine dans une autre base.
On ramène le produit d’entiers au produit de polynômes. Idée :
entiers → polynômes : un entier a de n bits est découpé en k
parties de longueur `
→ a = Ak−1 2`(k−1) + Ak−2 2`(k−2) + · · · + A0 , avec
0 ≤ Ai ≤ 2` − 1. Le nombre a est donc la valeur en 2` du
polynôme

A(X ) = Ak−1 X k−1 + Ak−2 X k−2 + · · · + A0 .

Multiplication de Toom-Cook

a et b → polynômes A et B. Si on connait la valeur de A(X )

et B(X ) en 2k − 1 points choisis de sorte que ces valeurs
soient simples à calculer et soient des entiers  a et b, en
multipliant ces valeurs 2 à 2, on connait la valeur du polynôme
R(X ) = A(X )B(X ) en ces 2k − 1 points ;
Multiplication de Toom-Cook

a et b → polynômes A et B. Si on connait la valeur de A(X )

et B(X ) en 2k − 1 points choisis de sorte que ces valeurs
soient simples à calculer et soient des entiers  a et b, en
multipliant ces valeurs 2 à 2, on connait la valeur du polynôme
R(X ) = A(X )B(X ) en ces 2k − 1 points ;
R est de degré 2k − 2 → connaitre sa valeur en 2k − 1 points
suffit pour le “reconstruire” (i.e., connaître ses coefficients—là
encore, les points doivent être choisis pour que ce soit facile) ;
Multiplication de Toom-Cook

a et b → polynômes A et B. Si on connait la valeur de A(X )

et B(X ) en 2k − 1 points choisis de sorte que ces valeurs
soient simples à calculer et soient des entiers  a et b, en
multipliant ces valeurs 2 à 2, on connait la valeur du polynôme
R(X ) = A(X )B(X ) en ces 2k − 1 points ;
R est de degré 2k − 2 → connaitre sa valeur en 2k − 1 points
suffit pour le “reconstruire” (i.e., connaître ses coefficients—là
encore, les points doivent être choisis pour que ce soit facile) ;
la calcul de AB(2` ), qui ne pose aucune difficulté particulière
donne le résultat.
Multiplication de Toom-Cook

a et b → polynômes A et B. Si on connait la valeur de A(X )

et B(X ) en 2k − 1 points choisis de sorte que ces valeurs
soient simples à calculer et soient des entiers  a et b, en
multipliant ces valeurs 2 à 2, on connait la valeur du polynôme
R(X ) = A(X )B(X ) en ces 2k − 1 points ;
R est de degré 2k − 2 → connaitre sa valeur en 2k − 1 points
suffit pour le “reconstruire” (i.e., connaître ses coefficients—là
encore, les points doivent être choisis pour que ce soit facile) ;
la calcul de AB(2` ), qui ne pose aucune difficulté particulière
donne le résultat.
Ce mécanisme peut s’utiliser une fois ou de manière récursive (pour
calculer AB aux 2k − 1 points choisis).
Multiplication de Toom-Cook

Les points choisis en pratique sont de tous petits entiers

(. . . , −2, −1, 0, 1, 2, . . .), ou des inverses de toutes petites
puissances de 2 (1/2, 1/4) ainsi que +∞ (abus de langage :
on considère le terme de tête du polynôme).
Par exemple, GMP utilise ∞, 2, −1, 1, 0. L’idée qui sous-tend
ceci est que les valeurs de A et B en ces points sont d’une
taille juste légèrement plus grande que `.
avec des petites √
modifications de l’algorithme, des valeurs
comme ±i ou ± 2 seraient certainements intéressantes à
utiliser.
Multiplication de Toom-Cook

Ici : cas le plus utile en pratique : k = 3 (algorithme

“Toom-3”).
Polynômes A et B : A2 X 2 + A1 X + A0 et B2 X 2 + B1 X + B0 ,
où A2 , A1 , A0 , B2 , B1 etB0 s’écrivent sur au plus ` = dn/3e
chiffres.
Multiplication de Toom-Cook

Ici : cas le plus utile en pratique : k = 3 (algorithme

“Toom-3”).
Polynômes A et B : A2 X 2 + A1 X + A0 et B2 X 2 + B1 X + B0 ,
où A2 , A1 , A0 , B2 , B1 etB0 s’écrivent sur au plus ` = dn/3e
chiffres.
Nous devons choisir 2k − 1 = 5 points. Nous choisirons ∞, 2,
−1, 1, 0. La valeur de A en ces points est donc :
en ∞ : A(∞) = A2 ;
en 2 : A(2) = 4A2 + 2A1 + A0 ;
en −1 : A(−1) = A2 − A1 + A0 ;
en 1 : A(1) = A2 + A1 + A0 ;
en 0 : A(0) = A0 .
Multiplication de Toom-Cook

Ici : cas le plus utile en pratique : k = 3 (algorithme

“Toom-3”).
Polynômes A et B : A2 X 2 + A1 X + A0 et B2 X 2 + B1 X + B0 ,
où A2 , A1 , A0 , B2 , B1 etB0 s’écrivent sur au plus ` = dn/3e
chiffres.
Nous devons choisir 2k − 1 = 5 points. Nous choisirons ∞, 2,
−1, 1, 0. La valeur de A en ces points est donc :
en ∞ : A(∞) = A2 ;
en 2 : A(2) = 4A2 + 2A1 + A0 ;
en −1 : A(−1) = A2 − A1 + A0 ;
en 1 : A(1) = A2 + A1 + A0 ;
en 0 : A(0) = A0 .
même chose pour B.
Multiplication de Toom-Cook

En base 2, calculer les termes A(i) ou B(i) (i = ∞, 2, ±1, 0) :

additions et décalages ;
Multiplication de Toom-Cook

En base 2, calculer les termes A(i) ou B(i) (i = ∞, 2, ±1, 0) :

additions et décalages ;
la plus grande valeur que puissent avoir ces termes est
(4 + 2 + 1)(2` − 1) = 3 · 2`+1 − 6 → ils sont donc
représentables sur au plus ` + 3 bits.
Multiplication de Toom-Cook

En base 2, calculer les termes A(i) ou B(i) (i = ∞, 2, ±1, 0) :

additions et décalages ;
la plus grande valeur que puissent avoir ces termes est
(4 + 2 + 1)(2` − 1) = 3 · 2`+1 − 6 → ils sont donc
représentables sur au plus ` + 3 bits.
le calcul des produits A(i) · B(i) se fait donc par des
multiplications de nombres d’au plus dn/3e + 3 bits.
Multiplication de Toom-Cook

Passons maintenant au calcul des coefficients du polynôme

P = AB, à partir des 5 valeurs P(i) = A(i) · B(i). Ce polynôme est
de degré 4, et si on appelle Pk son coefficient de degré k, on a bien
évidemment :
en ∞ : P(∞) = P4 ;
en 2 : P(2) = 16P4 + 8P3 + 4P2 + 2P1 + P0 ;
en −1 : P(−1) = P4 − P3 + P2 − P1 + P0 ;
en 1 : P(1) = P4 + P3 + P2 + P1 + P0 ;
en 0 : P(0) = P0 .
Multiplication de Toom-Cook

Passons maintenant au calcul des coefficients du polynôme

P = AB, à partir des 5 valeurs P(i) = A(i) · B(i). Ce polynôme est
de degré 4, et si on appelle Pk son coefficient de degré k, on a bien
évidemment :
en ∞ : P(∞) = P4 ;
en 2 : P(2) = 16P4 + 8P3 + 4P2 + 2P1 + P0 ;
en −1 : P(−1) = P4 − P3 + P2 − P1 + P0 ;
en 1 : P(1) = P4 + P3 + P2 + P1 + P0 ;
en 0 : P(0) = P0 .
On connaît les P(i), on cherche les Pk .
Multiplication de Toom-Cook
c’est un système linéaire de matrice
 
1 0 0 0 0
 
 16 8 4 21 
 
 
 1 −1
M= 1 −1 1 
.
 
 1 1 1 1 1 
 
0 0 0 0 1
Cette matrice est inversible (Vandermonde), et
 
1 0 0 0 0
 
 −2 1/6 −1/6 −1/2 1/2 
 
M−1
 
 −1
= 0 1/2 1/2 −1 .
 
 2 −1/6 −1/3 1 −1/2 
 
0 0 0 0 1
Multiplication de Toom-Cook

On a donc



 P4 = P(∞)


P = −2P(∞) + 16 P(2) − 16 P(−1) − 12 P(1) + 21 P(0)

 3



P2 = −P(∞) + 12 P(−1) + 12 P(1) − P(0)


2P(∞) − 16 P(2) − 13 P(−1) + P(1) − 21 P(0)




 P1 =


 P = P(0)
0

divisions par 2 : des décalages.

divisions par 3 ou 6 : divisions par petits entiers (linéaires—pas
besoin de nombreux chiffres fractionnaires du quotient : on sait
que le résultat final pi est un entier !).
Il est possible de montrer qu’on ne peut pas éviter ces divisions.
Multiplication de Toom-Cook

produit de 2 nombres de n bits : 5 produits de nombres de

dn/3e + 3 bits. En oubliant lâchement le d.e et le +3, on
trouve alors que la complexité est en
 
O nlog 5/ log 3 = O n1.46··· .

Le cas général de l’algorithme de Toom-Cook (découpage en k

entiers de longueur `), appelé parfois “Toom-k” donnera une
complexité en  
O nlog(2k−1)/ log k .
Multiplication de Toom-Cook

Complexité en
 
O nlog 5/ log 3 = O n1.46··· .


Multiplication de Toom-Cook

Complexité en
 
O nlog 5/ log 3 = O n1.46··· .

on peut donc pour tout  > 0 trouver un algorithme de

multiplication en 0(n1+ ).
Multiplication de Toom-Cook

Complexité en
 
O nlog 5/ log 3 = O n1.46··· .

on peut donc pour tout  > 0 trouver un algorithme de

multiplication en 0(n1+ ).
la “constante cachée” croît très rapidement avec k, ce qui fait
que seuls Toom-3 et Toom-5 ont un intérêt pratique
Multiplication de Toom-Cook

Complexité en
 
O nlog 5/ log 3 = O n1.46··· .

on peut donc pour tout  > 0 trouver un algorithme de

multiplication en 0(n1+ ).
la “constante cachée” croît très rapidement avec k, ce qui fait
que seuls Toom-3 et Toom-5 ont un intérêt pratique
Raison : on n’a pas un nombre arbitrairement grand de points
d’interpolation “simples” (i.e., pour lesquels les matrices A et
A−1 sont toutes les deux très simples).
Multiplication de Toom-Cook

Complexité en
 
O nlog 5/ log 3 = O n1.46··· .

on peut donc pour tout  > 0 trouver un algorithme de

multiplication en 0(n1+ ).
la “constante cachée” croît très rapidement avec k, ce qui fait
que seuls Toom-3 et Toom-5 ont un intérêt pratique
Raison : on n’a pas un nombre arbitrairement grand de points
d’interpolation “simples” (i.e., pour lesquels les matrices A et
A−1 sont toutes les deux très simples).
Faire mieux ? prendre comme points les racines ne de l’unité (dans
C ou dans un anneau) → FFT.
La transformée de Fourier discrète

très utilisée en traitement du signal

Au départ : séries de Fourier, représentation d’une fonction
périodique f — ici de périodeP2π — comme somme d’une
+∞
série trigonométrique : a0 + n=1 (an cos nx + bn sin nx)
A soulevé de nombreuses questions intéressantes.
premières idées : Bernoulli (cordes vibrantes) ;
Fourier (travaillant sur l’équation de la chaleur) donne les
relations
1
R 2π
a0 = 2π 0
f (t)dt
2π
ap = π1 0 f (t) cos(pt)dt
R
R 2π
bp = π1 0 f (t) sin(pt)dt
Plusieurs questions se sont posées : pour quelles fonctions f la
série obtenue converge-t-elle, et a-t-elle pour somme f ?
Problème étudié par Dirichlet, Riemann, Lebesgue. . .
La transformée de Fourier discrète dans C
2i π
Soit ω = e − , et définissons la matrice
n

 
1 1 1 ··· 1
 1 ω ω 2 ··· ω n−1 
 
 1 ω2 ω 4 ··· ω 2(n−1) 
 
 .. .. .. .. .. 
 . .
Fω =  . . . 

 1 ωi ω 2i ··· ω i(n−1) 
 
 .. .. .. .. .. 
 . . . . . 
1 ω n−1 ω 2(n−1) · · · ω (n−1)(n−1)

i.e., Fω;ij = ω (i−1)(j−1) , alors la TFD d’ordre n est l’application

x ∈ Cn → X ∈ Cn définie par

X = Fω x.

On note parfois X = Fn (x).

TFD inverse

···
 
1 1 1 1

 1 ω −1 ω −2 ··· ω −(n−1) 


 1 ω −2 ω −4 ··· ω −2(n−1) 

 .. .. .. .. .. 
Fω−1 =
 . . . . . 


 1 ω −i ω −2i ··· ω −i(n−1) 

 .. .. .. .. .. 
 . . . . . 
1 ω −(n−1) ω −2(n−1) · · · ω −(n−1)(n−1)
vérifie

Fω Fω−1 = nIn .
On a donc :
1
(Fω )−1 = F −1 .
n ω
TFD inverse

D’où la définition de la transformée de Fourier inverse :

n−1
1X
x` = Xk (ω −1 )k` .
n
k=0

On notera x = Fn−1 (X ).
Algorithme FFT (Fast Fourier Transform) de Cooley-Tukey

probablement inventé vers 1805 par Gauss ;

Algorithme FFT (Fast Fourier Transform) de Cooley-Tukey

probablement inventé vers 1805 par Gauss ;

redécouvert (puis oublié) en 1914 ;
Algorithme FFT (Fast Fourier Transform) de Cooley-Tukey

probablement inventé vers 1805 par Gauss ;

redécouvert (puis oublié) en 1914 ;
formulation moderne : article de 1965 de Cooley et Tukey
Algorithme FFT (Fast Fourier Transform) de Cooley-Tukey

probablement inventé vers 1805 par Gauss ;

redécouvert (puis oublié) en 1914 ;
formulation moderne : article de 1965 de Cooley et Tukey
TFD, en supposant que n est une puissance de 2 (sinon, ∃
algorithmes basés sur factorisation de n).
Algorithme FFT (Fast Fourier Transform) de Cooley-Tukey

probablement inventé vers 1805 par Gauss ;

redécouvert (puis oublié) en 1914 ;
formulation moderne : article de 1965 de Cooley et Tukey
TFD, en supposant que n est une puissance de 2 (sinon, ∃
algorithmes basés sur factorisation de n).
On veut calculer, pour k = 0, . . . , n − 1,
n−1
X
Xk = x` ω k` .
`=0

Idée : se ramener à des TFD d’ordre n/2.

Algorithme FFT (Fast Fourier Transform) de Cooley-Tukey

n est une puissance de 2,

ω 2 est une racine (n/2)ième de l’unité
On s’intéresse aux TFD d’ordre n/2 (avec ω 2 comme racine de
l’unité) des vecteurs

(x0 , x2 , x4 , . . . , xn−2 ) → (P0 , P1 , . . . , Pn/2−1 )

et
(x1 , x3 , x5 , . . . , xn−1 ) → (I0 , I1 , . . . , In/2−1 ).
On a :
Pn−1 k`
Xk = `=0 x` ω
Pn/2−1
x2m ω k·(2m) + x2m+1 ω k·(2m+1)


= m=0
Pn/2−1 2 km + ω k
Pn/2−1 2 km
= m=0 x2m (ω ) m=0 x2m+1 (ω ) .
Algorithme FFT (Fast Fourier Transform) de Cooley-Tukey

pour k ≤ n/2 − 1, on retombe immédiatement sur la FFT

d’ordre n/2 :
Xk = Pk + ω k Ik

pour k compris entre n/2 et n − 1, notons k 0 = k − n2 . Il vient

0 +n/2)m 0m
(ω 2 )km = (ω 2 )(k = (ω 2 )k
0+ n 0
et ω k = ω k 2 = −ω k , d’où l’on tire

Xk = Pk−n/2 − ω k−n/2 Ik−n/2 .

On ramène donc le calcul d’une TFD d’ordre n à celui de 2 TFD

d’ordre n/2. L’algorithme FFT est constitué par ce calcul récursif.
Algorithme FFT (Fast Fourier Transform) de Cooley-Tukey

on suppose que les ω i (pour i < n) sont précalculés,

la FFT d’ordre 2 consomme 2 opérations (ω n/2 = −1 est
racine carrée de l’unité, donc X0 = x0 + x1 et X1 = x0 − x1 ),
la FFT d’ordre n demande 2 FFT d’ordre n/2 et 2n opérations
dans C.
Conclusion : Par récurrence, la FFT d’ordre n consomme
2n log2 n − n = O(n log n) opérations dans C.
Algorithme FFT (Fast Fourier Transform) de Cooley-Tukey

on suppose que les ω i (pour i < n) sont précalculés,

la FFT d’ordre 2 consomme 2 opérations (ω n/2 = −1 est
racine carrée de l’unité, donc X0 = x0 + x1 et X1 = x0 − x1 ),
la FFT d’ordre n demande 2 FFT d’ordre n/2 et 2n opérations
dans C.
Conclusion : Par récurrence, la FFT d’ordre n consomme
2n log2 n − n = O(n log n) opérations dans C.
S’adapte trivialement à la TFD inverse (ω → ω −1 , et division par
n).
Remarque importante
de
   
1 1 1 ··· 1 x0

 1 ω ω2 ··· ω n−1  
  x1 


 1 ω2 ω4 ··· ω 2(n−1)  
  x2 

 .. .. .. .. ..   .. 
Fω · x = 
 . . . . . ·
  . ,


 1 ωi ω 2i ··· ω i(n−1)  
  xi 

 .. .. .. .. ..   .. 
 . . . . .   . 
1 ω n−1 ω 2(n−1) · · · ω (n−1)(n−1) xn−1

on tire que Xk (i.e., la k + 1ième composante de F .x) est la valeur

du polynôme

Px (t) = x0 + x1 t + x2 t 2 + · · · + xn−1 t n−1

au point ω k .
Remarque importante

FFT : algorithme de calcul de la valeur d’un polynôme aux

points 1, ω, ω 2 , . . . , ω n−1 en ne consommant que O(n log n)
opérations dans C.
Remarque importante

FFT : algorithme de calcul de la valeur d’un polynôme aux

points 1, ω, ω 2 , . . . , ω n−1 en ne consommant que O(n log n)
opérations dans C.
FFT inverse : moyen rapide d’interpoler les coefficients d’un
polynôme de degré n − 1 à partir de sa valeur aux points
1, ω, ω 2 , . . . , ω n−1 .
Remarque importante

FFT : algorithme de calcul de la valeur d’un polynôme aux

points 1, ω, ω 2 , . . . , ω n−1 en ne consommant que O(n log n)
opérations dans C.
FFT inverse : moyen rapide d’interpoler les coefficients d’un
polynôme de degré n − 1 à partir de sa valeur aux points
1, ω, ω 2 , . . . , ω n−1 .
→ même principe évaluation-interpolation que l’algorithme de
Toom-Cook.
Application à la multiplication

Idée générale d’algorithme pour multiplier deux grands entiers.

on veut multiplier deux entiers (représentés en base β) :

x = xp−1 β p−1 + xp−2 β p−2 + · · · + x0

et
y = yp−1 β p−1 + yp−2 β p−2 + · · · + y0
on leur associe les deux polynômes

x(t) = xp−1 t p−1 + xp−2 t p−2 + · · · + x0

et
y (t) = yp−1 t p−1 + yp−2 t p−2 + · · · + y0
Application à la multiplication

on construit des vecteurs de dimension 2p, en complétant les

coefficients de x(t) et y (t) avec des zéros
(0, · · · , 0 xp−1 , xp−2 , . . . , x0 ) et (0, · · · , 0 yp−1 , yp−2 , . . . , y0 )
| {z } | {z }
p termes p termes
on calcule par l’algorithme FFT la TFD d’ordre 2p de ces deux
vecteurs, on obtient (X2p−1 , . . . , X0 ) et (Y2p−1 , . . . , Y0 )
on calcule les produits termes à termes des Xi et des Yi , i.e.,
on calcule
Ci = Xi · Yi
pour i = 0, 1, . . . , 2p − 1. Suivant la taille de Xi et Yi pourra
se faire récursivement ou par un produit de deux petits entiers.
Application à la multiplication

Coefficients Ci : valeurs du polynôme x(t)y (t) en les racines

2pièmes de l’unité.
on calcule la TFD inverse d’ordre 2p de
(C2p−1 , C2p−2 , · · · , C0 ), qui donne un vecteur
(c2p−1 , c2p−2 , . . . , c0 ) ;
le polynôme x(t)y (t) vaut c2p−1 t 2p−1 + c2p−2 t 2p−2 + · · · + c0 ,
et l’entier xy vaut c2p−1 β 2p−1 + c2p−2 β 2p−2 + · · · + c0
Application à la multiplication

Attention :
Il ne faut pas en déduire que la complexité de la multiplication
d’entiers est en O(n log n), ce sera un peu plus, car dans C, les
calculs devront être menés avec une précision qui dépendra de
n.
Pour que ceci soit intéressant, il faut que les produits Xi · Yi
n’aient besoin de se faire qu’avec une précision très petite
devant la taille de x et y (on verra que c’est de l’ordre du
logarithme de cette précision)
Algorithme de Schönhage et Strassen (il y en a 2)

calculs approchés en utilisant des nombres VF ;

 : erreur relative d’une opération complexe ;
pour simplifier, on suppose que tous les arrondis se font dans
le sens qui minimise le module du résultat.
En partant de (xp−1 , xp−2 , . . . , x0 ) et (yp−1 yp−2 , . . . , y0 ) où
|xi |, |yi | ≤ β − 1 sont les chiffres des entiers que l’on veut
multiplier, on va calculer comme expliqué plus haut les coefficients
(c2p−1 , c2p−2 , . . . , c0 ). Ils sont calculés de manière approchée mais
ce sont des entiers → si on s’arrange pour que l’erreur absolue sur
chaque ci soit < 1/2, on les retrouvera sans difficulté
Algorithme de Schönhage et Strassen (il y en a 2)

Une opération élémentaire de l’algorithme FFT :

opération exacte : a ← b + ρc, où ρ est une puissance de ω ;
calcul effectivement réalisé : a0 ← ◦(b0 + ρ0 c 0 ), où ◦ désigne
l’arrondi.
On a |ρ| = 1, et on suppose que les valeurs précalculées des ρ0 sont
telles que |ρ0 | ≤ 1 (arrondis dans le “bon sens” lors du précalcul), et
avec une erreur relative ≤ , i.e.,

ρ0 − ρ
≤ .
ρ

On note ` = log2 (2p) = log2 (p) + 1).

Algorithme de Schönhage et Strassen (il y en a 2)
L’étude d’erreur n’est pas difficile, mais fastidieuse. À la fin, on
trouve une erreur finale (sur les ci ) majorée par

6 · ` · 22` ·  · (β − 1)2 .

On veut que cette erreur soit < 1/2, ce qui donne la condition
1
< (1)
12 · ` · 22`· (β − 1)2
Avoir une erreur relative  demande à manipuler des nombres d’au
moins 1 − log2 () bits (vous verrez partout − log2 (), mais il est
nécessaire d’avoir un bit de plus pour arrondir d’une façon qui
minimise le module : on n’arrondit pas forcément au plus près). Il
faut donc que l’on représente les nombres sur

2 log(β − 1) + 2` + log ` + log 12 + 1 bits

Algorithme de Schönhage et Strassen (il y en a 2)
Donc on se ramène de la manipulation de nombres de

N = p log2 (β) bits

à celle de nombres de (on rappelle : ` = log2 (p) + 1))

2 log(β − 1) + 2 log2 p + 2 + log2 (log2 (p) + 1) + log2 3 + 3 bits

Par exemple, avec N = 220 ≈ 106 bits, et p = 32768 (donc

β = 2N/p = 232 ), on se retrouve à manipuler des nombres de 104
bits.
De manière générale, on se retrouve donc de la manipulation de
nombres de N bits à celle de nombres d’environ
N
2 + 2 log2 (p) bits.
p
Algorithme de Schönhage et Strassen (il y en a 2)
Si on ne fait qu’une seule étape de l’algorithme, en faisant les
multiplications élémentaires des “petits” nombres avec l’algorithme
 2
en O(n2 ), on fait O(p log p) opérations de coût O Np + log p ,
ce qui donne un coût total (en temps) en
 2 !
N
O p · log p · + log p .
p
Algorithme de Schönhage et Strassen (il y en a 2)
Si on ne fait qu’une seule étape de l’algorithme, en faisant les
multiplications élémentaires des “petits” nombres avec l’algorithme
 2
en O(n2 ), on fait O(p log p) opérations de coût O Np + log p ,
ce qui donne un coût total (en temps) en
 2 !
N
O p · log p · + log p .
p

En prenant p = N/ log2 (N) pour équilibrer les termes (en pratique,

recherche systématique de l’optimum autour de cette valeur), on
obtient un coût en temps en
 
O N log2 N .
Algorithme de Schönhage et Strassen (il y en a 2)
Si on ne fait qu’une seule étape de l’algorithme, en faisant les
multiplications élémentaires des “petits” nombres avec l’algorithme
 2
en O(n2 ), on fait O(p log p) opérations de coût O Np + log p ,
ce qui donne un coût total (en temps) en
 2 !
N
O p · log p · + log p .
p

En prenant p = N/ log2 (N) pour équilibrer les termes (en pratique,

recherche systématique de l’optimum autour de cette valeur), on
obtient un coût en temps en
 
O N log2 N .

On s’est ramené d’une taille N à une taille log N : pas de cas réels
où ceci demande à nouveau une très grande précision → dans la
vraie vie cet algorithme ne s’utilise pas récursivement
Algorithme de Schönhage et Strassen (il y en a 2)
En théorie, si on réutilise cet algorithme récursivement pour calculer
les produits terme à terme Xi · Yi et ceux qui apparaissent dans les
calculs élémentaires a ← b + ρc de le FFT, on arrive à un coût en
∗
N log N log log N log log log N · · · × 2O(log N)
.

En effet, si τ (N) est le temps de calcul d’un produit de nombres de

N bits, alors en choisissant à nouveau de découper nos nombres en
p = N/ log2 N blocs, on effectue λ(p log p) opérations sur des
nombres de taille O(log N), par conséquent
 
N N
τ (N) = λ log × τ (log N) ≈ λ(N × τ (log N)),
log N log N
ce qui donne
∗
τ (N) = N log(N) log log(N) log log log(N) · · · × 2O(log N)
.
L’autre algorithme de Schönhage et Strassen

utilise la TFD dans l’anneau

k
Z/(22 + 1)Z;

complexité en
O(N log(N) log log(N));
L’autre algorithme de Schönhage et Strassen

utilise la TFD dans l’anneau

k
Z/(22 + 1)Z;

complexité en
O(N log(N) log log(N));

a été implanté. . . pour multiplier 2 nombres de 784141 mots de

64 bits, ≈ 1 seconde sur un opteron 4GHz ;
L’autre algorithme de Schönhage et Strassen

utilise la TFD dans l’anneau

k
Z/(22 + 1)Z;

complexité en
O(N log(N) log log(N));

a été implanté. . . pour multiplier 2 nombres de 784141 mots de

64 bits, ≈ 1 seconde sur un opteron 4GHz ;
a été l’algorithme de meilleure complexité connue jusqu’en
2007 ;
L’autre algorithme de Schönhage et Strassen

utilise la TFD dans l’anneau

k
Z/(22 + 1)Z;

complexité en
O(N log(N) log log(N));

a été implanté. . . pour multiplier 2 nombres de 784141 mots de

64 bits, ≈ 1 seconde sur un opteron 4GHz ;
a été l’algorithme de meilleure complexité connue jusqu’en
2007 ;
algorithme de Martin Fürer en
∗ (N))
N log(N) · 2O(log .