Cours C++ : les opérateurs, Variables et constantes

Cours C++ : les opérateurs, Variables et constantes

...

1.2          Analyse asymptotique des algorithmes

Le but de l’algorithmique en général, et de l’algorithmique numérique en par- ticulier, est de concevoir l’algorithme le plus efficace possible pour résoudre un problème particulier posé. Pour cela, il est nécessaire de mesurer son efficacité afin de pouvoir les comparer entre eux. Cette section introduit quelques notions de base relatives à l’analyse des algorithmes.

1.2.1      Quelques concepts théoriques de l’algorithmique

Le choix de l’algorithme le plus efficace devient crucial lorsque la taille des données à traiter, notée n, devient grande. Aussi, nous nous intéressons à la croissance asymptotique du temps de calcul et de la place mémoire lorsque n tend vers l’infini.

Les fonctions usuelles 1, log n, n, n log n, n2, n3,. . . , 2n, dont les ordres de grandeur forment une suite croissante, peuvent constituer une échelle de com- paraison pour les algorithmes. Si, pour réaliser un même calcul, un premier algorithme a un temps de calcul T1(n) = c1n2 et un second algorithme un temps T2(n) = c2n, où c1 et c2 sont deux constantes strictement positives, alors pour n assez grand on aura T1(n) “ T2(n). Dans ce cas, nous dirons que le second algorithme est asymptotiquement plus rapide que le premier. Ceci reste vrai même si c1 est beaucoup plus grand que c2.

Ceci reste encore vrai, si le second algorithme est exécuté sur une machine lente, avec un langage de programmation peu performant et a été codé par un programmeur débutant, il sera quand même, pour des données assez grandes, plus rapide que le premier algorithme, même si il est exécuté sur un super-calculateur, écrit en langage assembleur très optimisé et codé par un program- meur expert. L’analyse asymptotique des performances permet donc de s’af- franchir de considérations particulières en comparant les algorithmes sur une machine abstraite.

Les performances d’un algorithme peuvent se mesurer en terme de temps de calcul et de place mémoire nécessaires à l’exécution de l’algorithme : plutôt que de considérer un programme particulier qui code cet algorithme sur une machine donnée, par un programmeur expérimenté ou non, avec un langage de programmation performant ou non, imaginons plutôt une machine abstraite idéale et définissons-lui son unité de temps et de mémoire.

Le modèle de calculateur sous-jacent à toutes les analyses de l’efficacité des algorithmes présenté dans ce livre est connu sous le nom de calculateur réel à accès aléatoire (voir [8] par exemple). Pour simplifier, nous supposons que la machine est capable de travailler en précision infinie : autrement dit, nous ne considérons pas la question de la gestion des arrondis du point de vue de la place mémoire ou des temps de calculs. Ce modèle suppose que chaque unité de mémoire peut contenir la représentation d’un nombre réel et que l’accès en lecture ou en écriture à une unité de mémoire est réalisé en un temps con- stant, c’est-à-dire indépendament de l’unité de mémoire accédée. L’ensemble des opérations élémentaires, qui seront exécutées en une unité de temps par notre machine abstraite comprend les quatre opérations algébriques +, , , / ainsi que la comparaison entre deux nombres, le calcul de toutes les fonctions mathématiques usuelles comme log, exp, sin, cos, et partie entière.

Afin de comparer les ordres de grandeur de différentes fonctions asymptotiques, nous définissons les notations suivantes. Soient f et g deux fonctions positives de la variable entière n et à valeurs réelles.

Nous dirons que f est un grand o de g, et nous noterons f (n) = (g(n)), si et seulement si il existe une constante C1 > 0 et un entier n0 tels que

f (n) ™ C1 g(n), ∀n “ n0

Nous dirons que f est un grand oméga de g, et nous noterons f (n) = Ω(g(n)), si et seulement si il existe une constante C2 > 0 et un entier n0 tels que

C2 g(n) ™ f (n), ∀n “ n0

