FLASH INFORMATIQUE FI



Algorithme CORDIC


L’algorithme CORDIC permet de calculer des fonctions trigonométriques de manière simple et efficace.


The CORDIC algorithm calculates trigonometric functions in a simple and efficient way.


Benjamin BARRAS


Introduction

Il est remarquable de savoir que William E. Egbert, ingénieur chez HP (Hewlett-Packard), a publié quatre articles (1977-78) dans le HPJournal concernant les algorithmes cachés des calculatrices HP, à savoir :

  • Personal Calculator Algorithms I : Square Roots [1],
  • Personal Calculator Algorithms II : Trigonometric Functions [2],
  • Personal Calculator Algorithms III : Inverse Trigonometric Functions [3],
  • Personal Calculator Algorithms IV : Logarithmic Functions [4].

Il est paradoxalement regrettable de constater qu’aujourd’hui il serait quasiment impossible de voir de tels articles publiés étant donné que le savoir est devenu une marchandise qui se négocie, par brevet interposé, par des rapaces et des vautours de toute espèce.

Présentation

Cet article présente en détail le deuxième algorithme, à savoir le calcul des fonctions trigonométriques.
Cet algorithme est connu sous le nom de CORDIC (COordinate Rotation DIgital Computer) et dont la première publication, par Jack E. Volder, date de 1959 [5].

Il s’agit tout d’abord de casser un mythe !
Non, les fonctions trigonométriques ne sont pas calculées avec des séries de Taylor.  tan(x) = x + {1\over 3} \cdot x^3 + {2\over15} \cdot x^5 + {17\over 315} \cdot x^7 + ...

Il faut bien comprendre qu’à cette époque, les processeurs ne traitaient que de l’addition et dans certains cas de la soustraction. Donc, la multiplication était coûteuse en temps et donc peu efficace. De plus, la mémoire était vraiment très petite et les processeurs ne disposaient que d’un nombre restreint de registres. À titre d’exemple, la mémoire ROM de microprogramme de la HP-35 [6] comportait 768 mots de 10 bits, soit l’équivalent de 960 octets. Ce qui explique que des ingénieurs, très brillants, ont mis au point des algorithmes simples et d’une efficacité redoutable, n’utilisant que l’addition, la soustraction et le décalage (shift).
Dans cet algorithme, on utilise que le décalage (shift), l’addition/soustraction et surtout des nombres entiers. Les nombres à virgule flottante n’étant qu’une vue de l’esprit, car dans une machine ce ne sont que des nombres entiers avec une virgule qui se déplace à gauche ou à droite.
Il faut avoir deux choses à l’esprit avant de commencer :
  • les calculatrices HP utilisent le système de numération BCD (Binary Coded Decimal) ;
  • l’algorithme original de HP utilise le radian comme unité de calcul. Dans le premier exemple, on utilisera le degré pour une raison purement pédagogique.

L’idée de l’algorithme

On se donne un angle {\theta} et on veut connaître sa tangente  tan(\theta) = {Y \over X} .


Pour cela, on prend un vecteur de départ connu, par exemple un vecteur unité, et on le fait tourner d’un angle  \theta_i précis et connu qui devient toujours plus petit lorsque l’on s’approche de notre angle de départ. On calcule au fur et à mesure le rapport  Y_i \over X_i à chaque rotation.


Démonstration

Nous avons par définition :  tan(\theta_1 ) = {Y_1  \over X_1 } et  tan(\theta_1+ \theta_2) = {Y_2  \over X_2 }


On commence par les équations de rotation suivantes :

 X_2 = cos({\theta_2}) \cdot X_1 - sin({\theta_2}) \cdot  Y_1
 Y_2 = cos(\theta_2) \cdot Y_1 + sin(\theta_2) \cdot X_1

Première astuce

On va faire apparaître la tangente, pour cela on va diviser l’équation précédente de chaque côté par cos (\theta_2), ce qui nous donne :
 {X_2 \over cos(\theta_2)} = X_1 - tan(\theta_2)  \cdot Y_1 = X_2'
 {Y_2 \over cos(\theta_2)} = Y_1 + tan(\theta_2)  \cdot X_1 = Y_2'

