- CORDIC est un algorithme qui transforme des calculs complexes en opérations principalement fondées sur des additions et des décalages de bits, afin de calculer des fonctions trigonométriques comme
sin, cos et tan sans FPU ni grande table de correspondance
- Cette approche est utile dans les environnements embarqués, en particulier sur les microcontrôleurs peu performants et les FPGA, plutôt que dans les systèmes haute performance ; il est difficile d’en juger la valeur à la seule vitesse
- En utilisant la virgule fixe plutôt que la virgule flottante, on peut diviser un
int32_t en 16 bits de poids fort pour la partie entière et 16 bits de poids faible pour la partie fractionnaire, afin de représenter environ de -32768.99997 à 32767.99997
- En faisant tourner un vecteur vers un angle cible par angles de plus en plus petits, avec une table de 16 valeurs
atan(2**-i) et une valeur initiale x=39796, on peut remplacer les multiplications à chaque itération par des décalages de bits
- Pour l’angle d’exemple
0.9152, après 16 itérations, l’erreur absolue de sin(0.9152) descend à environ 0.00000956, et celle de cos(0.9152) à environ 0.0000434
Les environnements de calcul où CORDIC est adapté
- CORDIC est un algorithme destiné à calculer des fonctions trigonométriques comme
sin, cos et tan sur du matériel basse consommation
- Il fonctionne même dans des environnements dépourvus de FPU, c’est-à-dire d’unité de calcul en virgule flottante, ou lorsqu’il est difficile d’utiliser de grandes tables de correspondance
- Les calculs effectifs reposent surtout sur de simples additions et décalages de bits
- Il combine des idées de calcul vectoriel, de trigonométrie, de convergence et d’informatique pour approximer des fonctions complexes avec des opérations simples
- Sur du matériel haute performance, cette technique n’est pas forcément nécessaire
- Sa principale cible est l’embarqué
- Elle convient particulièrement aux microcontrôleurs peu performants et aux FPGA
- Il peut exister du matériel ou des périphériques plus rapides, mais la vitesse n’est pas le seul critère d’utilité
Représentation en virgule fixe pour éviter la virgule flottante
- Même les fonctions qui renvoient des valeurs entre
-1.0 et 1.0, comme sin(x), n’ont pas nécessairement besoin d’être représentées en virgule flottante
- La virgule fixe représente des nombres rationnels à l’intérieur d’un type entier en fixant la position de la virgule
- L’exemple divise un
int32_t en 16 bits de poids fort pour la partie entière et 16 bits de poids faible pour la partie fractionnaire
- Dans ce cas, la plage va approximativement de
-32768.99997 à 32767.99997
- Selon l’endroit où l’on place la virgule, on peut arbitrer entre l’étendue de la partie entière et la précision de la partie fractionnaire
- La valeur elle-même reste un
int32_t, auquel le programmeur attribue une signification supplémentaire dans l’agencement des bits
Conversion en virgule fixe et opérations de base
- Si la précision fractionnaire est de 16 bits, on peut convertir une valeur float comme
42.01 en virgule fixe en la multipliant par (1 << 16)
42.01 * (1 << 16), converti en int32_t, donne 2753167- Pour revenir en float, on calcule
2753167 / (1 << 16), ce qui donne environ 42.0099945
- On peut aussi encoder directement des valeurs comme
1.5 sans utiliser du tout la virgule flottante
- La partie entière
1 est portée par (1 << 16)
- La moitié fractionnaire peut être placée à
0x7fff, la valeur médiane entre 0x0000 et 0xffff
- Le résultat de cette méthode est
98303 en décimal
- Entre des valeurs utilisant le même facteur d’échelle, l’addition et la soustraction fonctionnent telles quelles
- Pour la multiplication, on multiplie deux valeurs en virgule fixe, puis on décale de nouveau le résultat vers la droite du facteur d’échelle
- Pour la division, on peut obtenir une précision supplémentaire en décalant d’abord le dividende vers la gauche du facteur d’échelle, puis en le divisant par le diviseur
Approximer les fonctions trigonométriques par rotation de vecteur
- CORDIC est l’abréviation de « co-ordinate rotation digital computer » et a été créé au milieu des années 1950
- L’idée centrale consiste à faire tourner un vecteur sur le cercle unité par angles de plus en plus petits, afin que, lorsqu’il atteint l’angle cible, ses composantes deviennent les valeurs du sinus et du cosinus
- Le processus ressemble à une recherche binaire
- On se déplace vers l’angle cible avec un grand angle
- On vérifie si l’on a dépassé la cible
- Puis on répète des rotations dans le sens horaire ou antihoraire avec des angles plus petits
- Pour calculer
sin(0.7), par exemple, on part du vecteur initial (1, 0) et de la cible 0.7 radian
- On commence par une rotation antihoraire de
0.7853 radian, soit 45˚
- La cible restante devient
0.7 - 0.7853 = -0.0853
- Comme la valeur est négative, on effectue ensuite une rotation horaire de
0.3926 radian, soit 22.5˚
- Ensuite, selon le signe de la cible restante, on change de direction avec des angles plus petits, comme
0.1963 radian
- Après 16 itérations, le vecteur est presque aligné sur l’angle cible initial,
y donne une approximation de sin(a) et x une approximation de cos(a)
Réduire les opérations coûteuses dans la matrice de rotation
- Une rotation vectorielle classique utilise une multiplication matricielle contenant des sinus et des cosinus
- CORDIC utilise des identités trigonométriques pour reformuler la matrice de rotation autour de
tan(a)
- Au départ, comme on utilise des angles de rotation fixes tels que
45˚, 22.5˚ et 11.25˚, les valeurs de tan(a) peuvent être placées dans une table précalculée
- Cette table ne nécessite que 16
uint32_t, soit 64 octets
- À titre de comparaison, une table
sin(x) non optimisée contenant 4096 valeurs de -1 à 1 nécessiterait 16 KiB, avec une précision considérée comme faible
- Le terme
cos(a) qui précède chaque rotation apparaît à chaque itération, mais le produit de tous ces termes converge vers une constante
- Avec des angles du type
45˚, 22.5˚, 11.25˚, ce produit vaut environ 0.6366
- Il suffit de multiplier par cette constante une seule fois après toutes les itérations
Choisir les angles pour ne garder que décalages et additions
- Pour éliminer les multiplications, on choisit les angles de sorte que le résultat de
tan(a) soit toujours une puissance inverse de 2
- Pour cela, on crée une table de 16 entrées contenant, pour chaque itération de
i=0 à 15, la valeur atan(2**-i)
- Les angles de rotation effectifs deviennent
45˚, 26.565˚, 14.036˚, 7.125˚, etc.
- Les angles ne diminuent pas exactement de moitié, mais le processus converge tout de même vers le bon résultat avec ces angles
- La multiplication par
tan(a) est remplacée par un décalage de bits correspondant au numéro d’itération i
- Le produit des termes
cos(a) est lui aussi recalculé en fonction de ce nouveau choix d’angles
- Sa valeur est d’environ
0.60725
- En virgule fixe sur 16 bits, cela donne
39796
- Au lieu de multiplier par cette valeur à la fin, on peut initialiser le
x du vecteur à 39796 plutôt qu’à 1
Procédure de l’algorithme
- Lors de l’étape de précalcul, on crée une table dont chaque entrée vaut
atan(2**-i), puis on convertit chaque valeur en virgule fixe
- La formule de conversion est
atan(2**-i) * (1 << 16)
- Pour calculer
sin ou cos, on convertit aussi l’angle d’entrée en virgule fixe
- L’exemple
0.9152 devient 0.9152 * (1 << 16) = 59978
- L’état initial est le suivant
x = 39796
y = 0
z = 59978
z ne fait pas partie du vecteur : il sert à suivre l’angle cible restant- Le signe de
z détermine le sens de rotation
- Si
z >= 0, on effectue une rotation antihoraire et on applique z -= table[i]
- Si
z < 0, on effectue une rotation horaire et on applique z += table[i]
- Chaque itération n’utilise que des additions, des soustractions et des décalages
>> i sur x et y
if z >= 0:
x_next = x - (y >> i)
y_next = y + (x >> i)
z -= table[i]
else:
x_next = x + (y >> i)
y_next = y - (x >> i)
z += table[i]
x = x_next
y = y_next
Exemple de convergence et sujets restants
- Dans l’exemple à
0.9152 radian, z est positif à la première itération, donc on effectue une rotation antihoraire d’environ 0.785 radian
- À la deuxième itération,
z est encore positif, donc on effectue une rotation antihoraire d’environ 0.436 radian, mais on dépasse la cible
- À la troisième itération,
z devient négatif, et on effectue une rotation horaire d’environ 0.244 radian
- À la quatrième itération aussi,
z est négatif, donc on effectue une rotation horaire d’environ 0.124 radian
- À mesure que les variations d’angle diminuent, le vecteur oscille autour du résultat réel et converge
- Après 16 itérations,
y est une approximation très proche de sin(0.9152)
- L’erreur absolue sur le sinus est de
0.00000956
- L’erreur absolue sur le cosinus donné par
x est de 0.0000434
- Certains sujets restent à traiter
- Le traitement spécial nécessaire lorsque l’angle d’intérêt se trouve en dehors des 1er ou 4e quadrants du cercle unité
- Les variantes de CORDIC permettant de calculer
tan, atan, asin, acos, sinh, cosh, tanh, sqrt, ln, e^x
- L’algorithme apparenté BKM, conçu pour le calcul des logarithmes et des exponentielles
- Ces sujets seront traités plus en détail sur la chaîne YouTube Low Byte Productions
1 commentaires
Avis de Hacker News
L’auteur dit que cela s’applique surtout à des environnements comme les FPGA, mais on peut aussi l’utiliser dans le développement de jeux ou les simulations physiques distribuées
Les calculs en virgule flottante rendent la déterminisme entre plateformes difficile à garantir, et une solution consiste à éviter complètement la virgule flottante pour implémenter un moteur physique en virgule fixe
Pour implémenter les fonctions trigonométriques, il faut quelque chose comme CORDIC
Il y a quelques années, j’avais commencé à construire ce genre de chose pour m’amuser, mais je ne l’ai jamais terminé, et j’aimerais m’y remettre un jour
https://randomascii.wordpress.com/2013/07/16/floating-point-...
En résumé, x87 avait ses bizarreries ; il faut aligner de manière cohérente des réglages comme le mode d’arrondi et le flush-to-zero ; les anciens processeurs n’ont pas de FMA ; les instructions d’approximation comme
mmsqrtpsn’ont pas de spécification cohérente ; et les compilateurs peuvent réassocier les expressionsPour de petites routines ou une bibliothèque écrite à la main, il est possible, même si c’est pénible, de garantir qu’on évite ces problèmes
IEEE-754 2008 a clarifié la spécification et supposait en pratique la mort de x87 ; en 2024, on peut clairement éviter x87
Le FMA fait aussi partie de la spécification IEEE-754 2008 et se trouve dans les processeurs modernes, y compris depuis Intel Haswell
Malgré tout, des différences d’architecture comme AVX2 en 8-wide et NEON en 4-wide peuvent poser problème, mais avec de l’assembleur, des intrinsèques, ou du C vérifié dans Compiler Explorer ou avec objdump, on peut examiner la sortie et se dire : « ceci sera cohérent »
« En fait, avant qu’IEEE 754 ne devienne un standard aussi populaire qu’aujourd’hui, la virgule fixe était utilisée partout. Demandez à un développeur de jeux ayant travaillé entre 1980 et environ 2000, il vous en parlera en détail »
Rapier, la nouvelle bibliothèque réécrite à partir de nphysics, s’appuie plutôt sur les garanties d’IEEE-754 2008 pour fournir un déterminisme entre plateformes
Elle ne fonctionne donc pas sur les anciennes plateformes, mais elle est déterministe sur les plateformes modernes, y compris wasm
Bien sûr, on ne peut pas dépendre des routines de fonctions transcendantes comme sin et cos fournies par chaque plateforme ; il faut les implémenter soi-même pour qu’elles se comportent partout de la même manière
Mais si l’on n’exécute pas sur des plateformes non conformes, c’est une approche possible
https://www.rustsim.org/blog/2020/06/01/this-month-in-rustsi...
https://rapier.rs/docs/user_guides/rust/determinism/
CORDIC peut servir non seulement à calculer et générer sinus et cosinus, mais aussi à de nombreuses opérations comme les logarithmes, exponentielles, racines carrées, normes de vecteurs, conversions coordonnées polaires-cartésiennes et rotations de vecteurs
L’auteur évoque aussi ces possibilités dans la conclusion
J’ai l’impression qu’en utilisant des quaternions au lieu des matrices orthonormées habituelles, les opérations basées sur CORDIC pourraient être exécutées plus efficacement, c’est-à-dire avec moins de cycles de calcul et de mémoire, et avec moins d’erreurs
https://core.ac.uk/works/8439118
En pré-calcul au lycée, j’ai appris les séries de Taylor, et mon professeur nous a dit que les fonctions trigonométriques des calculatrices étaient réellement implémentées comme ça
En cherchant, j’ai découvert qu’en réalité c’était CORDIC, et je me suis amusé à l’implémenter en TI Basic
Ce n’était pas CORDIC, mais l’algorithme présente des similarités
http://files.righto.com/calculator/sinclair_scientific_simul...
Articles sur les implémentations matérielles :
https://arxiv.org/pdf/2211.04053
https://hal.science/hal-01327460/document
https://archive.ll.mit.edu/HPEC/agendas/proc05/Day_1/Abstrac...
J’aimerais voir comment cela se compare aux implémentations de fonctions trigonométriques logicielles et matérielles classiques sur divers matériels selon les époques.
Avec la croissance de l’IoT et des communications entre machines, et compte tenu des implémentations de CORDIC et de leur efficacité de calcul, son utilisation va probablement fortement augmenter ; il faut donc de bonnes références pour des implémentations correctes et optimisées.
Les livres du Pr Omondi et du Pr Deschamps font exception.
https://www.worldscientific.com/worldscibooks/10.1142/p1054
http://www.arithmetic-circuits.org/guide2fpga/vhdl_codes.htm
sinetcossont souvent utilisés pour les rotations de vecteurs.Dans ce cas, l’astuce avec CORDIC consiste à éviter le calcul traditionnel sin/cos/multiplication, et à fournir comme entrée à CORDIC le vecteur à faire tourner lui-même.
CORDIC produit alors directement le vecteur tourné, sans calculer
sin/cosni effectuer de multiplication de nombres complexes.CORDIC brille particulièrement lorsque la latence n’est pas cruciale.
En pipelineant chaque étape du calcul, on peut obtenir un débit élevé, ce qui convient bien au mixage numérique dans les systèmes radio.
En 2023, certains MCU modernes sont peu coûteux tout en disposant d’une FPU.
Le STM32G4 en est un bon exemple et, contrairement à des MCU comme les M0, si l’on ne veut pas utiliser de virgule fixe, on peut employer librement
f32.On peut se procurer ce type de puce pour environ 1 à 2 dollars par MCU.
Cela dit, le G4 possède aussi un périphérique matériel CORDIC qui implémente cet algorithme pour les usages en virgule fixe.
Je me demande si c’est principalement destiné à éviter les pertes de précision en virgule flottante.
On le programme via des registres, mais ce n’est pas le CPU qui implémente directement CORDIC : c’est un matériel dédié à l’intérieur de l’IC qui s’en charge.
Le STM32G4 le moins cher est le STM32G441KBT6, qui s’arrondit à 4 dollars https://www.digikey.com/en/products/detail/microchip-technol...
Je me demande où l’on peut en trouver à moins de 2 dollars.
Chez Digi-Key, la puce Nuvoton passe tout juste sous les 2 dollars pour une quantité de 500 pièces.
Il est rapide et gère des produits intermédiaires sur 64 bits, ce qui donne une précision suffisante pour la division et les fonctions trigonométriques dans la plupart des usages.
Si nécessaire, on peut encore augmenter la précision en logiciel.
J’ai découvert CORDIC tardivement ; avant cela, j’utilisais beaucoup la virgule fixe dans le monde de l’assembleur 8 et 16 bits, pour les performances et le déterminisme.
Quand je l’ai découvert, j’ai été surpris.
C’était rapide, et les notions mathématiques nécessaires pour s’en servir utilement étaient seulement basiques.
Cela me rappelle un petit morceau de code assez mignon auquel j’avais participé il y a longtemps.
Il fallait trouver les coordonnées de la bissectrice de l’angle formé par un arc du cercle unité, et on avait déjà les coordonnées
(x,y)des deux côtés.L’implémentation existante était un gros bloc de trigonométrie qui convertissait les coordonnées
(x,y)en coordonnées polaires(r,θ), vérifiait que leθcalculé était dans le bon quadrant, divisaitθpar deux, puis reconvertissait en(x,y).Résultat : beaucoup d’appels à des fonctions trigonométriques et inverses.
Comme c’était en Python et que les nombres complexes y sont des objets de première classe, il suffisait de définir deux nombres complexes
z1à partir de(x1,y1)etz2à partir de(x2,y2), puis de prendre la moyenne géométrique du produit,√(z1*z2).Le nouveau code ne contenait ni fonction trigonométrique explicite, ni conversion et conversion inverse explicites.
https://fgiesen.wordpress.com/2010/10/21/finish-your-derivat...
Il est écrit : « il est assez évident que tourner de 22,75˚ revient à tourner de 45˚ puis de -22,5˚ ». Mais dans ce cas, ne s’agit-il pas d’une rotation de 22,5° ?
Je me demande si c’est une erreur dans l’article, ou si c’est moi qui ai mal compris.
Le système d’octree de Meagher est connu pour n’utiliser que de l’arithmétique entière, sans multiplication ni division entières.
« Des algorithmes efficaces en temps linéaire ont été développés pour les opérations booléennes (union, intersection, différence), les opérations géométriques (translation, redimensionnement, rotation), la détection d’interférences en N dimensions et l’affichage incluant l’élimination des surfaces cachées depuis n’importe quel point de l’espace. Ces algorithmes ne nécessitent ni opérations en virgule flottante, ni multiplication entière, ni division entière. »
https://doi.org/10.1016/0146-664X(82)90104-6
Grâce à cela, il était plus facile de créer du matériel d’accélération graphique VLSI rapide et spécialisé pour la représentation par octree.
Je me demande quelles performances offre CORDIC par rapport à une interpolation cubique utilisant une petite table, ou à d’autres interpolations polynomiales.
J’ai appris que les synthétiseurs aux ressources limitées utilisaient parfois l’interpolation cubique, probablement à une époque où CORDIC était relativement récent.
À vue de nez, CORDIC gagne 1 bit de précision à chaque itération, donc le calcul semble plus coûteux, mais il devrait utiliser moins d’espace que les polynômes.
Cela dit, côté espace, il faut souligner qu’on peut faire encore moins cher que la table de correspondance de 4096 entrées pour
sin(x)présentée dans l’article.Grâce aux symétries, seul 1/4 du cercle complet est nécessaire.
Avec des angles sur un octet, le bouclage se faisait automatiquement, ce qui était pratique, et
2^8suffisait largement pour les rotations dans les jeux 2D.En revanche, si l’on voulait des mouvements fluides, cela n’allait pas très loin en 3D.