Précédent Remonter Suivant

Chapitre 10  Polynômes et transformée de Fourier

Nous allons donner quelques idées sur la réalisation de bibliothèques de fonctions s'appliquant à un domaine commun, en l'illustrant sur un exemple, celui des calculs sur les polynômes à coefficients entiers. Une bonne référence pour les algorithmes décrits dans ce chapitre est [Knu81].

Comment écrit-on une bibliothèque ? On commence d'abord par choisir les objets de base, puis on leur adjoint quelques prédicats, des primitives courantes (fabrication, entrées sorties, test d'égalité, etc.). Puis dans un deuxième temps, on construit des fonctions un peu plus complexes, et on poursuit en assemblant des fonctions déjà construites.

10.1  La classe Polynome

Nous décidons de travailler sur des polynômes à coeffients entiers, que nous supposerons ici être de type long1. Un polynôme P(X) = pd Xd + ⋯ + p0 a un degré d, qui est un entier positif ou nul si P n'est pas identiquement nul et -1 sinon (par convention).

10.1.1  Définition de la classe

Cela nous conduit à définir la classe, ainsi que le constructeur associé qui fabrique un polynôme dont tous les coefficients sont nuls :
class Polynome{
    int deg;
    long[] coeff;

    Polynome(int d){
        this.deg = d;
        this.coeff = new long[d+1];
    }
}
Nous faisons ici la convention que les arguments d'appel d'une fonction correspondent à des polynômes dont le degré est exact, et que la fonction retourne un polynôme de degré exact. Autrement dit, si P est un paramètre d'appel d'une fonction, on suppose que P.deg contient le degré de P, c'est-à-dire que P est nul si P.deg ==-1 et P.coeff[P.deg] n'est pas nul sinon.

10.1.2  Création, affichage

Quand on construit de nouveaux objets, il convient d'être capable de les créer et manipuler aisément. Nous avons déjà écrit un constructeur, mais nous pouvons avoir besoin par exemple de copier un polynôme :
  static Polynome copier(Polynome P){
      Polynome Q = new Polynome(P.deg);

      for(int i = 0; i <= P.deg; i++)
          Q.coeff[i] = P.coeff[i];
      return Q;
  }

On écrit maintenant une fonction toString() qui permet d'afficher un polynôme à l'écran. On peut se contenter d'une fonction toute simple :
  public String toString(){
      String s = "";

      for(int i = this.deg; i >= 0; i--){
          s = s.concat("+"+this.coeff[i]+"*X^"+i);
      }
      if(s == "") return "0";
      else return s;
  }

