- CORDIC ist ein Algorithmus, der zur Berechnung trigonometrischer Funktionen wie
sin, cos und tan ohne FPU oder große Lookup-Tabellen komplexe Operationen in solche umwandelt, die vor allem auf Additionen und Bit-Shifts beruhen
- Dieser Ansatz ist eher in Embedded-Umgebungen nützlich als in Hochleistungssystemen, insbesondere auf leistungsschwachen Mikrocontrollern und FPGAs, und sein Wert lässt sich schwer allein an der Geschwindigkeit messen
- Verwendet man statt Gleitkommazahlen Fixed-Point-Arithmetik, kann man bei
int32_t die oberen 16 Bits als Ganzzahlteil und die unteren 16 Bits als Nachkommateil nutzen und so ungefähr Werte von -32768.99997 bis 32767.99997 darstellen
- Indem man einen Vektor mit immer kleineren Winkeln in Richtung des Zielwinkels rotiert und eine Tabelle mit 16 Werten von
atan(2**-i) sowie den Initialwert x=39796 verwendet, lässt sich die Multiplikation in jeder Iteration durch Bit-Shifts ersetzen
- Wiederholt man das Beispiel mit dem Winkel
0.9152 16-mal, sinkt der absolute Fehler von sin(0.9152) auf 0.00000956 und der absolute Fehler von cos(0.9152) auf etwa 0.0000434
Rechenumgebungen, in denen CORDIC gut passt
- CORDIC ist ein Algorithmus zur Berechnung trigonometrischer Funktionen wie
sin, cos und tan auf stromsparender Hardware
- Er funktioniert auch in Umgebungen ohne FPU, also ohne Gleitkommaeinheit, oder dort, wo große Lookup-Tabellen schwer einsetzbar sind
- Die eigentliche Berechnung besteht vor allem aus einfachen Additionen und Bit-Shifts
- Er kombiniert Vektormathematik, Trigonometrie, Konvergenz und informatische Ideen, um komplexe Funktionen durch einfache Operationen zu approximieren
- Auf Hochleistungshardware ist diese Technik nicht unbedingt nötig
- Die wichtigsten Einsatzgebiete sind Embedded-Umgebungen
- Besonders geeignet ist sie für leistungsschwache Mikrocontroller und FPGAs
- Es mag schnellere Hardware oder Peripherie geben, aber Geschwindigkeit allein ist nicht der einzige Maßstab für den praktischen Nutzen
Fixed-Point-Darstellung als Alternative zu Gleitkomma
- Auch Funktionen wie
sin(x), die Werte zwischen -1.0 und 1.0 liefern, müssen nicht zwingend als Gleitkommazahlen dargestellt werden
- Bei Fixed-Point-Arithmetik wird die Position des Kommas innerhalb eines Integer-Typs festgelegt, um rationale Zahlen darzustellen
- Im Beispiel wird
int32_t in obere 16 Bits für den Ganzzahlteil und untere 16 Bits für den Nachkommateil aufgeteilt
- In diesem Fall liegt der Bereich ungefähr zwischen
-32768.99997 und 32767.99997
- Je nachdem, wo das Komma gesetzt wird, tauscht man Bereich des Ganzzahlteils gegen Präzision des Nachkommateils ein
- Der Wert selbst bleibt weiterhin ein
int32_t; der Programmierer gibt dem Bitmuster lediglich zusätzliche Bedeutung
Umwandlung in Fixed-Point und Grundrechenarten
- Wenn der Nachkommateil 16 Bit Präzision hat, lässt sich ein Float-Wert wie
42.01 in einen Fixed-Point-Wert umwandeln, indem man ihn mit (1 << 16) multipliziert
42.01 * (1 << 16) ergibt nach dem Cast zu int32_t den Wert 2753167- Um wieder einen Float zu erhalten, berechnet man
2753167 / (1 << 16) und bekommt ungefähr 42.0099945
- Man kann auch Werte wie
1.5 direkt codieren, ganz ohne Gleitkomma zu verwenden
- Der Ganzzahlteil
1 wird mit (1 << 16) nach oben verschoben
- Die Hälfte des Nachkommateils kann als Mittelwert zwischen
0x0000 und 0xffff, also 0x7fff, gesetzt werden
- Das Ergebnis dieser Darstellung ist dezimal
98303
- Zwischen Werten mit demselben Skalierungsfaktor funktionieren Addition und Subtraktion unverändert
- Bei der Multiplikation multipliziert man zwei Fixed-Point-Werte und verschiebt das Ergebnis anschließend um den Skalierungsfaktor nach rechts
- Bei der Division kann man den Dividend vorab um den Skalierungsfaktor nach links verschieben und dann durch den Divisor teilen, um zusätzliche Präzision zu gewinnen
Trigonometrische Funktionen über Vektorrotation approximieren
- CORDIC steht für „co-ordinate rotation digital computer“ und wurde Mitte der 1950er Jahre entwickelt
- Die Kernidee ist, einen Vektor auf dem Einheitskreis mit immer kleineren Winkeln zu rotieren, sodass die Vektorkomponenten beim Erreichen des Zielwinkels den Sinus- und Kosinuswert ergeben
- Der Ablauf ähnelt einer binären Suche
- Zunächst bewegt man sich mit einem großen Winkel in Richtung des Zielwinkels
- Dann prüft man, ob man das Ziel überschritten hat
- Anschließend wiederholt man Rotationen im oder gegen den Uhrzeigersinn mit immer kleineren Winkeln
- Wenn man zum Beispiel
sin(0.7) berechnen will, startet man mit dem Vektor (1, 0) und dem Zielwinkel 0.7 Radiant
- Zuerst rotiert man um
0.7853 Radiant, also 45˚, gegen den Uhrzeigersinn
- Der verbleibende Zielwert ist dann
0.7 - 0.7853 = -0.0853
- Da der Wert negativ ist, rotiert man als Nächstes um
0.3926 Radiant, also 22.5˚, im Uhrzeigersinn
- Danach wechselt man je nach Vorzeichen des verbleibenden Zielwerts die Richtung und rotiert mit kleineren Winkeln wie
0.1963 Radiant weiter
- Nach 16 Iterationen ist der Vektor fast genau auf den ursprünglichen Zielwinkel ausgerichtet, und
y ist eine Approximation von sin(a), x eine Approximation von cos(a)
Teure Operationen in der Rotationsmatrix reduzieren
- Eine gewöhnliche Vektorrotation verwendet eine Matrixmultiplikation mit Sinus und Kosinus
- CORDIC formt die Rotationsmatrix mithilfe trigonometrischer Identitäten so um, dass
tan(a) im Zentrum steht
- Zunächst werden feste Rotationswinkel wie
45˚, 22.5˚ und 11.25˚ verwendet, sodass sich die Werte von tan(a) vorab in einer Tabelle berechnen lassen
- Diese Tabelle benötigt nur 16
uint32_t und damit 64 Byte
- Zum Vergleich: Eine nicht optimierte
sin(x)-Tabelle mit 4096 Werten von -1 bis 1 würde 16 KiB benötigen und zudem eine geringere Genauigkeit liefern
- Der bei jeder Rotation davorstehende
cos(a)-Term tritt zwar in jeder Iteration auf, aber das Produkt aller dieser Faktoren konvergiert gegen eine Konstante
- Bei Winkeln wie
45˚, 22.5˚ und 11.25˚ liegt dieses Produkt bei ungefähr 0.6366
- Diese Konstante muss nur einmal nach allen Iterationen multipliziert werden
Winkelwahl, damit nur Shifts und Additionen übrig bleiben
- Um die Multiplikation ganz zu eliminieren, wählt man Winkel so, dass das Ergebnis von
tan(a) immer eine inverse Zweierpotenz ist
- Dafür erstellt man für jede Iteration
i=0 bis 15 eine Tabelle mit 16 Einträgen der Form atan(2**-i)
- Die tatsächlichen Rotationswinkel werden dann
45˚, 26.565˚, 14.036˚, 7.125˚ usw.
- Die Winkel halbieren sich nicht exakt, doch mit ihnen konvergiert das Verfahren trotzdem korrekt
- Die Multiplikation mit
tan(a) wird durch einen Bit-Shift um i Stellen ersetzt
- Auch das Produkt der
cos(a)-Terme wird an diese neue Winkelwahl angepasst und neu berechnet
- Der Wert beträgt ungefähr
0.60725
- In 16-Bit-Fixed-Point entspricht das
39796
- Statt am Ende zu multiplizieren, kann man den Anfangswert von
x direkt auf 39796 statt auf 1 setzen
Ablauf des Algorithmus
- Im Vorberechnungsschritt erstellt man eine Tabelle, deren Einträge jeweils
atan(2**-i) sind, und wandelt jeden Wert in Fixed-Point um
- Die Umrechnungsformel lautet
atan(2**-i) * (1 << 16)
- Um
sin oder cos zu berechnen, wandelt man auch den Eingabewinkel in Fixed-Point um
- Im Beispiel wird
0.9152 zu 0.9152 * (1 << 16) = 59978
- Der Anfangszustand sieht so aus
x = 39796
y = 0
z = 59978
z ist kein Teil des Vektors, sondern verfolgt den verbleibenden Zielwinkel- Das Vorzeichen von
z bestimmt die Rotationsrichtung
- Wenn
z >= 0, rotiert man gegen den Uhrzeigersinn und führt z -= table[i] aus
- Wenn
z < 0, rotiert man im Uhrzeigersinn und führt z += table[i] aus
- Jede Iteration verwendet für
x und y nur Additionen, Subtraktionen und den Shift >> i
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
Beispiel für die Konvergenz und weitere Themen
- Im Beispiel mit
0.9152 Radiant ist z in der ersten Iteration positiv, daher erfolgt eine Rotation gegen den Uhrzeigersinn um ungefähr 0.785 Radiant
- Auch in der zweiten Iteration ist
z positiv, daher rotiert der Vektor gegen den Uhrzeigersinn um ungefähr 0.436 Radiant und überschreitet dabei das Ziel
- In der dritten Iteration ist
z negativ und der Vektor rotiert um ungefähr 0.244 Radiant im Uhrzeigersinn
- In der vierten Iteration ist
z ebenfalls negativ, daher folgt eine weitere Rotation im Uhrzeigersinn um ungefähr 0.124 Radiant
- Je kleiner die Winkeländerungen werden, desto mehr bewegt sich der Vektor nahe am tatsächlichen Ergebnis vor und zurück und konvergiert schließlich
- Nach 16 Iterationen ist
y eine sehr genaue Approximation von sin(0.9152)
- Der absolute Fehler des Sinus beträgt
0.00000956
- Der absolute Fehler des Kosinus in
x beträgt 0.0000434
- Es bleiben noch weitere Themen offen
- Spezielle Behandlung, wenn der interessierende Winkel außerhalb des ersten oder vierten Quadranten des Einheitskreises liegt
- Mit CORDIC-Varianten berechenbare Funktionen wie
tan, atan, asin, acos, sinh, cosh, tanh, sqrt, ln, e^x
- Der verwandte Algorithmus BKM, der für Logarithmen und Exponentialfunktionen entworfen wurde
- Weitere Inhalte dazu sollen auf dem Low Byte Productions YouTube channel ausführlicher behandelt werden
1 Kommentare
Hacker-News-Kommentare
Der Autor sagte, dass es hauptsächlich in Bereichen wie FPGAs eingesetzt wird, aber es lässt sich auch in der Spieleentwicklung oder für verteilte Physiksimulationen nutzen.
Bei Gleitkommarechnungen ist es schwierig, Determinismus über Plattformen hinweg sicherzustellen; eine Lösung besteht darin, Gleitkomma ganz zu vermeiden und eine Fixed-Point-Physik-Engine zu implementieren.
Um trigonometrische Funktionen zu implementieren, braucht man etwas wie CORDIC.
Vor ein paar Jahren habe ich aus Spaß angefangen, so etwas zu bauen, es aber nie fertiggestellt; irgendwann würde ich es gern noch einmal versuchen.
https://randomascii.wordpress.com/2013/07/16/floating-point-...
Kurz gesagt: x87 hatte seine Eigenheiten; Einstellungen wie Rundungsmodus und Flush-to-zero müssen konsistent gesetzt werden; ältere Prozessoren haben kein FMA; Näherungsbefehle wie
mmsqrtpshaben keine einheitliche Spezifikation; und Compiler können Ausdrücke neu assoziieren.Bei kleinen Routinen oder selbst geschriebenen Bibliotheken ist es zwar mühsam, aber möglich sicherzustellen, dass man diese Dinge vermeidet.
IEEE-754 2008 hat die Spezifikation präzisiert und faktisch den Tod von x87 vorausgesetzt; 2024 kann man x87 sicher vermeiden.
FMA ist ebenfalls Teil der IEEE-754-2008-Spezifikation und steckt in modernen Prozessoren, unter anderem seit Intel Haswell.
Trotzdem können Architekturunterschiede wie 8-wide AVX2 und 4-wide NEON zum Stolperstein werden; aber wenn man Assembly oder Intrinsics verwendet – oder C, dessen Ausgabe man mit Compiler Explorer oder objdump prüft –, kann man sich das Ergebnis ansehen und beurteilen: „Das wird konsistent sein.“
„Tatsächlich wurde Fixed-Point immer verwendet, bevor IEEE 754 zu dem heute so populären Standard wurde. Fragt einen Spieleentwickler, der etwa zwischen 1980 und 2000 gearbeitet hat; er wird euch Details erzählen.“
Rapier, die neu geschriebene Nachfolgebibliothek von nphysics, stützt sich stattdessen auf die Garantien von IEEE-754 2008, um plattformübergreifenden Determinismus zu bieten.
Deshalb läuft sie nicht auf alten Plattformen, ist aber auf modernen Plattformen einschließlich wasm deterministisch.
Natürlich kann man sich nicht auf die von den jeweiligen Plattformen bereitgestellten Routinen für transzendente Funktionen wie sin und cos verlassen; man muss sie selbst so implementieren, dass sie überall gleich funktionieren.
Wenn man sie aber nicht auf nicht konformen Plattformen ausführt, ist dieser Ansatz möglich.
https://www.rustsim.org/blog/2020/06/01/this-month-in-rustsi...
https://rapier.rs/docs/user_guides/rust/determinism/
CORDIC kann nicht nur zum Berechnen und Erzeugen von Sinus und Kosinus verwendet werden, sondern auch für viele andere Operationen wie Logarithmen, Exponentialfunktionen, Quadratwurzeln, Vektorbeträge, Umwandlungen zwischen Polar- und kartesischen Koordinaten sowie Vektorrotationen.
Der Autor deutet diese Möglichkeiten auch im Fazit an.
Ich habe das Gefühl, dass man CORDIC-basierte Operationen effizienter ausführen könnte, wenn man statt der üblichen orthonormalen Matrizen Quaternionen verwendet – also mit weniger Rechenzyklen und Speicherbedarf sowie geringeren Fehlern.
https://core.ac.uk/works/8439118
In der Analysis-Vorbereitung in der Highschool haben wir Taylorreihen gelernt, und der Lehrer sagte, dass die trigonometrischen Funktionen von Taschenrechnern tatsächlich so implementiert seien.
Als ich nachsah, stellte sich heraus, dass es in Wirklichkeit CORDIC war, und ich hatte Spaß daran, es in TI Basic nachzubauen.
Es war kein CORDIC, aber der Algorithmus hat gewisse Ähnlichkeiten.
http://files.righto.com/calculator/sinclair_scientific_simul...
Beiträge zur Hardware-Implementierung:
https://arxiv.org/pdf/2211.04053
https://hal.science/hal-01327460/document
https://archive.ll.mit.edu/HPEC/agendas/proc05/Day_1/Abstrac...
Ich würde gern sehen, wie es sich auf unterschiedlicher Hardware aus verschiedenen Epochen mit allgemeinen Software- und Hardware-Implementierungen trigonometrischer Funktionen vergleicht.
Da IoT und Machine-to-Machine-Kommunikation wachsen und man die Implementierung von CORDIC sowie seine Recheneffizienz betrachtet, wird die Nutzung vermutlich stark zunehmen; deshalb braucht es gute Referenzen für korrekte und optimierte Implementierungen.
Ausnahmen sind die Bücher von Prof. Omondi und Prof. Deschamps.
https://www.worldscientific.com/worldscibooks/10.1142/p1054
http://www.arithmetic-circuits.org/guide2fpga/vhdl_codes.htm
sin und cos werden häufig für Vektorrotationen verwendet.
In diesem Fall besteht der Trick bei CORDIC darin, die traditionelle Berechnung von sin/cos/Multiplikation zu vermeiden und den zu rotierenden Vektor selbst als Eingabe für CORDIC zu verwenden.
Dann erzeugt CORDIC den rotierten Vektor direkt, ohne sin/cos zu berechnen oder eine komplexe Multiplikation auszuführen.
CORDIC glänzt besonders dann, wenn Latenz nicht besonders wichtig ist.
Wenn man jeden Berechnungsschritt pipelinet, lässt sich ein hoher Durchsatz erzielen; das passt gut zu digitalem Mischen in Funksystemen.
Stand 2023 haben einige moderne MCUs trotz niedriger Kosten eine FPU.
Der STM32G4 ist ein gutes Beispiel; anders als etwa bei M0-MCUs kann man
f32frei verwenden, wenn man kein Fixed-Point nutzen möchte.Solche Chips bekommt man ungefähr für 1 bis 2 Dollar pro MCU.
Allerdings hat der G4 auch ein Hardware-CORDIC-Peripheriemodul, das diesen Algorithmus für Fixed-Point-Zwecke implementiert.
Ich frage mich, ob das hauptsächlich dazu dient, Präzisionsverluste bei Floating-Point zu vermeiden.
Es wird über Register programmiert, aber nicht die CPU implementiert CORDIC direkt, sondern dedizierte Hardware im IC übernimmt die Arbeit.
Der günstigste STM32G4 ist der STM32G441KBT6 und aufgerundet kostet er 4 Dollar https://www.digikey.com/en/products/detail/microchip-technol....
Ich frage mich, wo man ihn für unter 2 Dollar bekommt.
Bei Digi-Key kommt der Nuvoton-Chip bei einer Menge von 500 Stück gerade so unter 2 Dollar.
Sie ist schnell und verarbeitet 64-Bit-Zwischenprodukte, sodass die Genauigkeit bei Divisionen und trigonometrischen Funktionen für die meisten Zwecke ausreicht.
Falls nötig, kann man die Genauigkeit per Software weiter erhöhen.
Ich habe CORDIC erst spät kennengelernt; davor habe ich in der Welt von 8-Bit- und 16-Bit-Assembler viel Fixed-Point verwendet, wegen Performance und Determinismus.
Als ich es kennenlernte, war ich überrascht.
Es war schnell, und die mathematischen Fähigkeiten, die man braucht, um es nützlich einzusetzen, waren nur grundlegend.
Mir fällt ein ziemlich nettes Stück Code ein, an dem ich früher beteiligt war.
Gesucht waren die Koordinaten der Winkelhalbierenden eines Winkels, der durch einen Bogen des Einheitskreises gebildet wurde, und die
(x,y)-Koordinaten der beiden Schenkel lagen bereits vor.Die alte Implementierung war ein Haufen Trigonometrie: Sie wandelte die
(x,y)-Koordinaten in Polarkoordinaten(r,θ)um, prüfte, ob das berechneteθim richtigen Quadranten lag, halbierteθund wandelte dann wieder zurück nach(x,y).Am Ende wurden viele trigonometrische Funktionen und Umkehrfunktionen aufgerufen.
Da es Python war und komplexe Zahlen dort First-Class nutzbar sind, reichte es, aus
(x1,y1)und(x2,y2)zwei komplexe Zahlenz1undz2zu definieren und einfach das geometrische Mittel des Produkts√(z1*z2)zu nehmen.Im neuen Code gab es weder explizite trigonometrische Funktionen noch explizite Hin- und Rücktransformationen.
https://fgiesen.wordpress.com/2010/10/21/finish-your-derivat...
Dort steht: „Es ist ziemlich offensichtlich, dass eine Rotation um 22,75˚ dasselbe ist wie eine Rotation um 45˚, gefolgt von einer Rotation um -22,5˚“ — wäre das dann nicht eine Rotation um 22,5°?
Ich frage mich, ob das ein Fehler im Artikel ist oder ob ich es falsch verstanden habe.
Meaghers Octree-System ist dafür bekannt, ausschließlich Ganzzahlarithmetik ohne Integer-Multiplikation oder -Division zu verwenden.
„Es wurden effiziente Algorithmen mit linearer Laufzeit für boolesche Operationen (Vereinigung, Schnittmenge, Differenz), geometrische Operationen (Verschiebung, Skalierung, Rotation), N-dimensionale Interferenzerkennung sowie die Darstellung einschließlich Hidden-Surface Removal an beliebigen Punkten im Raum entwickelt. Diese Algorithmen benötigen keine Gleitkommaoperationen, keine Integer-Multiplikation und keine Integer-Division.“
https://doi.org/10.1016/0146-664X(82)90104-6
Dadurch war es einfacher, schnelle, maßgeschneiderte VLSI-Grafikbeschleuniger-Hardware für Octree-Darstellungen zu bauen.
Ich frage mich, wie CORDIC im Vergleich zu kubischer Interpolation mit kleinen Tabellen oder anderer polynomialer Interpolation abschneidet.
Ich habe gelernt, dass Synthesizer mit begrenzten Ressourcen manchmal kubische Interpolation verwendeten; vermutlich war das zu einer Zeit, als CORDIC noch vergleichsweise neu war.
Grob betrachtet gewinnt CORDIC pro Iteration 1 Bit an Präzision, dürfte also rechnerisch teurer sein, aber weniger Speicherplatz benötigen als Polynome.
Allerdings sollte man beim Speicherbedarf betonen, dass es günstiger sein kann als die im Artikel vorgeschlagene Lookup-Tabelle mit 4096 Einträgen für
sin(x).Dank Symmetrie braucht man nur 1/4 des gesamten Kreises.
Mit Byte-großen Winkeln war das praktisch, weil sie automatisch umliefen, und für Rotationen in 2D-Spielen war
2^8ziemlich ausreichend.Wenn man jedoch flüssige Bewegungen möchte, kommt man damit in 3D nicht besonders weit.