Lipsum.dev

Transformée de Fourier discrète et multiplications

27 décembre 2022 — FFT, Python

Dans cet article, nous allons introduire la transformée de Fourier discrète, qui associe à une séquence de nombres complexes une autre séquence de complexes de même longueur nn.

L’algorithme de transformation de Fourier rapide (FFT) permet de calculer cette transformée en O(nlogn)O(n \log n). En particulier, il peut être utilisé pour multiplier rapidement des polynômes, ou des entiers.

Le code source des fonctions présentées est disponible ici.

Transformée de Fourier discrète

La transformation de Fourier discrète associe à une séquence (a0,...,an1)(a_0, ..., a_{n-1}) de nombres complexes sa transformée (A0,...,An1)(A_0, ..., A_{n-1}).

Pour ce faire, on considère le polynôme P=a0+a1X+a2X2+...+an1Xn1{P = a_0 + a_1 X + a_2 X^2 + ... + a_{n-1} X^{n-1}} associé à la séquence initiale et on évalue ce polynôme en nn points bien choisis pour obtenir la transformée.

Plus précisément, on choisit les nn racines nn-èmes de l’unité 1,ω,ω2,...ωn1{1, \omega, \omega^2, ... \omega^{n-1}}.

On définit ici ω\omega comme une racine nn-ème primitive de l’unité, par exemple ω=e2iπn{\omega = e^{-\frac{2i\pi}{n}}}.

On a donc : Ak=P(ωk)=m=0n1ame2iπmkn{A_k = P(\omega^k) = \displaystyle\sum_{m=0}^{n-1} a_m e^{-2i\pi\frac{mk}{n}}}

Cette relation s’écrit aussi matriciellement :

