Algorithme cordic

Rappels de trigonométrie.

Le cercle de rayon \(1\) centré en l'origine \(O:=(0,0)\) du plan euclidien est appelé cercle tri­go­no­mé­tri­que. Soit \({\color{yellow}M}\) un point de ce cercle. On note respectivement \(\color{yellow}{C}\) et \(\color{yellow}{S}\) les projections orthogonales de \(M\) sur l'axe des abscisses et sur l'axe des ordonnées et \(\theta\) l'angle \(\widehat{COM}\) qui est défini en première ap­pro­che comme la mesure de l'arc* entre le point \(I:=(1,0)\) et le point \(M:=(x,y)\)*La mesure du cercle en­tier est donc \(2\pi\) radians. Par­fois nous uti­li­se­rons le degré comme unité de mesure des angles avec \(2\pi\equiv 360°\)..

L'abscisse \(x\) et l'ordonnée \(y\) du point \(M\) définissent respectivement le cosinus et le sinus de l'angle \(\theta\), le ratio \(y/x\) définit sa tangente : \begin{equation}\label{eq:deftan} \cos \theta:=x=\overline{OC}\qquad \sin \theta:=y=\overline{OS} \qquad \tan \theta:=\frac{y}{x}=\frac{\sin \theta}{\cos \theta} \end{equation}

Par construction, le triangle \(OCM\) est rectangle en \(C\) et le théorème de Pythagore fournit un premier résultat : \begin{equation}\label{eq:pythagore} \def\R{{\mathbb R}} \def\Z{{\mathbb Z}} \cos^2\theta+\sin^2\theta=1. \end{equation}

Le domaine de définition des fonctions cosinus et sinus est étendu à l'ensemble des réels \(\R\), l'angle étant calculé modulo \(2\pi\) et \(\cos\theta\in[-1,1]\), \(\sin\theta\in[-1,1]\). La fonction cosinus s'annule en tout angle \(\theta=k\pi+\frac{\pi}{2}\) où \(k\in\Z\), le domaine de définition de la fonction tangente est donc \(\R\setminus\frac{\pi}{2}+\pi\Z\) et elle est à valeur dans \(\R\). Les fonctions réciproques sont notées \(\arccos\), \(\arcsin\) et \(\arctan\).

On ne représente ci-dessous que la part du cercle trigonométrique dans le premier quadrant et l'angle \(\theta\) est mesuré en degrés pour plus de simplicité :

Cercle trigonométrique dans le premier quadrant.

Avec un angle \(\color{yellow}{\theta}\) de degrés, le point \(\color{yellow}{M}\) a pour coor­don­nées \(x=\cos\theta=\) et \(y=\sin\theta=\) . La tangente de l'angle est donc égale à . Notons que la valeur de la tangente peut éga­le­ment se lire sur le graphique. En effet, en notant \(T\) l'intersection (quand elle existe) de la droite \((OM)\) et de la droite d'équation \(x=1\), les triangles rectangles \(OCM\) et \(OIT\) sont ho­mo­thé­ti­ques donc \(\tan\theta=|IT|\) d'après le théorème de Thalès.

Si \(ABC\) désigne un triangle rectangle en \(C\) quelconque, le théorème de Thalès permet également d'affirmer que \(\cos\widehat{ACB}=\frac{AC}{AB}\) (côté adjacent sur hypothénuse) et que \(\sin\widehat{ACB}=\frac{BC}{AB}\) (côté opposé sur hypothénuse), que l'hy­po­thé­nu­se \(AB\) soit ou non égale à \(1\).
Soit \(\theta\in\R\). À l'aide des symétries sur le cercle, calculez \(f(-\theta)\), \(f(\theta\pm\pi)\), \(f(\pi-\theta)\), \(f(\theta+\frac{\pi}{2})\), \(f(\theta-\frac{\pi}{2})\) et \(f(\frac{\pi}{2}-\theta)\) où \(f\) désigne les trois fonctions \(\cos\), \(\sin\) et \(\tan\). Exprimez \(\sin a\cos b\) sous forme de somme de sinus.

On rappelle quelques formules utiles : \begin{align*} \cos(a+b)&=\cos a\cos b-\sin a\sin b & \cos a\cos b&=\frac{1}{2}\left(\cos(a-b)+\cos(a+b)\right)\\ \sin(a+b)&=\sin a\cos b+\cos a\sin b &\sin a\sin b&=\frac{1}{2}\left(\cos(a-b)-\cos(a+b)\right) \\ \cos a+\cos b&=2\cos\frac{a+b}{2}\cos\frac{a-b}{2} & \sin a+\sin b&=2\sin\frac{a+b}{2}\cos\frac{a-b}{2}\\ \end{align*}