Si on veut tenir compte des simplifications habituelles (pas d'affichage des coefficients nuls de P sauf si P=0, 1 X1 est généralement écrit X), il vaut mieux écrire la fonction de la figure 10.1.

    public String toString(){
        String s = "";
        long coeff;
        boolean premier = true;

        for(int i = this.deg; i >= 0; i--){
            coeff = this.coeff[i];
            if(coeff != 0){
                // on n'affiche que les coefficients non nuls
                if(coeff < 0){
                    s = s.concat("-");
                    coeff = -coeff;
                }
                else
                    // on n'affiche "+" que si ce n'est pas
                    // premier coefficient affiché
                    if(!premier) s = s.concat("+");
                // traitement du cas spécial "1"
                if(coeff == 1){
                    if(i == 0)
                        s = s.concat("1");
                }
                else{
                    s = s.concat(coeff+"");
                    if(i > 0)
                        s = s.concat("*");
                }
                // traitement du cas spécial "X"
                if(i > 1)
                    s = s.concat("X^"+i);
                else if(i == 1)
                    s = s.concat("X");
                // à ce stade, un coefficient non nul
                // a été affiché
                premier = false;
            }
        }
        // le polynôme nul a le droit d'être affiché
        if(s == "") return "0";
        else return s;
    }

Figure 10.1 : Fonction d'affichage d'un polynôme.


10.1.3  Prédicats

Il est commode de définir des prédicats sur les objets. On programme ainsi un test d'égalité à zéro :
    static boolean estNul(Polynome P){
        return P.deg == -1;
    }
De même, on ajoute un test d'égalité :
    static boolean estEgal(Polynome P, Polynome Q){
        if(P.deg != Q.deg) return false;
        for(int i = 0; i <= P.deg; i++)
            if(P.coeff[i] != Q.coeff[i])
                return false;
        return true;
    }

10.1.4  Premiers tests

Il est important de tester le plus tôt possible la bibliothèque en cours de création, à partir d'exemples simples et maîtrisables. On commence par exemple par écrire un programme qui crée le polynôme P(X)=2 X+1, l'affiche à l'écran et teste s'il est nul :
class TestPolynome{
  public static void main(String[] args){
      Polynome P;

      // création de 2*X+1
      P = new Polynome(1);
      P.coeff[1] = 2;
      P.coeff[0] = 1;
      System.out.println("P="+P);
      System.out.println("P == 0 ? " + Polynome.estNul(P));
  }
}
L'exécution de ce programme donne alors :
P=2*X+1
P == 0? false
Nous allons avoir besoin d'entrer souvent des polynômes et il serait souhaitable d'avoir un moyen plus simple que de rentrer tous les coefficients les uns après les autres. On peut décider de créer un polynôme à partir d'une chaîne de caractères formattée avec soin. Un format commode pour définir un polynôme est une chaîne de caractères s de la forme "deg s[deg] s[deg-1] ... s[0]" qui correspondra au polynôme P(X) = sdeg Xdeg + ⋯ + s0. Par exemple, la chaîne "1 1 2" codera le polynôme X+2. La fonction convertissant une chaîne au bon format en polynôme est alors :
  static Polynome deChaine(String s){
      Polynome P;
      long[] tabi = TC.longDeChaine(s);

      P = new Polynome((int)tabi[0]);
      for(int i = 1; i < tabi.length; i++)
          P.coeff[i-1] = tabi[i];
      return P;
  }
(la fonction TC.longDeChaine appartient à la classe TC décrite en annexe) et elle est utilisée dans TestPolynome de la façon suivante :
    P = Polynome.deChaine("1 1 2"); // c'est X+2
Une fois définis les objets de base, il faut maintenant passer aux opérations plus complexes.

10.2  Premières fonctions

10.2.1  Dérivation

La dérivée du polynôme 0 est 0, sinon la dérivée de P(X) = ∑i=0d pi Xi est :
P'(X) =
d
i=1
i pi Xi-1.
On écrit alors la fonction :
  static Polynome deriver(Polynome P){
      Polynome dP;

      if(estNul(P)) return copier(P);
      dP = new Polynome(P.deg - 1);
      for(int i = P.deg; i >= 1; i--)
          dP.coeff[i-1] = i * P.coeff[i];
      return dP;
  }

10.2.2  Évaluation ; schéma de Horner

Passons maintenant à l'évaluation du polynôme P(X) = ∑i=0d pi Xi en la valeur x. La première solution qui vient à l'esprit est d'appliquer la formule en calculant de proche en proche les puissances de x. Cela s'écrit :
  // évaluation de P en x
  static long evaluer(Polynome P, long x){
      long Px, xpi;

      if(estNul(P)) return 0;
      // Px contiendra la valeur de P(x)
      Px = P.coeff[0];
      xpi = 1;
      for(int i = 1; i <= P.deg; i++){
          // calcul de xpi = x^i
          xpi *= x;
          Px += P.coeff[i] * xpi;
      }
      return Px;
  }
Cette fonction fait 2d multiplications et d additions. On peut faire mieux en utilisant le schéma de Horner :
P(x) = (⋯ ((pd x + pd-1) x + pd-2) x+ ⋯ ) x + p0.
La fonction est alors :
  static long Horner(Polynome P, long x){
      long Px;

      if(estNul(P)) return 0;
      Px = P.coeff[P.deg];
      for(int i = P.deg-1; i >= 0; i--){
          // à cet endroit, Px contient:
          // p_d*x^(d-i-1) + ... + p_{i+1}
          Px *= x;
          // Px contient maintenant
          // p_d*x^(d-i) + ... + p_{i+1}*x
          Px += P.coeff[i];
      }
      return Px;
  }
On ne fait plus désormais que d multiplications et d additions. Notons au passage que la stabilité numérique serait meilleure, si x était un nombre flottant.

10.3  Addition, soustraction

Si P(X) = ∑i=0n pi Xi, Q(X) = ∑j=0m qj Xj, alors
P(X)+Q(X) =
min(n, m)
k=0
(pk+qk) Xk +
n
i = min(n, m)+1
pi Xi +
m
j = min(n, m)+1
qj Xj.
Le degré de P+Q sera inférieur ou égal à max(n, m) (attention aux annulations).

Le code pour l'addition est alors :
  static Polynome plus(Polynome P, Polynome Q){
      int maxdeg = (P.deg >= Q.deg ? P.deg : Q.deg);
      int mindeg = (P.deg <= Q.deg ? P.deg : Q.deg);
      Polynome R = new Polynome(maxdeg);

      for(int i = 0; i <= mindeg; i++)
          R.coeff[i] = P.coeff[i] + Q.coeff[i];
      for(int i = mindeg+1; i <= P.deg; i++)
          R.coeff[i] = P.coeff[i];
      for(int i = mindeg+1; i <= Q.deg; i++)
          R.coeff[i] = Q.coeff[i];
      trouverDegre(R);
      return R;
  }

Comme il faut faire attention au degré du résultat, qui est peut-être plus petit que prévu, on a dû introduire une nouvelle primitive qui se charge de mettre à jour le degré de P (on remarquera que si P est nul, le degré sera bien mis à -1) :
  // vérification du degré
  static void trouverDegre(Polynome P){
      while(P.deg >= 0){
          if(P.coeff[P.deg] != 0)
              break;
          else
              P.deg -= 1;
      }
  }

On procède de même pour la soustraction, en recopiant la fonction précédente, les seules modifications portant sur les remplacements de + par - aux endroits appropriés.

Il importe ici de bien tester les fonctions écrites. En particulier, il faut vérifier que la soustraction de deux polynômes identiques donne 0. Le programme de test contient ainsi une soustraction normale, suivie de deux soustractions avec diminution du degré :
static void testerSous(){
    Polynome P, Q, S;

    P = Polynome.deChaine("1 1 1"); // X+1
    Q = Polynome.deChaine("2 2 2 2"); // X^2+X+2
    System.out.println("P="+P+" Q="+Q);
    System.out.println("P-Q="+Polynome.sous(P, Q));
    System.out.println("Q-P="+Polynome.sous(Q, P));

    Q = Polynome.deChaine("1 1 0"); // X
    System.out.println("Q="+Q);
    System.out.println("P-Q="+Polynome.sous(P, Q));
    System.out.println("P-P="+Polynome.sous(P, P));
}
dont l'exécution donne :
P=X+1 Q=2*X^2+2*X+2
P-Q=-2*X^2-X-1
Q-P=2*X^2+X+1
Q=1
P-Q=X
P-P=0

10.4  Deux algorithmes de multiplication

10.4.1  Multiplication naïve

Soit P(X) = ∑i=0n pi Xi, Q(X) = ∑j=0m qj Xj, alors



Le code correspondant en Java est :
static Polynome mult(Polynome P, Polynome Q){
    Polynome R;

    if(estNul(P)) return copier(P);
    else if(estNul(Q)) return copier(Q);
    R = new Polynome(P.deg + Q.deg);
    for(int i = 0; i <= P.deg; i++)
        for(int j = 0; j <= Q.deg; j++)
            R.coeff[i+j] += P.coeff[i] * Q.coeff[j];
    return R;
}

10.4.2  L'algorithme de Karatsuba

Nous allons utiliser une approche diviser pour résoudre de la multiplication de polynômes.

Comment fait-on pour multiplier deux polynômes de degré 1 ? On écrit :
P(X) = p0 + p1 X,    Q(X) = q0 + q1 X,
et on va calculer
R(X) = P(X) Q(X) = r0 + r1 X + r2 X2,
avec
r0 = p0 q0,    r1 = p0 q1 + p1 q0,    r2 = p1 q1.
Pour calculer le produit R(X), on fait 4 multiplications sur les coefficients, que nous appelerons multiplication élémentaire et dont le coût sera l'unité de calcul pour les comparaisons à venir. Nous négligerons les coûts d'addition et de soustraction.

Si maintenant P est de degré n-1 et Q de degré n-1 (ils ont donc n termes), on peut écrire :
P(X) = P0(X) + Xm P1(X),    Q(X) = Q0(X) + Xm Q1(X),
m = ⌈ n/2 ⌉, avec P0 et Q0 de degré m-1 et P1, Q1 de degré n-1-m. On a alors :
R(X) = P(X) Q(X) = R0(X) + Xm R1(X) + X2 m R2(X),
R0 = P0 Q0,    R1 = P0 Q1 + P1 Q0,    R2 = P1 Q1.
Notons M(d) le nombre de multiplications élémentaires nécessaires pour calculer le produit de deux polynômes de degré d-1. On vient de voir que :
M(21) = 4 M(20).
Si n = 2t, on a m = 2t-1 et :
M(2t) = 4 M(2t-1) = O(22t) = O(n2).

L'idée de Karatsuba est de remplacer 4 multiplications élémentaires par 3, en utilisant une approche dite évaluation/interpolation. On sait qu'un polynôme de degré n est complètement caractérisé soit par la donnée de ses n+1 coefficients, soit par ses valeurs en n+1 points distincts (en utilisant par exemple les formules d'interpolation de Lagrange). L'idée de Karatsuba est d'évaluer le produit P Q en trois points 0, 1 et ∞. On écrit :
R0 = P0 Q0,   R2 = P1 Q1,   R1 = (P0+P1) (Q0+Q1) - R0 - R1
ce qui permet de ramener le calcul des Ri à une multiplication de deux polynômes de degré m-1, et deux multiplications en degré n-1-m plus 2 additions et 2 soustractions. Dans le cas où n = 2t, on obtient :
K(2t) = 3 K(2t-1) = O(3t) = O(nlog2 3) = O(n1.585).
La fonction K vérifie plus généralement l'équation fonctionnelle :



et son comportement est délicat à prédire (on montre qu'elle a un comportement fractal).

Première implantation

Nous allons implanter les opérations nécessaires aux calculs précédents. On a besoin d'une fonction qui récupère P0 et P1 à partir de P. On écrit donc une fonction :
// crée le polynôme 
//     P[début]+P[début+1]*X+...+P[fin]*X^(fin-début)
static Polynome extraire(Polynome P, int début, int fin){
    Polynome E = new Polynome(fin-début);

    for(int i = début; i <= fin; i++)
        E.coeff[i-début] = P.coeff[i];
    trouverDegre(E);
    return E;
}
Quel va être le prototype de la fonction de calcul, ainsi que les hypothèses faites en entrée ? Nous décidons ici d'utiliser :
// ENTRÉE: deg(P) = deg(Q) <= n-1, 
//     P.coeff et Q.coeff sont de taille >= n;
// SORTIE: R tq R = P*Q et deg(R) <= 2*(n-1).
static Polynome Karatsuba(Polynome P, Polynome Q, int n){
Nous fixons donc arbitrairement le degré de P et Q à n-1. Une autre fonction est supposée être en charge de la normalisation des opérations, par exemple en créant des objets de la bonne taille.

On remarque également, avec les notations précédentes, que P0 et Q0 sont de degré m-1, qui est toujours plus grand que le degré de P1 et Q1, à savoir n-m-1. Il faudra donc faire attention au calcul de la somme P0+P1 (resp. Q0+Q1) ainsi qu'au calcul de R1.

La fonction complète est donnée dans la table 10.2.

static Polynome Karatsuba(Polynome P, Polynome Q, int n){
    Polynome P0, P1, Q0, Q1, SP, SQ, R0, R1, R2, R;
    int m;

    if(n <= 1)               // (cf. remarque 1)
        return mult(P, Q);
    m = n/2;
    if((n % 2) == 1) m++;
    // on multiplie P = P0 + X^m * P1 avec Q = Q0 + X^m * Q1
    // deg(P0), deg(Q0) <= m-1
    // deg(P1), deg(Q1) <= n-1-m <= m-1
    P0 = extraire(P, 0, m-1);
    P1 = extraire(P, m, n-1);
    Q0 = extraire(Q, 0, m-1);
    Q1 = extraire(Q, m, n-1);

    // R0 = P0*Q0 de degré 2*(m-1)
    R0 = Karatsuba(P0, Q0, m);

    // R2 = P2*Q2 de degré 2*(n-1-m)
    R2 = Karatsuba(P1, Q1, n-m);

    // R1 = (P0+P1)*(Q0+Q1)-R0-R2
    // deg(P0+P1), deg(Q0+Q1) <= max(m-1, n-1-m) = m-1
    SP = plusKara(P0, P1, m-1);       // (cf. remarque 2)
    SQ = plusKara(Q0, Q1, m-1);
    R1 = Karatsuba(SP, SQ, m);
    R1 = sous(R1, R0);
    R1 = sous(R1, R2);
    // on reconstruit le résultat
    // R = R0 + X^m * R1 + X^(2*m) * R2
    R = new Polynome(2*(n-1));
    for(int i = 0; i <= R0.deg; i++)
        R.coeff[i] = R0.coeff[i];
    for(int i = 0; i <= R2.deg; i++)
        R.coeff[2*m + i] = R2.coeff[i];
    for(int i = 0; i <= R1.deg; i++)
        R.coeff[m + i] += R1.coeff[i];
    trouverDegre(R);
    return R;
  }

Figure 10.2 : Algorithme de Karatsuba.


Expliquons la remarque 1. On décide pour l'instant d'arrêter la récursion quand on doit multiplier deux polynômes de degré 0 (donc n=1).

La remarque 2 est justifiée par notre invariant de fonction : les degrés de SP et SQ (ou plus exactement la taille de leurs tableaux de coefficients), qui vont être passés à Karatsuba doivent être m-1. Il nous faut donc modifier l'appel plus(P0, P1); en celui plusKara(P0, P1, m-1); qui retourne la somme de P0 et P1 dans un polynôme dont le nombre de coefficients est toujours m, quel que soit le dégré de la somme (penser que l'on peut tout à fait avoir P0 = 0 et P1 de degré m-2).
// ENTREE: deg(P), deg(Q) <= d.
// SORTIE: P+Q dans un polynôme R tel que R.coeff a taille
//         d+1.
static Polynome plusKara(Polynome P, Polynome Q, int d){
    int mindeg = (P.deg <= Q.deg ? P.deg : Q.deg);
    Polynome R = new Polynome(d);

    //PrintK("plusKara("+d+", "+mindeg+"): "+P+" "+Q);
    for(int i = 0; i <= mindeg; i++)
        R.coeff[i] = P.coeff[i] + Q.coeff[i];
    for(int i = mindeg+1; i <= P.deg; i++)
        R.coeff[i] = P.coeff[i];
    for(int i = mindeg+1; i <= Q.deg; i++)
        R.coeff[i] = Q.coeff[i];
    return R;
  }

Comment teste-t-on un tel programme ? Tout d'abord, nous avons de la chance, car nous pouvons comparer Karatsuba à mul. Un programme test prend en entrée des couples de polynômes (P, Q) de degré n et va comparer les résultats des deux fonctions. Pour ne pas avoir à rentrer des polynômes à la main, on construit une fonction qui fabrique des polynômes (unitaires) ``aléatoires'' à l'aide d'un générateur créé pour la classe :
static Random rd = new Random();

static Polynome aleatoire(int deg){
    Polynome P = new Polynome(deg);

    P.coeff[deg] = 1;
    for(int i = 0; i < deg; i++)
        P.coeff[i] = rd.nextLong();
    return P;
}
La méthode rd.nextLong() retourne un entier ``aléatoire'' de type long fabriqué par le générateur rd.

Le programme test, dans lequel nous avons également rajouté une mesure du temps de calcul est alors :
// testons Karatsuba sur n polynômes de degré deg
static void testerKaratsuba(int deg, int n){
    Polynome P, Q, N, K;
    long tN, tK, totN = 0, totK = 0;

    for(int i = 0; i < n; i++){
        P = Polynome.aleatoire(deg);
        Q = Polynome.aleatoire(deg);
        TC.demarrerChrono();
        N = Polynome.mult(P, Q);
        tN = TC.tempsChrono();

        TC.demarrerChrono();
        K = Polynome.Karatsuba(P, Q, deg+1);
        tK = TC.tempsChrono();

        if(! Polynome.estEgal(K, N)){
            System.out.println("Erreur");
            System.out.println("P*Q(norm)=" + N);
            System.out.println("P*Q(Kara)=" + K);
            for(int i = 0; i <= N.deg; i++){
                if(K.coeff[i] != N.coeff[i])
                    System.out.print(" "+i);
            }
            System.out.println("");
            System.exit(-1);
        }
        else{
            totN += tN;
            totK += tK;
        }
    }
    System.out.println(deg+" N/K = "+totN+" "+totK);
}
Que se passe-t-il en pratique ? Voici des temps obtenus avec le programme précédent, pour 100 ≤ deg ≤ 1000 par pas de 100, avec 10 couples de polynômes à chaque fois :
Test de Karatsuba
100 N/K = 2 48
200 N/K = 6 244
300 N/K = 14 618
400 N/K = 24 969
500 N/K = 37 1028
600 N/K = 54 2061
700 N/K = 74 2261
800 N/K = 96 2762
900 N/K = 240 2986
1000 N/K = 152 3229
Cela semble frustrant, Karatsuba ne battant jamais (et de très loin) l'algorithme naïf sur la plage considérée. On constate cependant que la croissance des deux fonctions est à peu près la bonne, en comparant par exemple le temps pris pour d et 2 d (le temps pour le calcul naïf est multiplié par 4, le temps pour Karatsuba par 3).

Comment faire mieux ? L'astuce classique ici est de décider de repasser à l'algorithme de multiplication classique quand le degré est petit. Par exemple ici, on remplace la ligne repérée par la remarque 1 en :
if(n <= 16)
ce qui donne :
Test de Karatsuba
100 N/K = 1 4
200 N/K = 6 6
300 N/K = 14 13
400 N/K = 24 17
500 N/K = 38 23
600 N/K = 164 40
700 N/K = 74 69
800 N/K = 207 48
900 N/K = 233 76
1000 N/K = 262 64
Le réglage de cette constante est critique et dépend de la machine sur laquelle on opère.

Remarques sur une implantation optimale

La fonction que nous avons implantée ci-dessus est gourmande en mémoire, car elle alloue sans cesse des polynômes auxiliaires. Diminuer ce nombre d'allocations (il y en a O(n1.585) également...) est une tâche majeure permettant de diminuer le temps de calcul. Une façon de faire est de travailler sur des polynômes définis par des extraits compris entre des indices de début et de fin. Par exemple, le prototype de la fonction pourrait devenir :
static Polynome Karatsuba(Polynome P, int dP, int fP,
                          Polynome Q, int dQ, int fQ, int n){
qui permettrait de calculer le produit de P' = PfPXfP-dP+⋯ +PdP et Q' = PfQXfQ-dQ+⋯ +QdQ. Cela nous permettrait d'appeler directement la fonction sur P0' et P1' (resp. Q0' et Q1') et éviterait d'avoir à extraire les coefficients.

Dans le même ordre d'idée, l'addition et la soustraction pourraient être faites en place, c'est-à-dire qu'on implanterait plutôt P := P-Q.

10.5  Multiplication à l'aide de la transformée de Fourier*

Quel est le temps minimal requis pour faire le produit de deux polynômes de degré n ? On vient de voir qu'il existe une méthode en O(n1.585). Peut-on faire mieux ? L'approche de Karatsuba consiste à couper les arguments en deux. On peut imaginer de couper en 3, voire plus. On peut démontrer qu'asymptotiquement, cela conduit à une méthode dont le nombre de multiplications élémentaires est O(n1+ε) avec ε>0 aussi petit qu'on le souhaite.

Il existe encore une autre manière de voir les choses. L'algorithme de Karatsuba est le prototype des méthodes de multiplication par évaluation/interpolation. On a calculé R(0), R(1) et R(+∞) et de ces valeurs, on a pu déduire la valeur de R(X). L'approche de Cooley et Tukey consiste à interpoler le produit R sur des racines de l'unité bien choisies.

10.5.1  Transformée de Fourier


Définition 1   Soit ω∈C et N un entier. La transformée de Fourier est une application
Fω : CN CN
  (a0, a1, …, aN-1) (a0*, a1*, …, aN-1*)
a i* =
N-1
j=0
ωi j aj
pour 0 ≤ iN-1.

Proposition 9   Si ω est une racine primitive N-ième de l'unité, (i.e., ωN=1 et ωi ≢ 1 pour 1 ≤ i < N), alors Fω est une bijection et
Fω-1 =
1
N
Fω-1.

Démonstration : Posons
αi =
1
N
N-1
k=0
ω-i k ak*.
On calcule
N αi =
 
k
ω-i k
 
j
ωk j aj =
 
j
aj
 
k
ωk (j-i) =
 
j
aj Si, j.
Si i = j, on a Si, j = N et si ji, on a
Si, j =
N-1
k=0
j-i)k =
1-(ωj-i)N
1-ωj-i
= 0.

10.5.2  Application à la multiplication de polynômes

Soient
P(X) =
n-1
i=0
pi Xi,    Q(X) =
n-1
i=0
qi Xi
deux polynômes dont nous voulons calculer le produit :
R(X) =
2n-1
i=0
ri Xi
(avec r2n-1 = 0). On utilise une transformée de Fourier de taille N = 2n avec les vecteurs :
p = (p0, p1, …, pn-1, 0, 0, …, 0n  termes),
q = (q0, q1, …, qn-1, 0, 0, …, 0n  termes).
Soit ω une racine primitive 2n-ième de l'unité. La transformée de p
Fω(p) = (p0*, p1*, …, p2n-1*)
n'est autre que :
(P0), P1), …, P2n-1)).
De même pour q, de sorte que le produit terme à terme des deux vecteurs :
Fω(p) ⊗ Fω(q) = (p0* q0*, p1* q1*, …, p2n-1* q2n-1*)
donne en fait les valeurs de R(X) = P(X) Q(X) en les racines de l'unité, c'est-à-dire Fω(R) ! Par suite, on retrouve les coefficients de P en appliquant la transformée inverse.