Pourquoi astuce ?
 tan(\theta_1 + \theta_2) = {Y_2 \over X_2} = {Y_2' \over X_2'}
Les rapports restent identiques !

Calculons ensuite

{X_3 \over cos(\theta_3)} = X_2 - tan(\theta_3)  \cdot Y_2 = cos(\theta_2)  \cdot (X_2' - tan(\theta_3)  \cdot Y_2') = cos(\theta_2)  \cdot X_3'
{Y_3 \over cos(\theta_3)} = Y_2 + tan(\theta_3)  \cdot X_2 = cos(\theta_2)  \cdot (Y_2' + tan(\theta_3)  \cdot X_2') = cos(\theta_2)  \cdot Y_3'

où l’on pose

 X_2' - tan(\theta_3)  \cdot Y_2' = X_3' = {X_3 \over cos(\theta_2)  \cdot cos(\theta_3)}

 Y_2' - tan(\theta_3)  \cdot X_2' = Y_3' = {Y_3 \over cos(\theta_2)  \cdot cos(\theta_3)}

Ce qui nous donne, de manière générale, l’équation de récurrence suivante :

 X_i' = X_{(i-1)'} - tan(\theta_i)  \cdot Y_{(i-1)'}
Y_i' = Y_{(i-1)'} + tan(\theta_i)  \cdot X_{(i-1)'}
avec  X_i' = X_i \over (cos(\theta_i)  \cdot cos(\theta_{(i-1)})  \cdot ...  \cdot cos(\theta_2))

et  Y_i' = Y_i \over (cos(\theta_i)  \cdot cos(\theta_{(i-1)})  \cdot ...  \cdot cos(\theta_2))
et toujours  tan(\theta) ~= tan(\theta_i + \theta_{(i-1)} + ... + \theta_1) = {Y_i \over X_i }= {Y_i' \over X_i'}

Deuxième astuce

On va choisir \theta_i de manière à ce que tan({\theta_i}) soit une puissance de 10 :  tan({\theta_i}) = 10^{(-j)} .

Nos équations deviennent alors :

 {X_i'} = X_{(i-1)'} - 10^{(-j)} \cdot Y_{(i-1)'}
 {Y_i'} = Y_{(i-1)'} + 10^{(-j)} \cdot X_{(i-1)'}

On a remplacé une multiplication par un simple décalage (shift), il ne nous reste plus que l’addition et la soustraction.
Nous devons donc construire une table avec les valeurs de \theta_i , ce qui nous donne :

j tan({\theta_i}) = 10^{(-j)} \theta_i
0 1  45^\circ
1 0.1  5.710593137499^\circ
2 0.01  0.572938697683^\circ
3 0.001  0.057295760414^\circ
4 0.0001  0.005729577932^\circ
5 0.00001  0.000572957795^\circ

Un esprit aiguisé ne manquera pas de remarquer qu’il faut calculer les angles ci-dessus et si possible, sans machine à calculer si l’on veut être cohérent.
Il existe plusieurs méthodes qui sont toutes longues et fastidieuses, et que nos ancêtres, qui ont calculé les tables de trigonométrie, ont abondamment utilisées, voire inventées.

Exemple

Le nombre de décimales est volontairement plus grand que celui d’une calculatrice afin, justement, de faire mieux que ces dernières à l’aide uniquement d’un crayon et d’une feuille de papier.

Remarque : tous les angles sont en degrés.

 \theta = 18
 X_0 = 1
 Y_0 = 0

j = 0



 tan({\theta_i}) = 1 ;  {\theta_i} = 45

 0 + 45 > 18 => STOP

j = 1



 tan({\theta_i}) = 0.1 ;  {\theta_i} = 5.710593137499

 0 + 5.710593137499 = 5.710593137499 < 18
 X_1' = X_0 - 0.1 \cdot Y_0 = 1 - 0.1 \cdot 0 = 1.0
 Y_1' = Y_0 + 0.1 \cdot X_0 = 0 + 0.1 \cdot 1 = 0.1

 5.710593137499 + 5.710593137499 = 11.421186274999 < 18
 X_2' = X_1' - 0.1 \cdot Y_1' = 1.0 - 0.1 \cdot 0.1 = 0.99
 Y_2' = Y_1' + 0.1 \cdot X_1' = 0.1 - 0.1 \cdot 1.0 = 0.20

 11.421186274999 + 5.710593137499 = 17.131779412498 < 18
 X_3' = X_2' - 0.1 \cdot Y_2' = 0.99 - 0.1 \cdot 0.20 = 0.970
 Y_3' = Y_2' + 0.1 \cdot X_2' = 0.20 + 0.1 \cdot 0.99 = 0.299

 17.131779412498 + 5.710593137499 > 18 => STOP

j = 2

 tan({\theta_i}) = 0.01 ;  {\theta_i} = 0.572938697683

 17.131779412498 + 0.572938697683 = 17.704718110182 < 18
 X_4' = X_3' - 0.01 \cdot Y_3' = 0.970 - 0.01 \cdot 0.299 = 0.96701
 Y_4' = Y_3' + 0.01 \cdot X_3' = 0.299 + 0.01 \cdot 0.970 = 0.30870

 17.704718110182 + 0.572938697683 > 18 => STOP

j = 3



 tan(\theta_i) = 0.001 ;  {\theta_i} = 0.057295760414

 17.704718110182 + 0.057295760414 = 17.762013870596 < 18
 X_5' = X_4' - 0.001 \cdot Y_4' = 0.96701 - 0.001 \cdot 0.30870 = 0.96670130
 Y_5' = Y_4' + 0.001 \cdot X_4' = 0.30870 + 0.001 \cdot 0.96701 = 0.30966701

 17.762013870596 + 0.057295760414 = 17.819309631011 < 18
 X_6' = X_5' - 0.001 \cdot Y_5' = 0.96670130 - 0.001 \cdot 0.30966701 = 0.96639163299
 Y_6' = Y_5' + 0.001 \cdot X_5' = 0.30966701 + 0.001 \cdot 0.96670130 = 0.31063371130

 17.819309631011 + 0.057295760414 = 17.876605391425 < 18
 X_7' = X_6' - 0.001 \cdot Y_6' = 0.96639163299 - 0.001 \cdot 0.31063371130 = 0.96608099927870
 Y_7' = Y_6' + 0.001 \cdot X_6' = 0.31063371130 + 0.001 \cdot 0.96639163299 = 0.31160010293299

 17.876605391425 + 0.057295760414 = 17.933901151840 < 18
 X_8' = X_7' - 0.001 \cdot Y_7' = 0.96608099927870 - 0.001 \cdot 0.31160010293299 = 0.96576939917576701
 Y_8' = Y_7' + 0.001 \cdot X_7' = 0.31160010293299 + 0.001 \cdot 0.96608099927870 = 0.31256618393226870

 17.933901151840 + 0.057295760414 = 17.991196912254 < 18
 X_9' = X_8' - 0.001 \cdot Y_8' = 0.96576939917576701 - 0.001 \cdot 0.31256618393226870 = 0.96545683299183474130
 Y_9' = Y_8' + 0.001 \cdot X_8' = 0.31256618393226870 + 0.001 \cdot 0.96576939917576701 = 0.31353195333144446701

 17.991196912254 + 0.057295760414 > 18 => STOP

...

j = 17

 tan({\theta_i}) = 0.00000000000000001 ;  {\theta_i} = 0.000000000000000572957795130823
...

 17.999999999999999359893381975846 + 0.000000000000000572957795130823
 = 17999999999999999932851177106669 < 18

 X_{80}' = X_{79}' - 0.00000000000000001 \cdot Y_{79}'
       = 0.965408654718813480084209861507 - 0.00000000000000001 \cdot 0.313680286831855612515132822077
       = 0.965408654718813476947406993189

 Y_{80}' = Y_{79}' + 0.00000000000000001 \cdot X_{79}'
       = 0.313680286831855612515132822077 + 0.00000000000000001 \cdot 0.965408654718813480084209861507
       = 0.313680286831855622169219369265

On vérifie à l’aide de la commande bc (bc - An arbitrary precision calculator language) de la fondation GNU) :

# 25 décimales et une division
echo \’scale=  25;0.313680286831855622169219369265/0.965408654718813476947406993189' | bc
 0.3249196962329063248601754

18 décimales,  a(1) = arctg(1) = \pi /4, {\theta} = (18/45) \star \pi /4
 #  tg(\theta) = sin(\theta)/cos(\theta)
echo \’scale=  25;s(18/45 \star a(1))/c(18/45 \star a(1))' | bc -l
 0.3249196962329063261558713

 HP-15C: tan(18) = 0.324919696

Remarques

Nous avons vu et appliqué l’algorithme CORDIC.

  • L’algorithme CORDIC original (Volder) permet de dépasser l’angle initial, mais change de sens de rotation lorsqu’il dépasse l’angle initial.

Maintenant, on va voir comment les ingénieurs de HP ont fait :

  • ils ne dépassent jamais l’angle initial (comme l’exemple ci-dessus) ;
  • ils utilisent le radian à la place du degré ;
  • ils utilisent l’approximation suivante : {tan(x)} \approx {x} pour x petit.

Méthode HP

Remarque : tous les angles sont en radian, et

 18^\circ=0.314159265358979323 [rad].
 {\theta}={0.314159265358979323}=
+0 \star 0.78539816339744830
+3 \star  0.09966865249116202
+1 \star 0.00999966668666523
+5 \star 0.00099999966666686
+1 \star 0.00009999999966666
+5 \star 0.00000999999999966
+0 \star 00000364286582877

On vérifie que tan(0.00000364286582877)=0.00000364286582878611418, où 0.00000364286582877 représente le reste de l’addition. On utilisera ce reste comme valeur de départ pour {Y_0}.

{\theta} = 0.314159265358979323
{X_0} = 1
{Y_0} = 0.00000364286582877

j = 0 (0 boucle)



 tan({\theta_i}) = 1  {\theta_i} = 0.78539816339744830

j = 1 (3 boucles)


 tan({\theta_i}) = 0.1  {\theta_i} = 0.09966865249116202

 0 + 0.09966865249116202 = 0.09966865249116202 < 314159265358979323
 {X_1'} = {X_0 - 0.1 \star Y_0} = {1 - 0.1 \star 0.00000364286582877} = 0.99999963571341713
 {Y_1'} = {Y_0 + 0.1 \star X_0} = {0.00000364286582877 + 0.1 \star 1} = 0.10000364286582877

 X_2' = {X_1' - 0.1 \star Y_1'}
 Y_2' = {Y_1' + 0.1 \star X_1'}

 X_3' = {X_2' - 0.1 \star Y_2'}
 Y_3' = {Y_2' + 0.1 \star X_2'}

j = 2 (1 boucle)




 tan({\theta_i}) = 0.01  {\theta_i} = 0.00999966668666523

 X_4' = X_5' - 0.01 \star Y_5'
 Y_4' = Y_5' + 0.01 \star X_5'

j = 3 (5 boucles)



 tan({\theta_i}) = 0.001  {\theta_i} = 0.00099999966666686

 X_5' = X_6' - 0.001 \star Y_6'
 Y_5' = Y_6' + 0.001 \star X_6'
...

 X_9' = X_8' - 0.001 \star Y_8'
 Y_9' = Y_8' + 0.001 \star X_8'

j = 4 (1 boucle)




 tan({\theta_i}) = 0.0001  {\theta_i} = 0.00009999999966666

 X_{10}' = X_9' - 0.0001 \star Y_9'
 Y_{10}' = Y_9' + 0.0001 \star X_9'

j = 5 (5 boucles)




 tan({\theta_i}) = 0.00001  {\theta_i} = 0.00099999966666686

 X_{11}' = X_{10}' - 0.00001 \star Y_{10}'
 Y_{11}' = Y_{10}' + 0.00001 \star X_{10}'

...

 X_{15}' = X_{14}' - 0.00001 \star Y_{14}' = 0.96540865472374197
 Y_{15}' = Y_{14}' + 0.00001 \star X_{14}' = 0.31368028683345698

On vérifie :

# 18 décimales et une division
echo 'scale=18;0.31368028683345698/0.96540865472374197' | bc
0.324919696232906317
# 18 décimales, a(1) = arctg(1) = {\pi}/4, {\theta} = (18/45) \star {\pi}/4
# tg({\theta}) = sin({\theta})/cos({\theta})
echo 'scale=18;s(18/45 ^\circ a(1))/c(18/45 \star a(1))' | bc -l
0.324919696232906325

16 décimales exactes pour 15 itérations : REMARQUABLE !

 HP-15C: tan(18) = 0.324919696

Le nombre de 9 et de 6 apparaissant dans les constantes ci-dessus n’aura pas échappé à un oeil exercé. Raison pour laquelle, Dave Cochran concepteur de la HP-35 a choisi, afin d’économiser de la place mémoire, de ne mettre qu’une constante en mémoire et de calculer les autres à l’aide d’une boucle et de chiffres donnant la position et le nombre de 9 (ou de 6).

Fonctions trigonométriques

À l’aide du calcul de tan({\theta}), on déduit les fonctions trigonométriques sin & cos à l’aide des équations suivantes :

sin({\theta}) = {tan({\theta})\over {\sqrt{1 + {tan({\theta})^2}}}
cos({\theta}) = {cot({\theta})\over {\sqrt{1 + {cot({\theta})^2}}}

 tan({\theta}) = {Y \over X} et  cot({\theta}) = {X \over Y}
On inverse simplement X et Y, puis on utilise la même équation et donc le même algorithme de calcul.

Il reste encore à savoir calculer une racine carrée, ce qui est justement le sujet du premier article de William E. Egbert.

Bonus

On peut aussi utiliser les nombres complexes pour effectuer des rotations, on a vu que l’on a tourné :

  • 0 fois pour j = 0
  • 3 fois pour j = 1
  • 1 fois pour j = 2
  • 5 fois pour j = 3
  • 1 fois pour j = 4
  • 5 fois pour j = 5

Ce qui nous donne :
 X + i \cdot Y \approx (1 + i \cdot0)  \cdot [{(1 + i \cdot 0.1)^3  \cdot (1 + i \cdot 0.01)  \cdot (1 + i \cdot0.001)^5  \cdot (1 + i \cdot 0.0001)  \cdot (1 + i \cdot 0.0001)^5} ]
 X + i \cdot Y  \approx 0.965409797 + i \cdot 0.313676769
 tan(18) = 0.324919  \approx {Y \over X} = {0.313676769 \over 0.965409797} = 0.324915

On retrouve facilement l’algorithme ci-dessus !

[1] HPJ-1977-05.pdf

[2] HPJ-1977-06.pdf

[3] HPJ-1977-11.pdf

[4] HPJ-1978-04.pdf

[5] VOLDER, Jack E., The CORDIC Trigonometric Computing Technique. IRE Transactions on Electronic Computers, septembre 1959

[6] ROULET François, La saga HP-35. FI7/07



Cherchez ...

- dans tous les Flash informatique
(entre 1986 et 2001: seulement sur les titres et auteurs)
- par mot-clé

Avertissement

Cette page est un article d'une publication de l'EPFL.
Le contenu et certains liens ne sont peut-être plus d'actualité.

Responsabilité

Les articles n'engagent que leurs auteurs, sauf ceux qui concernent de façon évidente des prestations officielles (sous la responsabilité du DIT ou d'autres entités). Toute reproduction, même partielle, n'est autorisée qu'avec l'accord de la rédaction et des auteurs.


Archives sur clé USB

Le Flash informatique ne paraîtra plus. Le dernier numéro est daté de décembre 2013.

Taguage des articles

Depuis 2010, pour aider le lecteur, les articles sont taggués:
  •   tout public
    que vous soyiez utilisateur occasionnel du PC familial, ou bien simplement propriétaire d'un iPhone, lisez l'article marqué tout public, vous y apprendrez plein de choses qui vous permettront de mieux appréhender ces technologies qui envahissent votre quotidien
  •   public averti
    l'article parle de concepts techniques, mais à la portée de toute personne intéressée par les dessous des nouvelles technologies
  •   expert
    le sujet abordé n'intéresse que peu de lecteurs, mais ceux-là seront ravis d'approfondir un thème, d'en savoir plus sur un nouveau langage.