Nous dirons que f est un grand thêta de g, et nous noterons f (n) = Θ(g(n)), si et seulement si il existe deux constantes C1 > 0 et C2 > 0 et un entier n0 tels que

C2 g(n) ™ f (n) ™ C1 g(n), ∀n “ n0

1.2.2      Exemple : le calcul de xn

L’idée la plus simple pour calculer xn est de dérouler une boucle et d’effectuer

n multiplications :

 pow_linear.cc

Pour exécuter ce programme, rien de plus simple :

make pow_linear

./pow_linear 2 8

./pow_linear 3.1 500

Le coût en temps de calcul est évidement Θ(n).

Pouvons-nous faire mieux ? Si n est pair, nous avons la possibilité de décom- poser le calcul en xn/2xn/2 : au lieu de n multiplications, nous n’en effectuons à présent plus que n/2 + 1 pour obtenir le même résultat. De même, si n/2 est pair, nous pouvons re-diviser le produit de xn/2 en xn/4xn/4, si bien que, ap- pliqué récursivement, nous effectuons beaucoup moins de multiplications. Si n est impair, il est possible de décomposer le produit en x(n−1)/2x(n−1)/2x. Nous obtenons une réduction importante du nombre d’opérations en appliquant récursivement cette méthode. L’implémentation correspondante en langague C++ est donnée par :

pow_recursive.cc

Nous pouvons tester cette nouvelle version :

make pow_linear

 ./pow_recursive 2 8

./pow_recursive 3.1 500

Déterminons maintenant le coût en temps de calcul T (n) de cet algorithme. Il se déduit de la formule de récurrence :

T (n) = T (|n/2∫) + Θ(1)

La notation ξ représente sa partie entière du nombre réel ξ. La relation précé- dente exprime que le nombre nécéssaire d’opérations pour calculer xn est celui pour calculer nn/2 plus un nombre constant d’opérations, indépendant de n. Ces dernières opérations incluent tout ce qui ne fait pas partie de l’appel récur- sif. Par récurrence, nous pouvons montrer sans difficulté que

T (n) = Θ(log n)

Nous avons par conséquent amélioré de façon très importante l’algorithme de

xn.

Cette approche de réduction du temps de calcul peut s’étendre à un très grand nombre de situations.

1.2.3      Résolution des récurrences

Plutôt que de résoudre, pour chaque algorithme, des récurrences similaires afin d’obtenir son comportement asymptotique, nous avons intérêt à le résoudre de manière générale. Supposons que nous disposons d’un algorithme dont le temps de calcul vérifie la relation de récurrence :

T (n) = aT (|n/b∫) + f (n)

où a “ 1 et b > 1 sont deux constantes, et f est une fonction réelle de l’entier

1. Ce type de relation est très courant dans les algorithmes récursifs et nous en verrons un exemple à la section suivante pour la multiplication de deux matrices. Cormen, Leiserson, Rivest et Stein [15, p. 86] proposent la résolution cette relation de récurrence en considérant différents cas possibles :

1)            Si  f (n)   =            (nlogb (a)+λ) pour une certaine constante λ > 0, alors

T (n) = Θ(nlogb (a)).

2)            Si que f (n) = Θ(nlogb (a)) alors T (n) = Θ(nlogb (a)).

3)            Si que f (n) = Ω(nlogb (a)+λ) pour une certaine constante λ > 0, et qu’il existe une constante C, 0 < C < 1 et un entier n0 tels que

af (n/b) ™ Cf (n),  ∀n “ n0

alors T (n) = Θ(f (n)).

Dans ce qui précède, nous avons noté logb la fonction logarithme en base b, qui est défini pour tout réel positif ξ par logb(ξ) = log(ξ)/ log(b).

 1.2.4     Algorithme de Strassen pour la multiplication ma- tricielle

La multiplication matricielle C = AB de deux matrices A = (ai,j)1™i,j™n et

B = (bi,j)1™i,j™n est définie par ses coefficients (ci,j)1™i,j™n :

n

