CORDIC

Der CORDIC-Algorithmus wurde 1959 von Jack E. Volder präsentiert. In der ursprünglichen Version war es damit möglich, trigonometrische Funktionen wie Sinus, Cosinus und Tangens sowie die Multiplikation und Division von Zahlen allein durch die in digitalen Schaltungen einfach realisierbaren Additionen und Schiebeoperationen (engl. shift-and-add operations) zu bilden. Schiebeoperationen zur Zahlenbasis 2 sind in digitalen Schaltungen sehr leicht durch entsprechende Verschaltung realisierbar.

Volders Motivation war der Ersatz der üblichen und fehleranfälligen analogen Navigationsrechner in Convair-B-58-Bombern durch digitale Rechner zur genauen Positionsbestimmung. Die Anforderung war die Positionsberechnung der mit Überschallgeschwindigkeit fliegenden Bomber in Echtzeit über einer vereinfacht als kugelförmig angenommenen Erdoberfläche.

Mitte der 1960er Jahre wurde der CORDIC-Algorithmus auch in zivilen Anwendungen eingesetzt. Vorläufer der heutigen Taschenrechner wie der Tischrechner 9100 von Hewlett-Packard aus dem Jahr 1968 setzten ihn zur Berechnung der trigonometrischen Funktionen ein.

Im Jahr 1971 wurde der CORDIC-Algorithmus von J. S. Walther auf die heute übliche Form erweitert und damit auch die effiziente Berechnung von Logarithmen, der Exponentialfunktion und der Quadratwurzel in digitalen Schaltungen möglich.

Digitalschaltung CORDIC

CORDIC-Algorithmen werden zur Berechnung der wichtigsten Elementarfunktionen in Mikrocontroller-Rechenwerken wie Taschenrechnern eingesetzt. So findet sich auch in arithmetischen x87-Koprozessoren von Intel der CORDIC-Algorithmus zur Berechnung mathematischer Operationen. Weitere Anwendungsbeispiele liegen in der Nachrichtenübertragung. Damit lassen sich beispielsweise effizient Betrag und Phase eines komplexen Signals bestimmen.

Da Multiplizierwerke vor allem in digitalen Schaltungen umfangreich und damit teuer zu realisieren sind, wird CORDIC oft genau da eingesetzt, wo Multiplizierer nicht effizient verfügbar sind. Dies umfasst vor allem den Bereich der digitalen Schaltungstechniken wie FPGAs oder ASICs.

CORDIC ist zwar nicht der schnellste Algorithmus, wird aber wegen seiner Einfachheit und Vielseitigkeit oft eingesetzt.

Vorweg wird eine Tabelle ${\displaystyle T}$ fester Länge ${\displaystyle L}$ angelegt mit ${\displaystyle T_{1}=\pi /4=\arctan \delta _{1}}$, wobei ${\displaystyle \delta _{1}=1}$ ist. Die folgenden Werte sind: ${\displaystyle T_{j}=\arctan \delta _{j}}$ mit ${\displaystyle \delta _{j}=\delta _{j-1}/2}$. (Die Werte des Arcustangens lassen sich mit der hier gut konvergierenden Potenzreihenentwicklung bestimmen.)

Die Länge ${\displaystyle L}$ der Tabelle bestimmt die erreichbare Genauigkeit. Führt man alle Drehungen eines Einheitsvektors ${\displaystyle (0,\;1)^{T}}$ mit den so berechneten Werten hintereinander in gleichem Drehsinn aus, erzielt man eine Gesamtdrehung von etwas mehr als ${\displaystyle \pm \pi /2}$. Der Skalenfaktor ${\displaystyle K}$ wird mit einem Aufruf im Vektormodus (s. u.) berechnet, indem man die Verlängerung des Einheitsvektors ${\displaystyle (0,1)^{T}}$ ohne Skalierung berechnet.

