L'ALGORITHME CORDIC Christophe rouen fr

Christophe Devalland

Description

Niveau: Secondaire, Lycée, Première

• mémoire

1 L'ALGORITHME CORDIC () Utiliser sa calculatrice pour déterminer la valeur d'un cosinus, d'un logarithme ou d'une racine carrée est un geste devenu tellement banal que plus personne ne se demande pourquoi et comment cela marche. Mais, pour reprendre une formule bien connue, a-t-on vraiment besoin de soulever le capot de sa voiture et de comprendre le fonctionnement du moteur pour s'en servir ? Evidemment non. Pourtant, un jour en classe de première, un élève de ma femme un peu curieux a été fier d'expliquer à ses camarades comment on pouvait extraire une racine carrée à la main. A cette occasion, il s'est ensuite demandé si la calculatrice utilisait la même méthode pour calculer une racine carrée. Nous nous sommes également posé cette question sans pouvoir y répondre. On pourrait même aller plus loin en se demandant comment la calculatrice peut donner la valeur d'un sinus ou d'un logarithme avec autant de précision. Quand il s'agit d'opérations simples telles que l'addition ou la multiplication, on peut assez bien imaginer comment elle procède parce que nous savons faire ces opérations à la main. Mais pour les autres fonctions, comment fait-elle réellement ? Utilise-t-elle des développements limités, des approximations de fonctions, ou d'autres mécanismes plus complexes ? C'est en fouillant à droite et à gauche que j'ai trouvé quelques réponses qui m'ont servies de base pour cet article.

• ?1


• besoins de calculs en temps réel


• ?? ??


• matrice de la rotation vectorielle d'angle ?


• ?n


• coordinate rotation


• série ∑


• angle

Sujets

• Première
• Mathématiques
• Mémoire
• Rotation (mathematics)

Informations

• Publié par: apmep
• Nombre de lectures: 40
• Langue: Français

Extrait

L’ALGORITHME CORDIC ([e-mail]) Utiliser sa calculatrice pour déterminer la valeur d’un cosinus, d’un logarithme ou d’une racine carrée est un geste devenu tellement banal que plus personne ne se demande pourquoi et comment cela marche. Mais, pour reprendre une formule bien connue, a-t-on vraiment besoin de soulever le capot de sa voiture et de comprendre le fonctionnement du moteur pour s’en servir ? Evidemment non. Pourtant, un jour en classe de première, un élève de ma femme un peu curieux a été fier d’expliquer à ses camarades comment on pouvait extraire une racine carrée à la main. A cette occasion, il s’est ensuite demandé si la calculatrice utilisait la même méthode pour calculer une racine carrée. Nous nous sommes également posé cette question sans pouvoir y répondre. On pourrait même aller plus loin en se demandant comment la calculatrice peut donner la valeur d’un sinus ou d’un logarithme avec autant de précision. Quand il s’agit d’opérations simples telles que l’addition ou la multiplication, on peut assez bien imaginer comment elle procède parce que nous savons faire ces opérations à la main. Mais pour les autres fonctions, comment fait-elle réellement ? Utilise-t-elle des développements limités, des approximations de fonctions, ou d’autres mécanismes plus complexes ? C’est en fouillant à droite et à gauche que j’ai trouvé quelques réponses qui m’ont servies de base pour cet article. Je me suis par ailleurs amusé à programmer certains algorithmes sur ma calculatrice pour vérifier que, malgré leur simplicité, ils donnent de bons résultats.

