Que calcule la calculette ?

Vous êtes-vous demandé comment fait votre calculette pour retourner instantanément avec 10 chiffres de précision n'importe quelle valeur d'une fonction usuelle ? À la lecture de ce qui précède, vous êtes peut-être tenté de croire à une approximation polynomiale, calculée par la méthode des différences. Non : votre calculette utilise probablement l'algorithme CORDIC (pour Coordinate Rotation Digital Computer), inventé par J. Volder en 1959, et généralisé par J. Walther en 1971⁵.

Voici le principe, pour les fonctions trigonométriques. Soit à calculer $\sin(\theta)$ et $\cos(\theta)$ , pour un certain angle $\theta$ . L'idée consiste à approcher $\theta$ par la suite $(\theta_n)_{n\in\mathbb{N}}$ , définie par $\theta_0=0$ et pour tout $n\geqslant 0$  :

$$\theta_{n+1} = \left\{ \begin{array}{lcl} \theta_{n}&\mbox{s... ...\arctan(2^{-n})&\mbox{si}&\theta<\theta_{n} \end{array}\right.$$

La série de terme général $\arctan(2^{-n})$ est convergente. De plus, pour tout $n\in\mathbb{N}$  :

$$\arctan(2^{-n})< \sum_{i=n+1}^{+\infty} \arctan(2^{-i})\;.$$

Démontrez-le, et déduisez-en que la suite $(\theta_n)$ converge vers $\theta$ , pour tout $\theta$ compris entre deux valeurs extrêmes :

$$-\sum_{n=0}^{+\infty} \arctan(2^{-n})\leqslant \theta \leqslant +\sum_{n=0}^{+\infty} \arctan(2^{-n})\;.$$

Ce n'est pas une limitation gênante : comme vous le savez, on peut toujours se ramener à un angle compris entre 0 et $\pi/2$ , quitte à changer ensuite le signe du résultat.

n=[1:100]; sum(atan(2^(-n)))/%pi

Les angles $\arctan(2^{-n})$ sont fixes : ils peuvent être tabulés et mémorisés. Peu de valeurs sont nécessaires, car pour $n$ assez grand, $2^{-n}$ est une bonne approximation de $\arctan(2^{-n})$ .

d=2^(-[0:10]); atan(d)-d

Notons $\alpha_n$ l'angle $\theta_{n+1}-\theta_n=\pm\arctan(2^{-n})$ . Les valeurs approchant le cosinus et le sinus de $\theta$ sont $x_n=\cos(\theta_n)$ , $y_n=\sin(\theta_n)$ . Le vecteur de coordonnées $(x_{n+1},y_{n+1})$ se déduit du précédent par une rotation d'angle $\alpha_{n}$ (d'où le nom de l'algorithme) :

$$\left(\begin{array}{r}x_{n+1}\\ y_{n+1}\end{array}\right) = \left... ...a_n)&1\end{array}\right) \left(\begin{array}{r}x_{n}\\ y_{n}\end{array}\right)$$

Pour tout $n\in\mathbb{N}$ , posons

$$K_n=\prod_{i=0}^{n-1} \cos(\alpha_i)=\prod_{i=0}^{n-1} \cos(\arctan(2^{-i})) =\prod_{i=0}^{n-1} \frac{1}{\sqrt{1+2^{-2i}}}\;.$$

Posons $(X_0,Y_0)=(x_0,y_0)=(1,0)$ , et pour tout $n\geqslant 1$ , $(X_n,Y_n)=K_{n}^{-1}(x_n,y_n)$ . Il vient :

$$\left(\begin{array}{r}X_{n+1}\\ Y_{n+1}\end{array}\right) = \left... ...a_n)&1\end{array}\right) \left(\begin{array}{r}X_{n}\\ Y_{n}\end{array}\right)$$

La puissance de l'algorithme réside dans cette formule de récurrence : elle ne comporte que des additions, et des multiplications par $\tan(\alpha_n)$ . Or par définition de $\alpha_n$ , $\tan(\alpha_n)=\pm 2^{-n}$ . Multiplier par $\tan(\alpha_n)$ en arithmétique binaire, c'est décaler la virgule de $n$ places, et éventuellement inverser le bit de signe. Le calcul de $(X_n,Y_n)$ est donc peu coûteux. Pour revenir à $(x_n,y_n)$ , il faut multiplier par $K_n$  : les valeurs de $K_n$ sont fixes et peuvent être mémorisées. Peu de valeurs sont nécessaires car le produit converge vite.

cumprod(cos(atan(2^(-[0:15]))))

Voici en Scilab le calcul de $\sin(1)$ et $\cos(1)$ (non optimisé bien sûr).

theta=1; tn=0; v=[1;0]; V=[]; n=40;
for i=0:n, 
    s=sign(theta-tn); tn=tn+s*atan(2^(-i)); 
    b=s*2^(-i); M=[1,-b;b,1]; v=M*v; V=[V,v]; 
end;
K=cumprod(cos(atan(2^(-[0:n])))); sc=V.*[K;K]
sc-[cos(theta);sin(theta)]*ones(1,n+1)

Pour les autres fonctions usuelles, il faut définir une autre notion de «rotation», mais le principe reste le même : une récurrence linéaire peu coûteuse, éventuellement suivie d'une multiplication par des valeurs mémorisées.

Notes

... 1971⁵

J. Walther : A unified algorithm for elementary functions, Joint Computer Conference Proceedings, 38 p. 379-385 (1971)