Un algorithme en pseudo-code pour calculer R est alors :

10.5.3  Transformée rapide

Transformée multiplicative

Si l'on s'y prend naïvement, le calcul des xi* définis par
x k* =
N-1
m=0
xm ωm k, 0 ≤ kN-1
prend N2 multiplications2.

Supposons que l'on puisse écrire N sous la forme d'un produit de deux entiers plus grands que 1, soit N = N1 N2. On peut écrire :
m = N1 m2 + m1,    k = N2 k1 + k2
avec 0 ≤ m1, k1 < N1 et 0 ≤ m2, k2 < N2. Cela conduit à récrire :
xk* =
N1-1
m1=0
ωN2 m1 k1
ωm1 k2
N2-1
m2=0
xN1 m2+m1 ωN1 m2 k2.

On peut montrer sans grande difficulté que ω1 = ωN2 est racine primitive N1-ième de l'unité, ω2 = ωN1 est racine primitive N2-ième. On se ramène alors à calculer :
xk* =
N1-1
m1=0
ω1m1 k1
ωm1 k2
N2-1
m2=0
xN1 m2+m1 ω2m2 k2.
La deuxième somme est une transformée de longueur N2 appliquée aux nombres
(xN1 m2+m1)0≤ m2<N2.
Le calcul se fait donc comme celui de N1 transformées de longueur N2, suivi de multiplications par des facteurs ωm1 k2, suivies elles-mêmes de N2 transformées de longueur N1.