Un peu d’histoire  A la fin des années 50, les calculateurs électroniques sont en plein essor. Les domaines où ils interviennent sont de plus en plus variés et on les trouve par exemple dans les avions où ils servent d’assistants à la navigation aérienne. De ce fait, les besoins de calculs en temps réel se font de plus en plus sentir, notamment pour les calculs trigonométriques. Ainsi, en 1959, Jack E. Volder met au point un algorithme qui permet d’approximer des fonctions trigonométriques à partir d’opérations élémentaires (additions, soustractions et multiplications). Cet algorithme appelé « algorithme CORDIC » (pour Coordinate Rotation DIgital Computer) repose, comme son nom l’indique, sur le calcul des coordonnées de vecteurs auxquels on applique une rotation bien choisie. Très vite, grâce à sa mise en œuvre matérielle très simple (l’ingéniosité de l’algorithme permet en effet un câblage électronique extrêmement simple), l’algorithme CORDIC est repris dans des domaines très diverses tels que : le traitement de signaux radars, les coprocesseurs mathématiques (intel i8087 par exemple) ou les calculatrices scientifiques (toutes marques confondues). Un autre avantage non négligeable de cet algorithme est qu’il permet d’obtenir une précision déterminée à l’avance en effectuant un nombre d’itération donné. Mais pour être certain du résultat, il faut restreindre l’intervalle d’utilisation au domaine de convergence de l’algorithme, comme on le verra plus bas (il existe toutefois des techniques permettant de contourner cette limitation).