ci,j =      ai,kbk,j

k=1

Il est clair que le calcul de chaque coefficient ci,j nécessite exactement n mul- tiplications et n 1 additions. Comme C possède n n coefficients, le nombre total de multiplications pour calculer C = AB est donc Θ(n3). Vous pouvez penser que cette façon de faire ne saurait être améliorée, c’est à dire effectuée asymptotiquement en moins d’itérations : vous avez tord. L’algorithme étudié ici, proposé par Strassen [54] a suscité un grand émoi lors de sa publication en 1969.

Supposons que n est pair et partitionnons les matrices en blocs :

A =         a11         a12

a21         a22

,   B =     b11        b12

b21        b22

,   C =     c11         c12

c21         c22

où les blocs sont de taille n/2. Lorsque n est impair, nous pouvons également partitionner en deux blocs de taille respective (n 1)/2 et (n + 1)/2. Aussi, pour simplifier l’exposé et sans perte de généralité, nous allons supposer que n est une puissance de 2. Le produit C = AB se développe alors en :

c11         =             a11b11 + a12b21

c12         =             a11b12 + a12b22

c21         =             a21b11 + a22b21

c22         =             a21b12 + a22b22

Ce calcul nécessite huit multiplications de matrices de taille n/2. L’idée nullement évidente de Strassen, est que ce calcul peut être effectué en sept multiplications au lieu de huit. Pour cela, calculons tout d’abord les sommes :

a12 − a22

Ensuite viennent

b11 + b12

 Nous pouvons vérifier, par un simple développement, que

c11         =             p3 + p4 − p2 + p6

c12         =             p1 + p2

c21         =             p3 + p4

c22         =             p5 + p1 − p3 − p7

Le temps de calcul T (n) de cet algorithme vérifie la relation de récurrence

T (n) = 7 T (n/2) + c n2

où c > 0 est une constante, indépendante de n. Ce dernier terme inclut toutes les opérations non-récursives telles que les sommes et soustractions de deux matrices. Sans difficulté, nous pouvons montrer par récurrence sur T (n) que

T (n) = Θ .nlog2 7Σ

Ce résultat peut également être obtenu directement en appliquant les formules du paragraphe 1.2.3. Autrement dit, il est possible de calculer le produit de deux matrices carrées de taille     en         log2 7     opérations, avec   ,

ce qui améliore l’algorithme en Θ .n3Σ opérations précédent.

1.2.5      Application à l’inversion d’une matrice

Supposons qu’une matrice C soit l’inverse d’une matrice A. La méthode la plus naturelle pour calculer C à partir de A est de résoudre n systèmes linéaires :

Axi = bi, 1 ™ i ™ n

avec pour second membre bi un vecteur ne contenant que des zéros sauf son i-ème élément qui vaut un. Alors le vecteur xi représente la i-ième colonne de la matrice C. Pour résoudre les n systèmes linéaires, on a recours à une factorisa- tion A = LU où L est une matrice triangulaire inférieure et U est triangulaire supérieure unitaire. Dans le cas où A est symétrique, les variantes sont la fac- torisation de Cholesky A = LLT , où L est triangulaire inférieure, ou bien la factorisation A = LDLT , avec D diagonale et L triangulaire inférieure unitaire (voir par exemple [35], page 266). Rappellons que cette dernière variante, très classique, s’écrit en pseudo-code :

algorithme : factorisation A = LDLT

entrée :  A matrice n      n symétrique définie positive

sortie : L triangulaire inférieure unitaire et D diagonale

début

pour i = 1 : n L(i, i) = 1

pour j = 1 : i − 1 L(i, j) = A(i, j) pour k = 1 : j − 1

L(i, j) = L(i, j) − L(i, k) ∗ D(k, k) ∗ L(j, k) L(i, j) = L(i, j)/D(j, j)

D(i, i) = A(i, i)

pour j = 1 : i − 1

D(i, i) = D(i, i) − L(i, j)2 ∗ D(j, j)

ftn