Démontrez que la surface d'un triangle \(ABC\) quelconque est égale à \(\frac{1}{2}|AB|\,|AC| \sin\widehat{BAC}\)

Étude du problème. Rotations.

Pour un angle donné \(\theta\), on cherche à calculer efficacement \(\sin\theta\), \(\cos\theta\) et \(\tan\theta\). Une première idée est d'utiliser un développement limité de ces fonctions à l'angle \(\theta\) considéré, mais cette approche gé­né­ri­que n'exploite pas les spécificités des fonctions trigonométriques.

L'algorithme CORDIC — ac­ro­ny­me de COordinate Rotation DIgital Computer — que nous pré­sen­te­rons à la section suivante est plus efficace. Il a été développé par Jack Volder en 1956 et consiste à approcher le point \(M=(\cos\theta,\sin\theta)\) par un point \(N\) du cercle unité ini­tialement placé en \(I=(1,0)\) et qui va subir une suc­cession de ro­ta­tions d'angles de plus en plus petits. À la manière de la di­cho­to­mie, le sens de la rotation dépend de la position du point \(N\) par rapport à \(M\) sur le cercle. Elle s'effectue dans le sens tri­gono­métrique si \(N\) est en-dessous de \(M\), dans le sens horaire dans le cas contraire. La première rotation se fait donc dans le sens trigonométrique puisque \(N\) commence son périple du point \(I=(1,0).\)

On désigne par \(\gamma_i\) l'angle de la \((i+1)\)-ème rotation, l'indexation débutant à 0 pour ne pas alourdir inu­ti­le­ment les futures expressions. Ces angles étant supposés positifs, le sens de la rotation \(i\) est noté \(s_i=\pm 1\), il est positif dans le sens trigonométrique, négatif sinon. On note \((N_n)_{n\in{\mathbf N}}\) la suite des points obtenus après**Le point \(N_0\) n'est donc pas le point \(I\). la \((n+1)\)-ème rotation (d'indice \(n\)) et \(\theta_n\) l'angle formé par le point \(N_n\). Par construction \begin{equation} \forall n\in{\mathbf N},\quad \theta_n=\sum_{i=0}^{n}s_i\gamma_i. \end{equation}

Après \(n+1\) rotations, les coordonnées \(N_n=(\cos\theta_n,\sin\theta_n)\) du point \(N\) fournissent donc une ap­pro­xi­ma­tion de \(\cos\theta\) et \(\sin\theta\). Une rotation s'exprime formellement à l'aide d'une application linéaire du corps des complexes \(\mathbf{C}\) dans lui même :

On appelle rotation d'angle \(\gamma\), l'application linéaire \(\text{rot}_\gamma:\mathbf{C}\rightarrow\mathbf{C}\) définie par \(u\mapsto e^{i\gamma}u\). En identifiant \(\mathbf{C}\) à \(\mathbf{R}\times\mathbf{R}\) et en notant \(u=(x,y)\), ceci se traduit matriciellement par \begin{equation}\label{eq:rotation} (x,y)\overset{\text{rot}_\gamma}{\longmapsto} \begin{pmatrix} \cos \gamma & -\sin\gamma\\ \sin \gamma & \cos\gamma \end{pmatrix} \begin{pmatrix} x\\ y\\ \end{pmatrix}. \end{equation}

Nous supposerons que \(\theta\in[0,\frac{\pi}{2}]\) puisque l'on peut déduire les cosinus et sinus pour les autres valeurs de \(\theta\) à l'aide des symétries usuel­les, cf. exercice suivant.

En supposant que l'on dispose d'un algorithme calculant le cosinus et le sinus pour un angle \(\theta\in[0,\pi/2]\), expliciter toutes les symétries nécessaires au calcul du cosinus et du sinus d'un angle \(\theta\) quelconque.

On note \(R_i\) la matrice de rotation d'angle \(s_i\gamma_i\). Nous allons réécrire l'expression \((\ref{eq:rotation})\) afin d'optimiser les calculs. On factorise par \(\cos\gamma_i\) :

\begin{equation}\label{eq:cosfact} R_i=\cos\gamma_i\begin{pmatrix} 1& -s_i\tan\gamma_i\\ s_i\tan\gamma_i & 1 \end{pmatrix}. \end{equation}