Première approche J’ai repris ici un article de Yves Suprin paru dans le bulletin n°10 de novembre 1996 de la ère régionale de Haute−Normandie. Il s’agit en fait d’un devoir à la maison proposé en classe de 1 S dont l’objectif est de permettre le calcul d’une valeur approchée du sinus d’un anglea appartenant à l’intervalle [0 ;π/2[ en s’inspirant de l’algorithme CORDIC. L’idée est la suivante :apeut se décomposer en une sommea 1 +a 2 +a 3 +...+a n . Les anglesa i étant de plus en plus petits, et suffisamment, pour assurer la convergence ; un anglea i  pouvant être utilisé plusieurs fois dans cette décomposition. On converge alors vers l’anglea de la façon suivante :

1

y

o

a

A

A 4

A 3

A 2

A 1

x

aveca=;OA)A A ) et (OA ) i (OA i i+1 i ⊥( i i+1 En posant : OA i =R i ;ω i =tana i eta=a 1 +a 2 +...+a n , on obtient les expressions : x n+1 y n+1 y n+1 cosa= , sina= et tana=. R n+1 R n+1 x n+1 On démontre ensuite les relations :  x=x− ω y  i+1i i i {1;...;n}y x yPouri∈  i = ω i i + i +1  2 R=R1+ ω  i+1i i Il suffit donc de connaîtrex 1 ,y 1 ,R 1  et la suite des ω i  pour pouvoir calculerx n ,y n ,R n ainsi en et déduire cosaet sina. -i Pour les calculs, on peut choisir les anglesa i dans l’intervalle [0 ;π/2[ tels que tana i =10 . Cela -i permettra des calculs assez simples dans les relationspuisqu’on aura alors ω i =10 . Les anglesa i déterminent une fois pour toute avec la calculatrice pour se iÂ4 et peuvent être -i confondus avec 10 pouriÃ5 (c’est d’ailleurs le but de la première partie du devoir que de -5 -14 démontrer que lorsque 0Âa i Â10 , tana i =a i à10près. On ne sera donc pas pénalisé par cette -8 approximation si l’on veut une valeur approchée de sinaet cosaà10près par exemple.). -i Pour faire tourner l’algorithme, il suffit donc de donner à la calculatrice la table des arctan(10 ), ces nombres étant des constantes prédéfinies dans le programme. Puis, pour un angle donnéa, on soustrait autant de fois qu’il est possible les anglesa,a,a, etc... jusqu’à ce que le reste soit 1 2 3 inférieur à la précision désirée. -8 Voici une version possible de cet algorithme pour la TI 92 Plus (précision 10 ) : Cordic1(a) Func Local angle,x,y,z,r,i {0.785398163,0.099668652,0.009999667,0.001,0.0001,1–ª005,1–ª006,1–ª007,1– ª008}»angle 1»x 0»y 1»r For i,0,8 While a>=angle[i+1]  x-y/10^i»z

2

 y+x/10^i»y  z»x  r*§(1+10^(ª2*i))»r  a-angle[i+1]»a EndWhile EndFor Return {x/r,y/r} EndFunc Il faut toutefois noter que cet algorithme n’est pas à proprement parler l’algorithme CORDIC. En effet, s’il en reprend le principe, il n’en a pas la simplicité originelle. Il faut se souvenir que dans les années 50, il aurait certainement été très pénalisant en temps d’exécution de réaliser un système électronique capable de réaliser une boucle « While » reposant sur la comparaison de deux nombres. De plus, le calcul d’une racine carrée à chaque itération est en contradiction avec les contraintes liées à l’utilisation d’opérations élémentaires. Mais puisque nous avons aujourd’hui des moyens plus performants à notre disposition, utilisons les et l’avantage de cette -i méthode est sa rapidité de convergence : avec 8 constantes ( les arctan(10 ) ) on obtient une -8 précision de l’ordre de 10 . Ce n’est pas le cas avec l’algorithme CORDIC car il en faut beaucoup plus, comme nous allons le voir. Ce que l’on gagne en simplicité, on le perd en rapidité de convergence.

L’algorithme CORDIC Afin de respecter les contraintes matérielles de l’époque et pour pouvoir effectuer des calculs en temps réel, il fallait donc trouver un algorithme n’utilisant pas de test (ou seulement une comparaison avec zéro, très facile à mettre en œuvre) et dont la durée d’exécution soit toujours constante. Ce n’est pas le cas avec l’algorithme précédent puisqu’on ne sait pas à l’avance combien de passages dans la boucle While seront nécessaires pour se rapprocher de l’angle de départ. La question que l’on peut se poser est la suivante : peut-on approcher un angleθune avec précision donnée à l’avance en le décomposant en une somme d’angles prédéterminés ? En d’autres termesθpeut-il s’écrire : n θó δ k ε k ∑ k=0 s et où ε n )pouvant prendre comme valeur que –1es ne (t une suite d’angles distincts prédéfini δ k ou +1 ? Plus précisément, comme l’on veut atteindre une précision donnée, peut-on écrireθ comme une combinaison d’angles ε tels que : k n ? θ − δ k ε k ≤ ε n ∑ k=0  approcher (Si tel était le cas, on serait certain d’effectuern+1exactement pour étapes θ à ε n près). Essayons de répondre à cette question. D’abord, on voud e c rait que la séri δ n ε n onverge versθ conformément à l’inégalité, ce qui impose àε n de tendre vers 0. Plus simplement nous utiliserons la condition : ( ε )est une suite de réels positifs décroissant vers 0 n Dans la pratique, nous choisirons même une s érie soit absolument uite( ε )telle que la s δ n ε n n convergente, c e que la ’est-à-dire tell série ε n converge. La convergence de la sérieθest-elle touj δ n ε n vers ours possible quelque soit le réelθ? En fait, non. ne pouvant prendre comme valeur que –1 ou +1 , si > ou < alors on voit δ k θ ε n θ- ε n bien qu’il sera impossible d’approcher l’angle avec la série θ δ n ε n :

θ

+∞ ∑ k ε k=0

3

Plus simplement, pour être certain de pouvoir approcher l’angleθà ε près il faut que : n n ou en se limitant θ ≤ ε n θ ≤ ε k auxn+1premiers termes. ∑ k=0 Voyons maintenant comment choisir la suite (δ n ). i−1 Posonssde la façon suivante : ∑ i = δ k ε k et construisons la suite (s n ) k=0 SiθÃ0 alorss 1 = ε 0 (on se rapproche deθen partant de ε 0 ; c’est-à-dire :δ 0 =+1)  sinons 1 =- ε 0 (on se rapproche deθen partant de - ε 0 ; c’est-à-dire :δ 0 =- 1) Siθ-δ 0 ε 0 Ã0 alorss 2 = ε 1 + ε 1  (on se rapproche deθajoutant en ε 1 ; ieδ 1 =+1)  sinons 2 = ε 1 - ε 1 (on se rapproche deθen retranchant ε 1 ; ieδ 1 =-1) etc... Plus généralement : - ε k ε k θ θ s k s k Sis k est supérieure à la valeur deθalors Sis k est inférieure à la valeur deθalors s k+1 =s k - ε k doncδ k =-1s k+1 =s k + ε k doncδ k =+1 Finalement : δ k =signe(θ-s k ) Il reste maintenant à étudier si l’inégalitéimplique certaines conditions sur la suite (ε n ). Supposons que cette inégalité soit vraie pour l’entiernet voyons quelle condition doit vérifierε n+1 pour qu’elle reste vraie pour l’entiern+1. n+1n θ ∑ δ k ε k θ ∑ δ k ε k δ n+1 ε n+1 On a, au rangn+1:− −− = k=0k=0 n 0 alors : Si θ − δ k ε k à ∑ k=0 n comme on suppose l’inégalitévraie au rangn: 0 θ − δ ε Âε n . ∑ k k k=0 n ∑ k k+1  - . De plusδ n+1 =signe(θ-s n+1 )=+1 donc : -ε n+1  θ − δ ε − ε n ε n ε n+1 k=0 n+1 ∑ Âmax(ε n+1 ;ε n -ε n+1 ). D’où : θ − δ k ε k k=0 n+1 ∑  Comme on veut avoir θ − δ k ε k ε n+1 , cela équivaut àε n -ε n+1 Âε n+1 . k=0 1 C’est-à-direε n+1 Ãε n . 2 n i θ − k Â0 ∑ je vous laisse le soin de vérifier que l’on trouve la même condition sur . S δ k ε ε n+1 k=0 Finalement, pour que l’inégalitévraie pour tout soit n∈É, il faut que la suite (ε n ) respecte la condition nécessaire et suffisante : 1  tnÉε n ε n+1 Âε n pour tou∈ 2  4

Sous cette condition, il est simple de démontrer par récurrence quevraie pour tout entier est naturel, sous réserve que respecte la condition de départ : . θ θ − δ 0 ε 0 ≤ ε 0 Il est maintenant possible de répondre à notre question : En respectant les conditions à, il est possible d’approcher l’angleθ àε n en faisant un près nombre déterminé d’additions ou de soustractions d’anglesε k prédéfinis. En revenant à l’exemple précédent, on peut remarquer que la suite (ε n ) définie par -k ε k pourarctan(10 ) k∈Éne vérifie pasla condition, ce qui oblige à répéter autant de fois que = possible la soustraction de l’angleε k (oua k pour respecter les notations). -i Cependant, l’intérêt d’avoir travaillé avecω i =10 est d’obtenir des calculs très simples dans la i relationpuisque diviser par 10 est effectivement très simple pour nous (humains). Mais il faut penser qu’un calculateur électronique ne compte bien souvent qu’en base 2. Il est donc beaucoup plus facile pour lui de diviser par des puissances de 2. Par exemple : 3 2 1 0 ce qui donne en base 2 :14 s’écrit 1×2 +1×2 +1×2 +0×2 14 10 =1110 2 2 1 0 14÷2=7 s’écrit 1×2 +1×2 +1×2 ce qui donne en base 2 : 7 10 =111 2 (un décalage à droite) 4 3 2 1 0 ce qui donne en base 2 :14×2=28 s’écrit 1×2 +1×2 +1×2 +0×2 +0×2 28 10 =11100 2 (un décalage à gauche) Tout comme nous le faisons en base 10, multiplier ou diviser par 2 en base 2 revient à décaler à droite ou à gauche les chiffres du nombre. Cette propriété est très intéressante car tout système électronique un peu élaboré est capable d’effectuer ce genre d’opération sur des cases mémoires, appelées aussi des registres. -k La base 2 est d’autant plus intéressante que si l’on considère les anglesε k pour=arctan(2 ) k∈É, on remarque qu’ils vérifient la condition. (Je laisse le soin au lecteur de faire l’étude fastidieuse 1 1 1 de la fonctionx֏arctan( )−arctan( )montrer qu’elle est bien toujours positive). pour x+1x 2 2 2 −(n+1) arctan(2 ) 1 Constatons expérimentalement que≥pourn∈É : −n 2 arctan(2 )

Ca marche en base 2

Ca ne marche pas en base 10

Application au calcul du sinus d’un réel Reprenons un réelθque l’on veut approcher par une combinaison dennombresε k de plus en plus petits. Ainsiθóδ 0 ε 0 +δ 1 ε 1 +...+δ n-1 ε n-1 oùδ k =±1. Considérons maintenant la matrice de la rotation vectorielle d’angleθobtenir la relation pour suivante entre les coordonnées des vecteursu Å et v Å.  x ′  v   y ′  x ′   cos θ −sin θ  x  =  y ′   sin θ cos θ  y       θ  x  u   y

5

ce qui revient à : cos sin cos sin cos sin  x ′   ε 0 − ε 0  ε 1 − ε 1   ε n−1 − ε n−1  x  =...   y sin cos y  ′  sin 0 cos 0  sin 1 cos 1   n−1 ε n−1     εε ε ε ε        L’idée maintenant est de faire apparaître tanε k en factorisant par cos ε k chaque matrice : n−1  x ′   1−tan ε k  x  =cos ε k  y ′  ∏   y    tan ε k 1   k=0   -k En prenantε k =arctan(2 ) et en définissantδ k =signe(θ-(δ 0 ε 0 +δ 1 ε 1 +...+δ k-1 ε k-1 )) on obtient : n−1−k  1 2   x ′  − δ k  x  =cos ε  y ′  ∏ δ 2  y  k  −k      1 k=  k  0 Compte tenu du choix des suites (ε n ) et (δ n ), on est assuré d’obtenir l’égalité : n −k , à condition que respecte la condition c’est−à−dire que Âarctan(2 )ó θ =lim δ k ε k θ θ ∑ n→+∞ k=0 1,7. Mais puisque les fonctions sinus et cosinus sont 2π−périodiques, on pourra se limiter à π θ ≤ 2 ∞ −k Ainsi, en appelant K le réelcos(arctan(2 )) (qui est une constante ó 0,60725), on a la ∏ k=0 relation : ∞ −k∞ −k  1 2   1 2   cos θ  − δ k  1  − δ k  K  =K=.  −k   −k   sin θ  ∏  0  ∏  0        δ 2 1 δ 2 1 =  k  =  k  k0k0 Cette étude préliminaire permet d’en déduire les relations de l’algorithme CORDIC définissant les suites (x n )et(y n )convergeant vers cosθet sinθ:  x=K;y=0 ;z= θ 0 0 0  δ k =signe(z k )  −k-k x y  x k+ = k − δ k k 2oùε k =arctan(2 ) 1 −k  xy y k+ = k + δ k k 2 1  z z k1 = k − δ k ε k  + En notant A(cosθ;sinθ) et A i (x i ;y i ), on obtient la suite de points suivant : y

o

θ

A 1

A 0

A 3

A 6 A 5 A 4 A 2

x

6

Voici ci-dessous un programme qui donne une valeur approchée de cosθ et sinθutilisant 15 en -14 -5 constantesε k . Cela permet une précision de l’ordre de 2 . C’est beaucoup moinssoit ó 6×10 -8 que la précision de 10 obtenue avec 9 constantes grâce à l’algorithme précédent en base 10 mais l’exécution est beaucoup plus rapide. Cordic2(z) Func Local angle,x,y,temp,s,i {.7854,.46365,.24498,.12435,.06242,.03124,.01562,.00781,.00391,.00195,9.8– ª4,4.9–ª4,2.4–ª4,1.2–ª4,6.–ª5}»angle .60725»x 0»y For i,0,14  when(z=0,1,signe(z))»s  x-s*y/2^i»temp  y+s*x/2^i»y  temp»x  z-s*angle[i+1]»z EndFor Return {x,y} EndFunc Notons que cet algorithme est « réversible ». En reprenant le systèmeen modifiant les et conditions suivantes : Siz 0 =0etδ k =-signe(y k ) alorsz n óarctan(y 0 /x 0 ). La preuve repose sur les égalités suivantes : y k + δ tan ε k k y x y k+1k k = =tan(arctan+ δ ε )(d’après une formule trigonométrique bien connue !) k k x y x k+1k k 1− δ tan ε k k x k y y k+1k d’où :arctan( )=arctan( )+ δ ε . k k x x k+1k n−1 y y n0 ( δ i ε i Et par une récurrence immédiate :arctan )=arctan( )+ ∑ x x n i = 0 0 y k plus,( ) (arctan( ))anger pour quexsoit toujours De δ k = −signe y k = −signe; il faudra donc s’arr k x k positif. k−1k−1 y y 0 0 ∑ ∑ (arctan( ) )signe( ( On en déduit : δ k = −signe+ δ i ε i = −arctan )− δ i ε i ) . x x 0 0 i=0i=0 Ainsi, comme les propriétés,etsont vérifiées, on a : n−1n−1 y y 00 arctan( ) , c’est−à lim arctan( ) . − − δ ε ≤ ε −dire δ i ε i = − ∑∑ i i n x x n→+∞ 0 0 i=0i=0 n−1n−1 y 0 ∑ ∑ Or, par construction, z z . La suite (z) converge donc bien vers arctan( ) . n = 0 − δ i ε i = − δ i ε i n x 0 i=0i=0 -5 Voici ci−dessous un programme qui donne une valeur approchée de arctan(y/x) à 6×10 près. Cordic3(x,y) Func Local angle,temp,z,s,i {.7854,.46365,.24498,.12435,.06242,.03124, .01562,.00781,.00391,.00195,9.8–ª4,4.9– ª4,2.4–ª4,1.2–ª4,6.–ª5}»angle 0»z For i,0,14  when(y=0,ª1,ªsigne(y))»s

 x-s*y/2^i»temp  y+s*x/2^i»y  temp»x  z-s*angle[i+1]»z EndFor Return z EndFunc

Extension à d’autres fonctions En 1971, John Walther a montré qu’il était possible d’étendre l’algorithme CORDIC à d’autres fonctions. En effet, en modifiant un peu le système (8), on arrive à exprimer les fonctions hyperboliques, les fonctions logarithme et exponentielle ainsi que la multiplication et la division ! Voici la généralisation de cet algorithme : −k  x=x−m δ y2 k+1kk k  −k y=y+x2oùm=0ou±1; k sont des constantes prédéfinies ; k =±1.  k+1k δ k k ε δ  z=z− δ ε k+1kk k  En choisissantm=1, on retrouve le cas des fonctions trigonométriques étudiées plus haut.

Cas des fonctions hyperboliques (m=−1) Comme pour les fonctions trigonométriques, on peut montrer qu’en prenant : +∞ -k−i ε k =argth(2 ) ;δ k =signe(z k );x 0 =ch(argth(2 ));y 0 =0 ;z 0 =θ ∏ i=0 alorsx n óchθety n óshθ. θ On obtient par conséquent quex n +y n óe. Cordic4(z) Func Local angle,x,y,temp,s,i {0.54931,0.25541,0.12566,0.06258,0.03126,0.01563,0.00781,0.00391,0.00195,0. 00097,0.00049,0.00024,0.00012,6–ª005}»angle 1.20513»x 0»y For i,1,14  when(z=0,1,sign(z))»s  x+s*y/2^i»temp  y+s*x/2^i»y  temp»x  z-s*angle[i]»z EndFor Return {x,y} EndFunc Faisons quelques remarques importantes pour ce programme : D’abord, il faut noter que les résultats fournis ne seront proches de chzet shzqu’à la condition où +∞ −i z<arg th(2 ) (soit environ 1,05). On peut alors se dire que cet algorithme n’a pas beaucoup ∑ i=1 d’intérêt puisque les fonctions ch et sh ne sont pas périodiques comme l’étaient les fonctions sinus et cosinus : on n’aura donc jamais accès avec ce programme à des valeurs de chzet shzavecz>2 par exemple. Mais il est possible de contourner ce problème. M. Vern m’a donné quelques indices : " quelque soit le type de développement ou d'algorithme que l'on utilise, il faut en général restreindre le domaine de calcul. Pour arctan, par exemple, on utilise : π 1x−k arctan(x)=−arctan  ou bien arctan(x)=arctan(k)+arctan  , 2  x   1+kx  pour ln,

8