Le pseudo-code précédent adopte les conventions de notations et de syntaxe du logiciel octave [20] qui implémente le langage matlab, et nous avons introduit l’indentation pour marquer la limite des blocs d’instructions afin d’alléger les notations. Il y a trois boucles imbriquées, et le coût de cette factorisation est clairement Θ(n3). Il est possible d’améliorer légèrement cet algorithme, en util- isant deux fois moins de multiplications et en effectuant l’algorithme sur place, c’est-à-dire en écrivant directement le résultat L et D dans le tableau servant à l’entrée A (voir par exemple [35], page 268). Le temps de calcul restera cepen- dant Θ(n3). La résolution de Ax = b équivaut à LLT x = b et s’effectue en deux étapes : tout d’abord Ly = b puis LT x = y. Chacune de ces deux étapes fait intervenir une matrice triangulaire. Ainsi, une fois L construite, par une descente, puis une remontée, la résolution d’un système linéaire peut s’effectuer en Θ(n2) opérations. Ayant n systèmes linéaires à résoudre, nous obtenons par cette méthode un coût total de la construction de l’inverse C = A−1 en Θ(n3)

opérations.

Revenons à présent à l’approche récursive adoptée pour le produit de deux matrices. En partitionnant A et C, et avec les notations de la section précédente, nous avons, par une généralisation de l’inverse d’une matrice 2 × 2

…

 Nous voyons apparaître un grand nombre de facteurs répétés, que nous pouvons ne calculer qu’une seule fois :

r1            =   a−1   c12         =    r3 r6

r2            =   a21 r1              c21         =    r6 r2

r3            =   r1 a12              r7            = r3 c21

r4            =   a21 r3              c11         = r1 − r7 r5          =    r4 − a22         c22         = −r6

r6            =    r5−1

Ainsi, pour calculer l’inverse d’une matrice de taille n, nous n’avons besoin que de deux inverses de matrices de taille n/2 et de six produits de matrices de taille n/2. Le coût d’une opération d’inversion est donné par :

T (n) = 2 T (n/2) + c nlog2 7

Par récurrence, nous montrons que :

T (n) = Θ .nlog2 7Σ

Autrement dit, le coût de l’inversion d’une matrice est dominé par le coût du produit entre deux matrices.

Le calcul du déterminant suit également cette réduction. De la formule

det(A) = det (a22) det .a11 − a12a−1a21Σ

nous déduisons que le calcul du déterminant d’une matrice d’ordre n se réduit à deux calculs de déterminants d’ordre n/2 plus une inversion et deux produits de matrices d’ordre n/2. Autrement dit, nous retrouvons un coût en Θ nlog2 7 .

1.2.6      La classe valarray(T)

Afin de manipuler des données de grande taille, nous allons avoir recours, tout au long de ce livre, à la classe valarray T . Cette classe, de type tableau,       est disponible dans la librairie standard du C++ et représente un vecteur de données numériques. Elle est définie sous l’espace de nom std dans l’entête (valarray) :

<valarray>

Cette classe est complétée par les opérations d’algèbre  linéaire  usuelles  pour les vecteurs de Rn.

Nous utiliserons souvent la fonction membre shift qui réalise un décalage logique (voir [55, p. 740]). L’instruction v=u.shift(d) renvoie dans v un objet de la classe valarray T de même taille que u, et dont le i-ème élément v[i] est u[i+d] si i+d est dans l’intervalle 0..u.size()-1, et zéro autrement. Ainsi, une valeur positive de d décale les éléments vers la gauche de d places, avec un remplissage par des zéro. Par exemple :

u = (u0, u1, . . . un−2, un−1)

u.shift(1) = (u1, u2, . . . un−1, 0)

u.shift(−1) = (0, u0, . . . un−3, un−2)

En s’inspirant de la syntaxe du logiciel octave [20] qui implémente le langage matlab, nous introduisons, pour notre confort, la fonction range qui opère sur la classe valarray. Voici la correspondance entre les deux syntaxes :