Pour que ce processus d'approximation de l'angle \(\theta\) soit valide avec cette succession de rotations d'angles décroissants prédéfinis \((\gamma_n)_{n\in{\mathbf N}}\), il faut qu'à chaque étape l'angle restant à combler entre \(N\) et \(M\) soit inférieur à l'angle de la rotation qui vient d'être réalisée, il est donc nécessaire que les iné­ga­li­tés suivantes soient satisfaites :

\begin{equation}\label{eq:approxtheta} \forall n\in{\mathbf N},\quad\Big|\theta-\theta_n\Big|\leq\gamma_{n}. \end{equation}

La suite \((\gamma_n)_{n\in{\mathbf N}}\) doit évidemment converger vers \(0\) puisque l'objectif est d'approcher \(\theta\). Cependant la convergence de la série \(\theta_n\) n'est pas nécessairement possible pour n'importe quelle valeur de \(\theta\). Les angles \(\gamma_i\) étant tous positifs, il est clair que les valeurs limites de la série \(\theta_n\) sont obtenues avec des signes \(s_i\) tous égaux à \(-1\) ou \(+1\), ce qui impose pour \(\theta\) que  :

\begin{equation}\label{ineq:encadrementtheta} -\lim_{n\rightarrow+\infty}\sum_{i=0}^{n}{\gamma_i}\leq\theta\leq \lim_{n\rightarrow+\infty}\sum_{i=0}^{n}{\gamma_i}. \end{equation}

Revenons à l'inégalité \((\ref{eq:approxtheta})\). Supposons qu'elle soit satisfaite au rang \(n\), observons ce que celà entraîne pour le rang suivant. Nous avons vu que le sens \(s_i\) de la rotation \(R_i\) est déterminé par la position du point \(N\) par rapport au point \(M\) sur le cercle unité. Il faut donc distinguer le cas \(\theta-\theta_n\geq 0\) du cas \(\theta-\theta_n < 0\). Nous n'étudierons que le premier cas, l'autre est laissé en exercice.

Si \(\theta-\theta_n\geq 0\), alors la rotation suivante se fait dans le sens trigonométrique donc \(s_{n+1}=1\). D'après \((\ref{eq:approxtheta})\), on a

\begin{equation}\label{ineq:Sn} 0\leq \theta-\underbrace{\sum_{i=0}^{n}s_i\gamma_i}_{\theta_n}\leq\gamma_n \end{equation}

et donc en soustrayant l'angle \(\gamma_{n+1}\) :

\begin{align*} -\gamma_{n+1}&\leq \theta-\sum_{i=0}^{n}s_i\gamma_i-s_{n+1}\gamma_{n+1}\leq\gamma_n-\gamma_{n+1}\\ -\gamma_{n+1}&\leq \theta-\sum_{i=0}^{n+1}s_i\gamma_i\leq\gamma_n-\gamma_{n+1}. \end{align*}

On en déduit que

\begin{equation} \Big|\theta-\theta_{n+1}\Big|\leq\max\{\gamma_{n+1},\gamma_{n}-\gamma_{n+1}\}. \end{equation}

Ainsi pour que l'encadrement \((\ref{ineq:Sn})\) soit satisfait au rang \(n+1\), c'est-à-dire \(0\leq\theta-\theta_{n+1}\leq\gamma_{n+1}\), il faut que \(\gamma_n-\gamma_{n+1}\leq\gamma_{n+1}\), ce qui nous donne finalement la condition nécessaire et suffisante sui­van­te :

\begin{equation}\label{ineq:exo} \forall n\in{\mathbf N},\quad \frac{1}{2}\gamma_n\leq\gamma_{n+1}\leq\gamma_n. \end{equation}

Autrement dit, la décroissance des angles pour l'approximation ne peut pas être trop rapide, chaque nouveau pas doit être inférieur au pas précédent mais doit également rester supérieur ou égal à la moitié de ce pas. Pour achever la preuve par récurrence de \((\ref{ineq:Sn})\), il faut que la condition initiale \(|\theta-s_0\gamma_0|\leq\gamma_0\) soit satisfaite.

Démontrez les inégalités \((\ref{ineq:exo})\) dans le cas où \(\theta-{\theta_n} < 0\).

L'algorithme cordic

Nous connaissons à présent les contraintes sur l'angle \(\theta\) à approximer et sur la suite des angles pré­dé­fi­nis \((\gamma_n)_{n\in{\mathbf N}}\) pour ce faire, mais quelle suite utiliser ?