Le nombre de multiplications élémentaires est alors :
N1 (N22) + N1 N2 + N2 (N12) = N1 N2 (N1+N2+1)
ce qui est en général plus petit que (N1 N2)2.

Le cas magique N = 2t

Appliquons le résultat précédent au cas où N1=2 et N2 = 2t-1. Les calculs que nous devons effectuer sont :
xk* =
N/2-1
m=0
x2m2)m k + ωk
N/2-1
m=0
x2m+12)m k
xk+N/2* =
N/2-1
m=0
x2m2)m k - ωk
N/2-1
m=0
x2m+12)m k
car ωN/2=-1. Autrement dit, le calcul se divise en deux morceaux, le calcul des moitiés droite et gauche du signe + et les résultats sont réutilisés dans la ligne suivante.

Pour insister sur la méthode, nous donnons ici le pseudo-code en Java sur des vecteurs de nombres réels :
static double[] FFT(double[] x, int N, double omega){
    double[] X = new double[N], xx, Y0, Y1;
    double omega2, omegak;

    if(N == 2){
        X[0] = x[0] + x[1];
        X[1] = x[0] - x[1];
        return X;
    }
    else{
        xx = new double[N/2];
        omega2 = omega*omega;
        for(m = 0; m < N/2; m++) xx[m] = x[2*m];
        Y0 = FFT(xx, N/2, omega2);
        for(m = 0; m < N/2; m++) xx[m] = x[2*m+1];
        Y1 = FFT(xx, N/2, omega2);
        omegak = 1.; // pour omega^k
        for(k = 0; k < N/2; k++){
            X[k]     = Y0[k] + omegak*Y1[k];
            X[k+N/2] = Y0[k] - omegak*Y1[k];
            omegak = omega * omegak;
        }
        return X;
    }
}