En  clair,  u  [range(first,last)] considère  l’intervalle  [first,last[ fermé à gauche et ouvert à droite, si bien que le nombre d’éléments est last-first. La fonction range est définie par

range.h

Lorsque le pas ne vaut pas un, nous avons

 range.h (suite) 

La classe slice, renvoyée ici par la fonction range, est  définie  avec  la classe valarray dans le fichier d’entête standard <valarray>.

Enfin, nous complétons la classe valarray par deux fonctions qui nous seront utiles tout au long de ce livre et qui ne font pas partie de la librairie standard.

Les fonctions dot(.,.) et norm(.) se définissent simplement par

valarray_util.h

1.2.7      Exercices

Exercice 3. (Factorisation LDLT ) La matrice

0   1

1   2

possède-t’elle une factorisation LDLT ?

Exercice 4. (Factorisation LDLT d’une matrice bande)

1)            Soit b un entier dans l’intervalle [0, n 2]. Montrer que si une matrice symétrique A est (2b + 1)-diagonale, c’est-à-dire que si Ai,j = 0 pour  i  j “ b + 1, alors la factorisation A = LDLT est telle que L est (b + 1)-diagonale, c’est-à-dire Li,j = 0 pour i j “ b + 1. L’entier b est appelé largeur de bande de la matrice A et on dit que la factorisation de Cholesky respecte la structure bande de A.

2)            Écrire une variante de l’algorithme de factorisation LDLT qui exploite le fait que les éléments de A hors de la bande sont nuls. Montrer que le temps de calcul de cet algorithme est

T (n) = Θ .b2nΣ

et que le temps de résolution des systèmes triangulaires, une fois la factorisation formée, est Θ(bn).

Exercice 5. (Algorithme de Strassen)

Déroulez à la main l’algorithme de Strassen pour le produit matriciel

. 1   3 Σ . 6   8 Σ

Exercice 6. (Multiplication de polynômes)

Montrez  comment  multiplier  deux  polynômes  réels  p1(x) = ax + b  et   p2(x) = cx + d en trois multiplications entre réels au lieu de quatre. Indica- tion : une des multiplications est (a + b)(c + d).

Exercice 7. (Multiplication de nombres complexes)

Montrez comment multiplier deux nombres complexes z1 = a + ib et z2 = c + id

en trois multiplications entre réels au lieu de quatre.

Exercice 8. (Matrices denses)

Nous nous proposons d’étudier la programmation la plus simple possible de matrices denses rectangulaires. Ces matrices ne font pas partie de la librairie standard actuelle. Cependant, différents projets récents de standardisation (tnt [40],blitz++ [59],eigen [29],flens [36],armadillo [48],boost/ublas [61]) proposent un interface du type :

dmatrix.h

Le tableau de type valarray dans la zone de données protégées contient les coefficients de la matrice, rangés colonne par colonne.

1)            Écrire le code du constructeur ainsi que les quatre fonctions d’accès.

2)            Écrire le code du produit matrice-vecteur suivant :

 valarray < T >

o p e r a t o r * ( const dmatrix < T >& a , const valarray < T >& x );

Quel est le coût en nombre de multiplications de cette opération ?

3)            Écrire le code du produit matrice-matrice suivant :

dmatrix < T >

o p e r a t o r * ( const dmatrix < T >& a , const dmatrix < T >& b );

L’algorithme sera le plus simple possible, sur la base de trois boucles imbriquées. Montrez que le coût en nombre de multiplication de cette opération est Θ(n3) pour une matrice carrée n × n.

1.2.8      Notes

Les ouvrage actuels de calcul scientifique n’abordent pas, ou que très succincte- ment, l’analyse asymptotique des algorithmes numériques : il est alors générale- ment difficile de comparer différentes méthodes pour résoudre un même pro- blème. Inversement, les ouvrages d’algorithmique à vocation pédagogique et qui abordent l’analyse asymptotique des algorithmes ne traitent pas de calcul numérique, c’est-à-dire avec des nombres à virgule flottante pour approcher les nombres réels : les algorithmes qui y sont présentés concernent alors plutôt les tris de tableaux ou de listes, ou encore des structures de données de type arbres pour représenter des dictionnaires.