La rotation \(R_i\) définie en \((\ref{eq:cosfact})\) engendre une multiplication de chaque coordonnée du point \(N\) par la valeur \(\tan\gamma_i\) avec éventuellement un changement de signe. Il nous faut à la fois contourner le calcul de la tangente (sans quoi nous serions revenus au point de départ, l'objectif initial étant précisément de calculer efficacement les fonctions tri­go­no­mé­tri­ques) et optimiser ces deux produits. La mul­ti­pli­ca­tion (resp. la division) par \(2\) se traduit en machine par un simple décalage des bits de la re­pré­sen­ta­tion binaire des nombres à gauche (resp. à droite) et une telle opération est réalisée en un seul cy­cle machine**c'est également le cas pour un décalage cir­cu­lai­re d'un nombre ar­bi­trai­re de positions avec un circuit spécialisé ap­pe­lé barrel shifter. On va donc fixer \(\tan\gamma_i=2^{-i}\), soit \(\gamma_i=\arctan 2^{-i}\) :

\begin{equation}\label{eq:cosarctan} R_i=\cos(\arctan 2^{-i})\begin{pmatrix} 1& -s_i2^{-i}\\ s_i2^{-i} & 1 \end{pmatrix}. \end{equation}

Avant de poursuivre la simplification des calculs, vérifions tout d'abord que la somme de la série de terme général \(\gamma_i:=\arctan 2^{-i}\) existe.

On note \(T'\) le point d'intersection entre la tangente au cercle unité en \(M\) et l'axe des abscisses. Démontrez que \(|MT'|=|IT|=\tan\theta\). En déduire que \(\forall\theta\in[0,\frac{\pi}{2}],\ \sin\theta < \theta < \tan\theta\).

D'après l'exercice précédent, on a \(\arctan 2^{-i}<2^{-i}\) et donc \begin{equation*} \sum_{i=0}^{n}\arctan 2^{-i}<\sum_{i=0}^{n}2^{-i} =\frac{1}{2^n}\sum_{i=0}^{n}2^i=\frac{1}{2^n}2^{n+1}=2. \end{equation*}

La série étant croissante et majorée, elle converge. Expérimentalement, on obtient

\begin{equation}\label{eq:bound} {\color{#0AF}\theta_\infty}:=\lim_{n\rightarrow+\infty}\sum_{i=0}^n\arctan 2^{-i}\simeq1.743286620472341 \end{equation}

il faut donc que \(\theta\in[-{\color{#0AF}\theta_\infty},{\color{#0AF}\theta_\infty}]\) ce qui est satisfait puisque nous avons supposé que \(\theta\in[0,\frac{\pi}{2}]\). Il faut encore s'assurer que la suite \((\gamma_n)_{n\in\mathbf N}\) satisfait \((\ref{ineq:exo})\). La fonction \(\arctan\) étant stric­te­ment croissante sur \([0,\frac{\pi}{2}]\), il ne reste qu'à vérifier que pour tout en­tier \(n\) on a \(\frac{1}{2}\gamma_n\leq\gamma_{n+1}\) autrement dit que \[\forall n\in{\mathbf N},\quad \arctan 2^{-n}\leq2\arctan 2^{-n-1}.\] La preuve est laissée en exercice.

Démontrez que la fonction \[2\arctan\frac{x}{2}-\arctan x\] est toujours strictement positive dans l'intervalle \([0,+\infty[\). Indication : on rappelle que la dérivée de la fonction \(\arctan x\) est la fonction \[\frac{1}{1+x^2}\] et que la dérivée d'une fonction composée \(g\circ f\) est la fonction \(f'.(g'\circ f)\).

Reprenons l'étude des rotations à partir de \((\ref{eq:cosarctan})\). La matrice de la composition de deux ro­ta­tions étant égale au produit des matrices de ces rotations, la matrice \(\Theta_n\) obtenue après \(n+1\) rotations est donnée par l'équation

\begin{equation}\label{eq:avant} \Theta_{n}:=\prod_{i=0}^nR_i=\prod_{i=0}^n \cos(\arctan 2^{-i})\begin{pmatrix} 1& -s_i2^{-i}\\ s_i2^{-i} & 1 \end{pmatrix}. \end{equation}
En utilisant les identités \((\ref{eq:deftan})\) et \((\ref{eq:pythagore})\), démontrez que \begin{equation}\label{eq:arctan} \forall \theta\in\R,\quad \cos(\arctan \theta)= \frac{1}{\sqrt{1+\theta^2}} \end{equation}

En notant \(K_n\) le produit

\begin{equation} K_n:=\prod_{i=0}^n\frac{1}{\sqrt{1+2^{−2i}}} \end{equation}

la matrice devient

\begin{equation} \Theta_{n}=K_n\prod_{i=0}^n\begin{pmatrix} 1& -s_i2^{-i}\\ s_i2^{-i} & 1 \end{pmatrix}. \end{equation}
Vérifiez que la première rotation \(R_0\) est une rotation d'angle \(\gamma_0=\frac{\pi}{4}\).

Pour \(n=\) , on obtient \(K_n\approx\) . Expérimentalement et à \(10^{-16}\) près, la limite

\[K:=\displaystyle\lim_{n\rightarrow+\infty}K_n\]

est égale à \(0.6072529350088814\), valeur atteinte dès que \(n\geq 27\). On constate qu'après seulement \(10\) itérations, \(K_n\) est égale à \(K\) à \(10^{-5}\) près. On peut donc légitimement remplacer \(K_n\) par la constante ci-dessus si l'on suppose que le nombre d'itérations \(n\geq 27\). Bien entendu, si l'on souhaite atteindre une précision supérieure à 16 chiffres après la virgule en base 10, il faut reéxaminer la situation.

Le point de départ étant le point \(I=(1,0)\), l'approximation des fonctions cosinus et sinus après \(n+1\) rotations sont fournies par

\begin{equation} \color{#EEF} \label{eq:finale} (\cos\theta,\sin\theta)\approx \prod_{i=0}^n\begin{pmatrix} 1& -s_i2^{-i}\\ s_i2^{-i} & 1 \end{pmatrix} \begin{pmatrix} K\\ 0 \end{pmatrix} \end{equation}

Maintenant que la méthode a été validée, il ne reste plus qu'à écrire l'algorithme. Le sens de la rotation courante étant fixé par la position de \(N\) par rapport à \(M\), c'est-à-dire la différence entre l'angle \(\theta\) et son approximation, il faut précalculer les valeurs \(\gamma_i=\arctan 2^{-i}\) nécessaires selon la précision souhaitée, ici nous avons supposé \(n=27\). 

Cordic(theta)
données
   theta: réel
constantes
   K ← 0.6072529350088814   
   ATAN ← table des arctan 2**(-i)
   n ← 27
variables
   i, s: entier
   x, y, aux, ecart: réels
DEBUT
   x ← K
   y ← 0
   ecart ← theta
   i ← 0
   TQ (i ≤ n) FAIRE
      s ← signe(ecart)
      aux ←  x - (s * (y >> k))
      y ← y + (s * (x >> k))
      x ← aux
      ecart ← ecart - (s * ATAN[i])
   FTQ
   RENVOYER (x,y)
FIN

L'animation ci-dessous montre l'évolution du point \(\color{#00AAFF}N\) obtenu avec cet algorithme pour un angle à approximer \(\color{yellow}{\theta}\) de . Après la rotation R, on obtient l'estimation suivante :

\(\cos\theta\approx\)
\(\sin\theta\approx\)

L'erreur commise en remplaçant chaque coefficient \(K_n\) par la limite \(K\) est bien visible pour les premières occurrences, mais s'estompe rapidement.

Estimation des fonctions sinus et cosinus avec l'al­go­ri­thme CORDIC.

Travaux pratiques

Écrivez l'algorithme CORDIC en langage \(C\). Comparez les résultats obtenus avec ceux de la bibliothèque mathématique fournie.
On rappelle que le développement en série entière des fonctions cosinus et sinus sont données par \begin{align*} \forall \theta\in\R,\quad \cos\theta&=\sum_{n=0}^{+\infty}(-1)^n\frac{x^{2n}}{(2n)!} \\ \forall \theta\in\R,\quad \sin\theta&=\sum_{n=0}^{+\infty}(-1)^n\frac{x^{2n+1}}{(2n+1)!} \\ \end{align*} Écrivez un algorithme qui utilise ces développements pour estimer les fonctions cosinus et sinus. Comparer les temps de calcul pour estimer ces fonctions entre cet algorithme et l'algorithme CORDIC pour la même précision.
Réécrivez l'algorithme CORDIC en représentant l'angle \(\theta\in[0,\frac{\pi}{2}]\) et les coordonnées \(x\) et \(y\) des points du cercle unité avec un entier codé sur \(64\) bits correspondant au développement binaire de \(\theta.\) On supposera que les bits de l'entier en question sont numérotés du poids faible au poids fort. Si on note \(T\) cet entier et \(b_i\) le \(i\)-ème bit de \(T\), alors \[\theta =\sum_{i=0}^{63}\frac{1}{2^i}.\] Donnez un encadrement des angles \(\theta\) que l'on peut ainsi représenter.