(A0A1A2...An1)=(111...11ωω2...ωn11ω2ω22...ω2(n1)1ωn1ω2(n1)...ω(n1)(n1))(a0a1a2...an1)\begin{pmatrix} A_0 \\ A_1 \\ A_2 \\ ... \\ A_{n-1} \end{pmatrix} = \begin{pmatrix} 1 & 1 & 1 & ... & 1 \\ 1 & \omega & \omega^2 & ... & \omega^{n-1} \\ 1 & \omega^2 & \omega^{2 \cdot 2} & ... & \omega^{2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{n - 1} & \omega^{2(n-1)} & ... & \omega^{(n-1)(n-1)} \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ ... \\ a_{n-1} \end{pmatrix}

La transformation de Fourier inverse consiste à inverser cette matrice, je vous laisse vérifier que la matrice ci-dessous convient :

(a0a1a2...an1)=1n(111...11ω1ω2...ω(n1)1ω2ω22...ω2(n1)1ω(n1)ω2(n1)...ω(n1)(n1))(A0A1A2...An1)\begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ ... \\ a_{n-1} \end{pmatrix} = \frac{1}{n} \begin{pmatrix} 1 & 1 & 1 & ... & 1 \\ 1 & \omega^{-1} & \omega^{-2} & ... & \omega^{-(n-1)} \\ 1 & \omega^{-2} & \omega^{-2 \cdot 2} & ... & \omega^{-2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{-(n - 1)} & \omega^{-2(n-1)} & ... & \omega^{-(n-1)(n-1)} \end{pmatrix} \begin{pmatrix} A_0 \\ A_1 \\ A_2 \\ ... \\ A_{n-1} \end{pmatrix}

Par exemple, dans le cas n=2n = 2 la transformée de Fourier discrète de (a,b)(a, b) est (A,B)=(a+b,ab)(A, B) = (a + b, a - b) constituée de la somme et de la différence. La transformée inverse s’obtient par demi-somme et demi-différence : (a,b)=(A+B2,AB2){(a, b) = (\frac{A + B}{2}, \frac{A - B}{2})}.

Produit de polynômes

Si PP et QQ sont deux polynômes de degrés inférieurs à d1d_1 et d2d_2, leur produit est un polynôme de degré inférieur à d1+d2d_1 + d_2, dont au plus les n=d1+d2+1n = d_1 + d_2 + 1 premiers coefficients sont non nuls.

La transformée de Fourier discrète (de taille nn) de PQPQ s’obtient en évaluant les (PQ)(ωk)=P(ωk)Q(ωk)(PQ)(\omega^k) = P(\omega^k)Q(\omega^k) pour 0k<n0 \leq k \lt n.

En particulier, on remarque que : la transformée de PQPQ est égale à la multiplication point par point des transformées de PP et de QQ.

Notons que toutes les transformées considérées ici sont de taille nn, une séquence de taille m<nm < n doit donc être complétée par nmn - m zéros pour calculer sa transformée.

La transformée de Fourier rapide

On suppose ici que nn est une puissance de 22 : n=2ln = 2^l pour un certain lNl \in \N.

On peut écrire

P=a0+a1X+a2X2+...+an1Xn1=P0(X2)+XP1(X2){P = a_0 + a_1 X + a_2 X^2 + ... + a_{n-1} X^{n-1}} = {P_0(X^2) + X P_1(X^2)}

P0=a0+a2X+...+an2Xn21{P_0 = a_0 + a_2 X + ... + a_{n-2} X^{\frac{n}{2}-1}} et P1=a1+a3X+...+an1Xn21{P_1 = a_1 + a_3 X + ... + a_{n-1} X^{\frac{n}{2} - 1}}.

En particulier, on a pour 0k<n0 \leq k \lt n : P(ωk)=P0(ω2k)+ωkP1(ω2k)P(\omega^k) = {P_0(\omega^{2k}) + \omega^k P_1(\omega^{2k})}

Dans l’expression précédente, ω2k=(ω2)k\omega^{2k} = (\omega^2)^kω2\omega^2 est une racine n2\frac{n}{2}-ème primitive de l’unité.

On a donc : P(ωk)=P0((ω2)kmodn2)+ωkP1((ω2)kmodn2)P(\omega^k) = {P_0((\omega^2)^{k \mod \frac{n}{2}}) + \omega^k P_1((\omega^2)^{k \mod \frac{n}{2}})}

Cette expression fait intervenir l’évaluation des polynômes P0P_0 et P1P_1 sur n2\frac{n}{2} racines n2\frac{n}{2}-èmes de l’unité.

C’est précisément ce que calculent les transformées de taille n2\frac{n}{2} de P0P_0 et P1P_1. On en déduit donc un algorithme de calcul de la transformée qui réduit le problème en deux sous-problèmes de taille n2\frac{n}{2} :

from math import e, pi
def fft(L, sign=-1):
n = len(L)
if n == 1:
return L
w_n = e**(sign * 2j * pi / n)
w = 1
F_even = fft(L[::2], sign)
F_odd = fft(L[1::2], sign)
F = [0] * n
for k in range(n // 2):
F[k] = F_even[k] + w * F_odd[k]
F[k + n // 2] = F_even[k] - w * F_odd[k]
w = w * w_n
return F

On désignera maintenant par FFT la transformée de Fourier discrète d’une séquence obtenue par cette méthode.

Note : l’algorithme implémenté ici est celui décrit dans le chapitre 30 de Introduction to Algorithms, que je vous recommande pour une présentation plus détaillée.

Complexité de l’algorithme de FFT

La fonction fft effectue un nombre proportionnel à nn d’opérations et exécute deux appels récursifs à elle-même sur des séquences de taille n2\frac{n}{2}.

Le nombre d’opération effectuée est donc de l’ordre de n+2n2+22n22+...+2ln2l=n(l+1)=O(nlogn){n + 2 \frac{n}{2} + 2^2 \frac{n}{2^2} + ... + 2^l \frac{n}{2^l}} = n(l+1) = O(n \log n)

Transformation inverse

On a vu plus haut que la matrice d’inversion s’obtenait en considérant la matrice de la transformée, en remplaçant ω\omega par son conjugué et en divisant par nn.

On en déduit la fonction inverse de la transformée ci-dessus :

def ifft(L):
return [c / len(L) for c in fft(L, sign=1)]

Produit de polynômes

En particulier, on en déduit un algorithme qui calcule le produit de deux polynômes en procédant ainsi :

  1. Déterminer le plus petit nn supérieur à degP+degQ+1\deg P + \deg Q + 1 qui soit une puissance de 22 (en O(logn)O(\log n) instructions).
  2. Calculer la FFT de PP et QQ, en ajoutant des zéros de padding pour avoir des séquences de tailles nn (en O(nlogn)O(n \log n) instructions).
  3. Multiplier point par point les deux FFT obtenues (en O(n)O(n) instructions).
  4. Prendre l’inverse de cette FFT : on obtient les coefficients du polynôme PQPQ (en O(nlogn)O(n \log n) instructions).

Cet algorithme a une complexité en O(nlogn)O(n \log n) alors que la multiplication naïve (en appliquant la formule) amènerait à une multiplication en O(n2)O(n^2) opérations.

En voici une implémentation en Python :

def multiply(P, Q):
n = 1
while n < len(P) + len(Q) - 1:
n <<= 1
P = P + [0] * (n - len(P))
Q = Q + [0] * (n - len(Q))
F = [a * b for a, b in zip(fft(P), fft(Q))]
return ifft(F)
>>> multiply([1, 2, 3], [4, 5])
[(4-6.661338147750939e-16j),
(13+6.123233995736765e-17j),
(22+6.661338147750939e-16j),
(15-6.123233995736765e-17j)]

Notons que cet algorithme a le défaut de propager des erreurs d’approximation des flottants comme ci-dessus. On peut nettoyer les coefficients obtenus à l’aide de Numpy :

>>> import numpy as np
>>> P = multiply([1, 2, 3], [4, 5])
>>> np.trim_zeros(np.round(np.real(P)).astype(np.int32), "b")
array([ 4, 13, 22, 15], dtype=int32)

Algorithme de Schönhage-Strassen

Une variante de l’algorithme décrit précédemment est utilisée dans un algorithme de multiplication d’entiers en O(nlognloglogn){O(n \log n \cdot \log \log n)}, l’algorithme de Schönhage-Strassen.

Nous avons vu précédemment que le calcul dans C\Complex de la transformée introduit inévitablement des erreurs d’approximations des flottants. Ceci n’est pas acceptable pour un algorithme de multiplication d’entiers qui doit obtenir un résultat exact.

Transformée sur un anneau quelconque

La transformée que nous avons décrite repose sur des opérations d’additions et multiplications et sur l’existence de racines primitives nn-èmes de l’unité.

Les résultats que nous avons obtenus, notamment la formule d’inversion et l’algorithme de calcul du produit de polynômes restent vrais en utilisant des séquences à valeurs dans d’autres anneaux que C\Complex.

Il est notamment possible de considérer des séquences dans des anneaux finis, où les calculs peuvent être effectués de façon exacte par un ordinateur.

C’est ce type d’anneau qui est utilisé dans l’algorithme de Schönhage-Strassen (voir ici).

Notons notamment l’utilisation de cet algorithme dans V8 (le moteur Javascript utilisé notamment par Google Chrome et Node.js) pour multiplier des entiers de plus de 15001500 chiffres (voir le code source).



Maths et applications, avec les mains et avec du code 💻
Suivre sur Twitter