La présentation des notions de l’analyse asymptotique des algorithmes de ce chapitre est volontairement réduite au strict minimum. Le lecteur qui sou- haiterait plus de détails est invité à se reporter à des ouvrages classiques. Une présentation récente et assez complète de l’algorithmique et de l’analyse asymp- totique est l’ouvrage de Cormen, Leiserson, Rivest et Stein [15] : le lecteur y trouvera également au chapitre 30 quelques algorithmes numériques comme la transformée de Fourier rapide et le calcul avec des polynômes. Knuth [33], plus ancien, reste une référence classique. Dans un second volume [34], l’auteur introduit l’analyse des calculs en virgule flottante ainsi que le calcul avec des polynômes. Enfin, Boissonat et Yvinec [8] proposent une analyse asymptotique de nombreux algorithmes pour la géométrie et les maillages.

L’algorithme de Strassen pour la multiplication matricielle a été proposé en 1969 [54] : il a suscité un grand émoi lors de sa publication. Jusque là, peu de gens pensaient qu’il pût exister un algorithme asymptotiquement plus rapide que la procédure en trois boucles imbriquées et avec un temps de calcul en Θ(n3). Strassen a proposé un algorithme en Θ(nlog2 (7)), avec log2(7)    2.81  <  3 et a montré que l’inversion d’une matrice n’est également pas une opération en Θ(n3) mais est dominé par le coût asymptotique de la multiplication de deux matrices. L’algorithme original de Strassen a été depuis amélioré a plusieurs reprise : l’algorithme le plus performant actuellement est celui proposé par Coppersmith et Winograd, en Θ(nω) avec ω 2.376, sachant que la meilleure borne inférieure connue pour le produit de deux matrices n’est autre que Θ(n2), car il faut bien remplir les n2 éléments du résultat. Le lecteur intéressé pourra lire l’article [44] qui fait une synthèse récente des différentes avancées sur le produit de deux matrices, ainsi que des commentaires dans [15] au chapitre 4 ou encore dans [41], paragraphe 2.11 pour l’inversion des matrices.

Dans la pratique, la constante cachée dans la notation Θ(.) reste relativement grande, même pour des implémentations très optimisées, et la procédure clas- sique en Θ(n3) est plus rapide pour des n suffisamment petits, disons n < n0. Dans les implémentations de la méthode essayent de déterminer ce n0, appelé point d’équilibre, et qui peut être finalement relativement grand, suivant le type de machine considéré. Dans les implémentations actuelles de la librairie d’algèbre linéaire atlas, ce n0 peut être évalué dynamiquement [16], par un étalonnage sur une machine spécifique lorsque la librairie est installée. Ensuite, à chaque appel et pour un n donné, la procédure la plus rapide est alors sélec- tionnée.

L’exercice 4 a présenté des optimisations possibles pour des matrices ayant une structure bande : ce type d’optimisation est disponible dans  la  li-  brairie lapack [2], écrite en Fortran, et qui dispose d’interfaces pour d’autres langages tels que le C ou le C++.

Aux chapitres suivants, et tout au long de cet ouvrage, nous aborderons des méthodes radicalement différentes de celle présentée ici. Ces méthodes nous permettront de diminuer encore de façon très importante le coût asympto- tique de la résolution d’un système linéaire. Elles nous permettrons de passer en dessous de la barre Θ(n2) pour la résolution d’un système matriciel. Elles s’appliqueront non plus à une matrice générale mais à des classes de matrices particulières. Les classes de matrices particulières qui seront considérées tout au long de cet ouvrage apparaissant très fréquemment en calcul scientifique. Elles s’appliquent efficacement à de nombreux problèmes, telles que la résolu- tions d’équations discrétisées par des méthodes différences ou éléments finis, et qui seront présentées dans les deux chapitres suivants du livre.