Le coût de l'algorithme est alors F(N) = 2 F(N/2) + N/2 multiplications élémentaires. On résout la récurrence à l'aide de l'astuce suivante :
F(N)
N
=
F(N/2)
N/2
+ 1/2 = F(1) + t/2 = t/2
d'où F(N) = 1/2 N log2 N. Cette variante a ainsi reçu le nom de transformée de Fourier rapide (Fast Fourier Transform ou FFT).

À titre d'exemple, si N=210, on fait 5× 210 multiplications au lieu de 220.

Remarques complémentaires

Nous avons donné ici une brève présentation de l'idée de la FFT. C'est une idée très importante à utiliser dans tous les algorithmes basés sur les convolutions, comme par exemple le traitement d'images, le traitement du signal, etc.

Il y a des milliards d'astuces d'implantation, qui s'appliquent par exemple aux problèmes de précision. C'est une opération tellement critique dans certains cas que du hardware spécifique existe pour traiter des FFT de taille fixe. On peut également chercher à trouver le meilleur découpage possible quand N n'est pas une puissance de 2. Le lecteur intéressé est renvoyé au livre de Nussbaumer [Nus82].

Signalons pour finir que le même type d'algorithme (Karatsuba, FFT) est utilisé dans les calculs sur les grands entiers, comme cela est fait par exemple dans la bibliothèque multiprécision GMP.


1
stricto sensu, nous travaillons en fait dans l'anneau des polynômes à coefficients définis modulo 264.
2
On remarque que ωmk = ω(m k)modN et le précalcul des ωi pour 0 ≤ i < N coûte N multiplications élémentaires.

Précédent Remonter Suivant