Der Ausgangsvektor ${\displaystyle (1,0)^{T}}$ wird in jedem der Schritte so gedreht, dass der Winkel ${\displaystyle z=\Theta }$ gegen Null geht. Es werden stets alle ${\displaystyle L}$ Teildrehungen ausgeführt, mit ggf. wechselndem Vorzeichen. Da der Kosinus eine gerade Funktion ist, spielt das Vorzeichen bei der Skalierung keine Rolle. Nach Reskalierung sind die Komponenten des erhaltenen Endvektors ${\displaystyle \cos \Theta }$ und ${\displaystyle \sin \Theta }$. Der Konvergenzbereich ergibt sich zu ${\displaystyle -\sum _{i=0}^{n-1}\arctan 2^{-i}\leq \Theta \leq +\sum _{i=0}^{n-1}\arctan 2^{-i}}$, also bei genügend großem ${\displaystyle n}$ etwa zu ${\displaystyle -1{,}74\leq \Theta \leq +1{,}74}$, d. h., er erstreckt sich über mehr als den vierten und ersten Quadranten.

Der vorgegebene Vektor, dessen Polarkoordinaten gesucht werden, wird immer so gedreht, dass sich der Betrag seiner ${\displaystyle y}$-Komponente verringert. Die Drehwinkel ${\displaystyle \Theta }$ werden dabei vorzeichenrichtig addiert. Die ${\displaystyle x}$-Komponente des Endvektors ist nach Reskalierung der Betrag des Ausgangsvektors. Dieser Modus wird auch benutzt zur Berechnung des Arcustangens aus zwei Argumenten, Start mit ${\displaystyle c(\cos \Theta ,\sin \Theta )^{T}}$. Der Konvergenzbereich ist derselbe wie oben. Aus ${\displaystyle \operatorname {arctan2} (y_{0},x_{0})}$ lassen sich die Funktionen ${\displaystyle \arcsin y_{0}}$ und ${\displaystyle \arccos x_{0}}$ unter Zuhilfenahme von ${\displaystyle x_{0}^{2}+y_{0}^{2}=1}$ leicht ableiten.

Der Startvektor ${\displaystyle (0,1)^{T}}$ bzw. ${\displaystyle (0,-1)^{T}}$ entspricht einer Vorwegdrehung von ${\displaystyle \pi /2}$ bzw. ${\displaystyle -\pi /2}$ (für den Rotationsmodus). Bei einem Startvektor mit negativer ${\displaystyle x}$-Komponente im Vektormodus bewirkt man entsprechende Drehungen durch Vertauschen der Komponenten und Änderungen der Vorzeichen.

Für ${\displaystyle m=0}$, ${\displaystyle \delta _{i}=t_{i}=[1,1/2,1/4,1/8,\dotsc ]}$ und ${\displaystyle k=1}$ erhält man

${\displaystyle {\begin{pmatrix}x_{i}\\y_{i}\\z_{i}\end{pmatrix}}={\begin{pmatrix}1&0&0&0\\\sigma _{i}\delta _{i}&1&0&0\\0&0&1&-\sigma _{i}\delta _{i}\end{pmatrix}}\cdot {\begin{pmatrix}x_{i-1}\\y_{i-1}\\z_{i-1}\\1\end{pmatrix}}}$,

womit sich Multiplikation und Division durchführen lassen. Eine Tabelle ${\displaystyle T}$ erübrigt sich hier.

Multiplikation

${\displaystyle x_{0}=a}$, ${\displaystyle z_{0}=b}$ ergibt im Rotationsmodus (${\displaystyle z}$ gegen 0) ${\displaystyle y_{n}=y_{0}+a\cdot b}$ für alle ${\displaystyle -2<b<2}$.

Division

${\displaystyle x_{0}=b}$, ${\displaystyle y_{0}=a}$ ergibt im Vektormodus (${\displaystyle y}$ gegen 0) ${\displaystyle z_{n}=z_{0}+a/b}$ für alle ${\displaystyle -2<a/b<2}$.