...

Chapitre 2 Transformée de Fourier et applications

L’objectif de ce second chapitre est d’introduire l’algorithme de la transformée de Fourier rapide et d’en présenter des applications à la méthode des différences finies pour résoudre des équations aux dérivées partielles. La première section présente l’algorithme à travers l’étude d’un problème simple, la déformation d’une corde de guitare. La méthode des différences finies est introduite dans une seconde section, pour le cas des problèmes en dimension un d’espace. En fin de section cette méthode est généralisé au cas de subdivisions non uniformes via la méthode des éléments finis. Nous considérons ici que le lecteur possède déjà quelques notions concernant les équations aux dérivées partielles. L’ex- cellent ouvrage de Raviart et Thomas [42] pourra être consulté sur ce sujet. La troisième section présente la résolution par transformé de Fourier des prob- lèmes en dimension deux d’espace et discrétisés par différences finies. Le cas de la dimension trois est traité en exercice.

2.1          Transformée de Fourier

Les transformées de Fourier servent intensivement dans le traitement du sig- nal. Le signal est défini dans un domaine de temps, à savoir une fonction qui à un temps associe une amplitude. L’analyse de Fourier permet de caractériser le signal en terme de fréquences, sous la forme de modes et de poids associés. Parmi les nombreuses applications récentes de la transformée de Fourier, citons les techniques de compression audio et vidéo, dont les populaires formats de fichier mp3 et jpeg. De nouvelles applications font leur apparition, comme le produit de deux polynômes ou celui de deux nombres entiers ou à virgule flot- tante en précision élevée, et dont les notes en fin de section donnent un aperçu.

Dans cet esprit, la section 2.3 présentera une autre application à la résolu- tion très rapide de solutions approchées d’équations aux dérivées partielles en dimension deux d’espace.

2.1.1      La transformée de Fourier discrète

Soit f une fonction de R dans C. La fonction f est supposée périodique de période T > 0 et est appelée signal. Ce signal se développe en série de Fourier de la manière suivante :

 f (x) = Σ cp e2iπpx/T avec cp = 1              ∫ T

f (x) e−2iπpx/T dx           (2.1)

p∈Z        0

Supposons qu’on ne connaisse qu’un nombre fini n de valeurs de f sur [0, T [ en des points régulièrement espacés : (f (kT /n))0™k™n−1 et qu’on appelle échan- tillon. Nous souhaitons déduire de cette information une approximation des coefficients de Fourier (cp)p∈Z de f :

f              kT           =             c n                                         p

p∈Z

e2iπpk/n =         cp

p∈Z

 wpk, 0 ™ k ™ n − 1,

où nous avons noté wn = e2iπ/n, appelé racine n-ième de l’unité, car il vérifie wn = 1. Ayant n données, il est logique de ne chercher à calculer que n coeffi- cients cp, p Z. De plus, les coefficients cp tendent vers zéro quand p devient grand, aussi nous chercherons à les approcher pour p n/2, . . . n/2 1 en supposant n pair. Dans le cas où n est impair, nous pourrions prendre un in- (n)p et caractérisés par :

. kT Σ −n/2™p™n/2−1 n/2−1

les coefficients de la série tronquée

Il s’agit d’un système de n équations à n inconnues, et nous allons voir que  ce système admet une solution unique. Pour une raison de commodité, réar- rangeons ce système en numérotant les inconnues de 0 à n 1 de la façon suivante :

 Fp = (n)p c(n)

p−n

si 0 ™ p ™ n/2 − 1

si n/2 ™ p ™ n − 1

et posons fk = f (kT /n). Le système devient :

n−1

n

p=0

Notons A la matrice du système. Elle est symétrique :

 1   1    1             . . .   1    

A =   1    w2    4             . . .    w2(n−1)   

   ...

 ...

 ...

...



1 wn−1

 w2(n−1)

 . . . w(n−1)2

Par un algorithme direct classique de type Gauss ou Cholesky, la résolution du système linéaire serait en Θ(n3) opérations. L’intérêt de ce système vient du fait que l’inverse de A admet une forme explicite connue :

 1   1    1             . . .   1    

1    wn−1

A−1  =   1    wn−2

 wn−2

wn−4

 . . .    wn−n+1

. . .    wn−2(n−1)

 

   ...

 ...

 ...

...

2   

2(n  1)   −(n−1)

1    wn−n+1        wn−       −             . . .    wn

Par conséquent, le système se résout explicitement :

n−1

F  = 1 Σ w−pkf ,  0 ™ k ™ n − 1.    (2.2)

Par un produit matrice-vecteur avec A−1, le calcul de (Fk)0™k™n−1 reste cepen- dant en Θ(n2) opérations. Bien que ce soit déjà bien plus compétitif que l’algo- rithme de Strassen vu précédemment, nous allons voir qu’il est possible de faire ici encore beaucoup mieux en exploitant la forme particulière de la matrice A et en faisant des hypothèses sur la parité de n.

2.1.2      Algorithme par bisection

Supposons que n est pair, et dans (2.2), séparons les termes d’indices pairs et impairs :

k             2             k             k Σ          − 1,

F   =  1   P + w−kI              , 0 ™ k ™ n

où Pk et Ik représentent la somme des indices pairs et impairs, respectivement :

P             =             2 F

k             n             0

I              =             2 F

k             n             1

+ wn−2k F2

+ wn−2k F3

+ . . . + wn−(n−2)k F

+ . . . + wn−(n−2)k F

 n−2Σ n−1Σ

 Supposons ensuite n divisible par 4. Ces sommes peuvent alors toutes deux  à leur tour être séparées en termes pairs et impairs. Supposant enfin n = 2r, nous pouvons poursuivre la dichotomie jusqu’à ne plus avoir qu’un seul terme dans la somme. Nous obtenons alors un algorithme récursif appelé transformée de Fourier  rapide. L’idée de base est similaire au autres algorithmes étudiés  au chapire précédent, comme le calcul de xn ou le produit de deux matrices (algorithme de Strassen) : il s’agit de découper le calcul en plusieurs sous- calculs de taille plus petite, ceci récursivement et jusqu’à obtenir des calculs très simples. Le découpage se fait ici, comme pour le calcul de xn, en deux problèmes de taille moitié. L’algorithme présente l’intérêt de se programmer simplement, en suivant exactement la démarche choisie.

fft.h

Un programme d’appel de cette fonction est :

fft.cc

Ce programme lit simplement les données à partir  de  l’entrée standard, puis appelle la transformée de Fourier rapide, et enfin écrit le résultat sur la sortie standard. Pour illustrer la transformée de Fourier rapide, calculons les amplitudes des harmoniques après un pincement de corde de guitare : nous pinçons la corde en son milieu pour avoir un beau son, bien propre.

guitar.cc

Observons maintenant le résultat du programme ftt avec un échantil-  lonage de 16 ou 32 valeurs pour le pincement de la corde de guitare :

make fft guitar

./guitar 16 | ./fft

./guitar 32 | ./fft

La première valeur F (n) est l’amplitude de la déformation de la corde, qui vaut 1/4 ici, indépendamment de n. Les autres valeurs de F (n), k “ 1, correspondent aux amplitudes des harmoniques successives de la corde de guitare. La quantité

F (n) correspondant à l’approximation de F1, est l’amplitude de la note fonda- mentale de la corde. Cette quantité est réelle. Pour ce problème particulier, la valeur exacte F1 peut être calculée directement (voir exercice 9). Nous pouvons donc calculer l’erreur commise F (n) F1 en tronquant la somme aux n pre- miers termes. Le calcul de cette amplitude est d’autant plus précis que n est grand. La pente sur la Fig. 2.2, en axes log-log, suggère que l’erreur commise

…

Figure 2.1 – Déformation de la corde de guitare et décroissance correspondante des amplitudes des harmoniques.