Mit ${\displaystyle m=-1}$ werden die Hyperbelfunktionen, ihre Umkehrungen (Areafunktionen), Exponentialfunktion und Logarithmus sowie die Quadratwurzel berechenbar. Einheitskreis bzw. -hyperbel werden durch ${\displaystyle {\sqrt {x^{2}+my^{2}}}=1}$ mit ${\displaystyle m=+1}$ bzw. ${\displaystyle m=-1}$ beschrieben. Das zu einem Vektor ${\displaystyle (x,y)^{T}}$ gehörende Winkel- bzw. Areaargument ist durch ${\displaystyle \Theta =\arctan({\sqrt {m}}\cdot y/x)/{\sqrt {m}}}$ gegeben, also

${\displaystyle m=1}$, Winkelfunktionen (s.o): ${\displaystyle \Theta =\arctan(y/x)}$ und

${\displaystyle m=-1}$, hyperbolische Fkt.:

${\displaystyle \Theta =\arctan(iy/x)/i}$

, hier ${\displaystyle i^{2}=-1}$; ${\displaystyle i\Theta =\arctan(iy/x)}$; und wegen ${\displaystyle \operatorname {tanh} \alpha =-i\tan(i\alpha )}$ auch ${\displaystyle \Theta =\operatorname {Atanh} (y/x)}$.

Das Verfahren ist analog zu dem eingangs gezeigten für die Winkelfunktionen. Erforderlich sind nur eine weitere Tabelle mit ${\displaystyle t_{i}=\operatorname {Atanh} \delta _{i}}$, ${\displaystyle \delta _{i}=[1/2,1/4,1/8,\dotsc ]}$ und die einmalige Berechnung des Skalenfaktors ${\displaystyle K}$.

${\displaystyle {\begin{pmatrix}x_{i}\\y_{i}\\z_{i}\end{pmatrix}}=k_{i}\cdot {\begin{pmatrix}1&\sigma _{i}\delta _{i}&0&0\\\sigma _{i}\delta _{i}&1&0&0\\0&0&1&-\sigma _{i}t_{i}\end{pmatrix}}\cdot {\begin{pmatrix}x_{i-1}\\y_{i-1}\\z_{i-1}\\1\end{pmatrix}}}$

Die Iterationen ${\displaystyle i=4,13,40,k,3k+1}$ müssen immer wiederholt werden, da der Areatangens hyperbolicus nicht die Bedingung ${\displaystyle t_{i+1}\leq {\tfrac {1}{2}}t_{i}}$ erfüllt, das somit für die Reihe ${\displaystyle \sum \sigma _{i}t_{i}}$ nicht konvergieren würde.

Rotation mit ${\displaystyle [x_{0},y_{0}]=[1,0]}$ liefert: _${\displaystyle x_{n}=\operatorname {cosh} z_{0},\quad y_{n}=\operatorname {sinh} z_{0}}$,

davon abgeleitet:

${\displaystyle \operatorname {tanh} z=\operatorname {sinh} z/\operatorname {cosh} z}$ und ${\displaystyle e^{z}=\operatorname {cosh} z+\operatorname {sinh} z}$

Vektormodus mit ${\displaystyle z_{0}=0}$ berechnet:

${\displaystyle z_{n}=\operatorname {Atanh} (y_{0}/x_{0})}$

und den hyperbolischen Betrag

${\displaystyle x_{n}={\sqrt {x_{0}^{2}-y_{0}^{2}}}}$

davon abgeleitet:

${\displaystyle \ln w=2\operatorname {Atanh} {\frac {w-1}{w+1}}}$

sowie

${\displaystyle {\sqrt {w}}}$

aus dem Betrag des Startvektors ${\displaystyle [w+1/4,w-1/4]}$

Der Konvergenzbereich ist in beiden Modi beschränkt durch die maximal mögliche Änderung von ${\displaystyle z}$. Alle mathematisch erlaubten Argumente können jedoch durch einfache Umstellungen und Shift-Operationen auf ihn abgebildet werden.