CORDIC-Algorithmus zur Auswertung elementarer Funktionen in Hardware

• Verfahren • Implementierung in VHDL

FH- Report

Marcel Beuler

II

Dipl.-Ing.(FH) Marcel Beuler Fachhochschule Gießen-Friedberg FB Elektro- und Informationstechnik Wiesenstraße 14 35390 Gießen

CORDIC-Algorithmus zur Auswertung elementarer Funktionen in Hardware

Version 1.2

Juni 2008

III

Inhaltsverzeichnis

1 Einführung...................................... ................1

2 Der CORDIC-Algorithmus.......................... ................2

2.1 Herleitung................................. ................2 2.2 Betriebsarten des CORDIC-Algorithmus....... ................5 2.3 Verallgemeinerung des CORDIC-Algorithmus... ................7 2.3.1 Hyperbolischer Modus (m=-1).......... ...............10 2.3.2 Linearer Modus (m=0)................. ...............13 2.3.3 Berechnung der Exponentialfunktion... ...............14 2.3.4 Erweiterung des zulässigen Argument-B ereiches für artanh(), ln() und sqrt()........ ...............15 2.3.5 Zusammenfassung CORDIC-Algorithmus... ...............18

3 Prinzipielle Möglichkeiten der Implementierung.. ...............19

3.1 Bitparallele iterative CORDIC-Struktur..... ...............19 3.2 Bitserielle iterative CORDIC-Struktur...... ...............19 3.3 Bitparallele CORDIC-Struktur, rein kombinat orisch.........20 3.4 Bitparallele CORDIC-Struktur mit Pipeline-R egistern.......21

4 Implementierung in VHDL......................... ...............22

4.1 Funktionaler Überblick..................... ...............22 4.2 Vorstufe................................... ...............23 4.3 Erweiterung des Argument-Bereiches der exp- Funktion.......25 4.4 CORDIC-Iterationen......................... ...............27 4.5 Skalierung – Konstantenmultiplizierer...... ...............30 4.6 Zusätzliche Operationen für ln-Berechnung.. ...............34 4.7 Zusätzliche Operationen für sqrt-Berechnung ...............35 4.8 Endstufe................................... ...............35 4.9 Simulation................................. ...............37

A Anhang.......................................... ...............40

A1 C-Quellcode zur Berechnung von sin() und co s()............40 A2 C-Quellcode zur Berechnung von arctan().... ...............40 A3 C-Quellcode zur Berechnung von arcsin().... ...............41 A4 C-Quellcode zur Berechnung von sinh() und c osh()..........41 A5 C-Quellcode zur Berechnung von sinh() und c osh()..........42 A6 C-Quellcode zur Berechnung von artanh().... ...............42 A7 C-Quellcode zur Multiplikation im linearen Modus..........43 A8 C-Quellcode zur Division im linearen Modus. ...............43 A9 C-Quellcode zur Berechnung von exp()....... ...............44 A10 VHDL-Quellcode für CORDIC-Algorithmus...... ...............45

Literaturverzeichnis............................... ...............66

1

1 Einführung Taylor-Polynome werden traditionell als die erste und häufig auch alleinige Möglichkeit dargestellt, elementare Funktionen auszuwerten [Riss04]. Dabei zeigen schon einfachste, ältere Taschenrechner, dass eine Potenzreihen-Entwicklung eher nicht die Methode dieser Berechnung sein kann. Eigenschaften von Intel x86 (Co-) Prozessoren deuten in dieselbe Richtung. Zunächst sei an die Taylor-Reihe einer Exponentialfunktion erinnert. Nach [Papu01] gilt hierfür:

1 2 3 4

0

11! 2! 3! 4! !

nx

n

x x x x xe

n

∞

== + + + + + =∑… (1.1)

Um zu zeigen, dass die Konvergenz-Geschwindigkeit dieser Reihe vom Argument abhängig ist, betrachten wir folgendes Zahlenbeispiel (die Taylor-Reihe wird bei n=8 abgebrochen):

0.5

5.7

1.6487212707 (Taschenrechner)

e 298.8674009671 (Taschenrechner)

e =

=

8

0.5

0

0.51.6487212650

!

n

n

en=

≈ =∑

8

5.7

0

5.7261.9926451175

!

n

n

en=

≈ =∑

Bei einer Hardware-Realisierung könnten die Nenner-Fakultäten als Konstanten in einem Speicher (Look-Up Table) abgelegt werden. Dann hat man immer noch das Problem, dass mit größer werdendem x immer mehr Glieder zur Erzielung einer vorgegebenen Genauigkeit erforderlich sind. Um die Berechnungszeit in Grenzen zu halten, müssten die Potenzen durch den Einsatz vieler Multiplizierer parallel ermittelt werden, und deren Umsetzung ist in digitalen Schaltungen bekanntlich aufwändig. Es wird deutlich, dass die Berechnung elementarer Funktionen mittels Taylor-Reihe für eine Hardware-Realisierung aus Aufwandsgründen und schlechter Abschätzung des Konvergenzverhaltens (Berechnungszeit der Exponentialfunktion würde mit dem Argument schwanken) nicht in Frage kommt. Abhilfe schaffen iterative Verfahren wie CORDIC- und Bitalgorithmus, die auf einem Tripel (x,y,z) bzw. einem Tupel (x,y) operieren. Diese Algorithmen werden üblicherweise in der ingenieurwissenschaftlichen Ausbildung nicht behandelt, obwohl sie in vielen Anwendungen wie beispielsweise Signal- oder Bildverarbeitung eine wichtige Rolle spielen. Grund genug, sich mit diesen Verfahren auseinanderzusetzen! Eine ganz besondere Bedeutung hat der hier behandelte CORDIC-Algorithmus erlangt. Er ist zwar nach [WEB01] nicht der schnellste Algorithmus, wird aber aufgrund des Bestrebens nach einer regulär und einfach strukturierten Hardware oft eingesetzt. Durch Änderung eines Parameters können hier trigonometrische und hyperbolische Funktionen sowie die Exponential-, Logarithmus- und die Wurzelfunktion berechnet werden.

2

β

( , )n nx y

0 0( , )x yϕ

y

x

r

r

β

( , )n nx y

0 0( , )x yϕ

y

x

r

r

2 Der CORDIC-Algorithmus Der CORDIC-Algorithmus wurde 1959 von J. E. Volder entwickelt und 1971 von J. S. Walther auf die heute übliche Form erweitert. CORDIC ist die Abkürzung für Coordinate Rotation Digital Computer [Pirs96]. Dieser Algorithmus ermöglicht es, elementare trigonometrische und hyperbolische Funktionen unter fast ausschließlicher Verwendung von schnellen Operationen effizient zu berechnen. Mit schnellen Operationen sind hier Additionen (adds) und Multiplikationen mit Zweier-Potenzen (shifts) gemeint. 2.1 Herleitung Ausgangspunkt für die Herleitung des CORDIC-Algorithmus ist die Drehung eines Vektors [x0,y0]

T um den Winkel φ in einem kartesischen Koordinatensystem, siehe Bild 2.1. Durch die Drehung wird der Betrag des Vektors nicht verändert.

Bild 2.1: Vektordrehung im kartesischen Koordinatensystem

00

00

cos( ) cos( )

sin( ) sin( )

xx r

ry

y rr

β β

β β

= ⇒ = ⋅

= ⇒ = ⋅

Unter Anwendung der Additionstheoreme nach [Papu01] ergibt sich:

0 0

cos( ) cos( ) cos( ) sin( ) sin( )n

x y

x r r rβ ϕ β ϕ β ϕ= ⋅ + = ⋅ ⋅ − ⋅ ⋅����� �����

0 0cos( ) sin( )nx x yϕ ϕ= ⋅ − ⋅ (2.1a)

0 0

sin( ) sin( ) cos( ) cos( ) sin( )n

y x

y r r rβ ϕ β ϕ β ϕ= ⋅ + = ⋅ ⋅ + ⋅ ⋅����� �����

0 0sin( ) cos( )ny x yϕ ϕ= ⋅ + ⋅ (2.1b) Die Gln. (2.1a) und (2.1b) lassen sich zur sog. Drehmatrix zusammenfassen. Der Ergebnisvektor [xn,yn]

T berechnet sich demnach wie folgt:

0

0

cos( ) sin( )

sin( ) cos( )n

n

x x

y y

ϕ ϕϕ ϕ

− = ⋅

(2.2)

Nach [Papu01] gilt die Identität:

3

2

1cos( )

1 tan ( )ϕ

ϕ=

+ (2.3)

Durch Ausklammern von cos(φ) lässt sich die Beziehung nach Gl. (2.2) derart modifizieren, dass innerhalb der Drehmatrix nur noch eine Abhängigkeit von einer Winkelfunktion besteht:

0 0

20 0

1 tan( ) 1 tan( )1cos( )

tan( ) 1 tan( ) 11 tan ( )

n

n

x x x

y y y

ϕ ϕϕ

ϕ ϕϕ

− − = ⋅ ⋅ = ⋅ ⋅

+ (2.4)

Bei dem CORDIC-Algorithmus wird die Drehung um den Winkel φ durch mehrere Rotationen mit bekanntem Teilwinkel αi realisiert [Pirs96]:

1

0

{ 1,1}n

i i ii

ϕ σ α σ−

== ⋅ ∈ −∑ (2.5)

Durch festzulegende Vorzeichen σi wird die Summe der n Teilwinkel αi den gegebenen Winkel φ annähern. Um n für eine vorgegebene Genauigkeit möglichst gering zu halten, wird der Betrag der Teilwinkel αi mit größerem Index abnehmen. Sollen nur einfachste Bauelemente zum Einsatz kommen, dann kann dies durch eine Folge von Zweierpotenzen erreicht werden:

tan( ) 2 iiα −= (2.6)

Weiterhin wird nach [Pirs96] eine Hilfsvariable zi eingeführt, die die aufakkumulierten Teilwinkel enthält und sich zur Steuerung der Vorzeichen der Teilwinkel verwenden lässt (Rotating-Mode):

1 arctan(2 ) sgn( ) arctan(2 )i ii i i i iz z z zσ − −+ = − ⋅ = − ⋅ (2.7)

mit: 1 0

sgn( )1 0

ii i

i

zz

zσ

+ ≥= = − <

(2.8)

CORDIC-Algorithmen operieren somit auf Tripeln (xi,yi,zi). Für das Matrixprodukt nach Gl. (2.4) ergibt sich nun:

1 1

0 0

0 00 0

1 2 1 2cos( )

2 1 2 1

i in nn i i

i i in i ii i

x x xK

y y y

σ σα

σ σ

− −− −

− −= =

− ⋅ − ⋅ = ⋅ ⋅ = ⋅ ⋅

⋅ ⋅ ∏ ∏ (2.9)

mit: 1 1 1

2 20 0 0

1 1cos( )

1 tan ( ) 1 2

n n n

i ii i ii

K αα

− − −

−= = == = =

+ +∏ ∏ ∏ (2.10)

K ist der sog. Skalierungsfaktor. Da der Kosinus eine gerade Funktion ist, spielt das Vorzeichen bei der Skalierung keine Rolle [WEB01]. Für n→∞ liefert K folgenden Wert:

1

20

1lim 0.607253

1 2

n

in i

K−

−→∞ == ≈

+∏ (2.11)

Ist n groß genug (n≥10), dann wird dieser Wert gut approximiert. Bei dem CORDIC-Algorithmus ist also am Ende des Iterationsablaufes zur Skalierung auf den korrekten Amplitudenwert noch eine abschließende Multiplikation mit dem in Gl. (2.11) angegebenen Faktor durchzuführen. Die einzelnen Teilschritte des Matrixproduktes nach Gl. (2.9) lassen sich (ohne Skalierung) durch Gl. (2.12) beschreiben. Eine solche CORDIC-Rotation besteht nach [Riss04] aus zwei SHIFTs und zwei ADDs und lässt sich daher mit minimalem Hardware-Aufwand durchführen.

4

( )

1

1

1 0

1 2 2

2 1 2

0 sgn( ) arctan 2

i ii i i i i i

i ii ii i i i

ii i i

x x x y

y y y x

z z z z

σ σσ σ

ϕ

− −+

− −+

−+

− ⋅ − ⋅ ⋅ = ⋅ =

⋅ + ⋅ ⋅

← = − ⋅ =

(2.12)

Die in [Riss04] angegebenen SHIFT- und ADD-Operationen berücksichtigen jedoch noch nicht den zu bestimmenden Arkustangens der Hilfsvariablen zi. Hierzu wird nach [WEB01] vor der eigentlichen Berechnung eine Tabelle T fester Länge L angelegt, wobei gilt:

1 1 1arctan( ) mit: 14

Tπδ δ= = = (2.13a)

1arctan( ) mit: 2 ; 2j j j jT j Lδ δ δ −= = ≤ ≤ (2.13b)

Die Länge L der Tabelle bestimmt die erreichbare Genauigkeit. Führt man alle CORDIC-Rotationen mit den so berechneten Tabellenwerten hintereinander im gleichen Drehsinn aus, so erzielt man nach Gl. (2.14) eine Gesamtdrehung von etwas mehr als ±π/2. Der Konvergenzbereich des Verfahrens erstreckt sich somit über mehr als den vierten und ersten Quadranten.

( )1

0

arctan 2 1.74n

i

i

ϕ−

−

== ± ≈ ±∑ (2.14)

Beispiel 2.1: Der Vektor [x0,y0]

T=[4,1]T soll um φ=52°=0,907571 gedreht werden (vgl. mit Bild 2.1). M ittels CORDIC-Algorithmus ist der skalierte Ergebnisvektor [x8,y8]

T mit K nach Gl. (2.11) zu bestimmen und mit der exakten Lösung zu vergleichen. Da φ>0 ist, ist das Start-Vorzeichen σ0 positiv. Mit Gl. (2.12) ergeben sich folgende Iterationsschritte:

( )

01 0 0 0

01 0 0 0

01 0 0 1

2 3

2 5

arctan 2 0.122173 0 1

x x y

y y x

z z

σ

σ

σ σ

−

−

−

= − ⋅ ⋅ =

= + ⋅ ⋅ =

= − ⋅ = > ⇒ =

( )

12 1 1 1

12 1 1 1

12 1 1 2

2 0.5

2 6.5

arctan 2 0.341475 0 1

x x y

y y x

z z

σ

σ

σ σ

−

−

−

= − ⋅ ⋅ =

= + ⋅ ⋅ =

= − ⋅ = − < ⇒ = −

( )

23 2 2 2

23 2 2 2

23 2 2 3

2 2.125

2 6.375

arctan 2 0.096496 0 1

x x y

y y x

z z

σ

σ

σ σ

−

−

−

= − ⋅ ⋅ =

= + ⋅ ⋅ =

= − ⋅ = − < ⇒ = −

( )

34 3 3 3

34 3 3 3

34 3 3 4

2 2.92188

2 6.10938

arctan 2 0.027859 0 1

x x y

y y x

z z

σ

σ

σ σ

−

−

−

= − ⋅ ⋅ =

= + ⋅ ⋅ =

= − ⋅ = > ⇒ =

( )

45 4 4 4

45 4 4 4

45 4 4 5

2 2.54004

2 6.292

arctan 2 0.03456 0 1

x x y

y y x

z z

σ

σ

σ σ

−

−

−

= − ⋅ ⋅ =

= + ⋅ ⋅ =

= − ⋅ = − < ⇒ = −

( )

56 5 5 5

56 5 5 5

56 5 5 6

2 2.73667

2 6.21262

arctan 2 0.00332 0 1

x x y

y y x

z z

σ

σ

σ σ

−

−

−

= − ⋅ ⋅ =

= + ⋅ ⋅ =

= − ⋅ = − < ⇒ = −

( )

67 6 6 6

67 6 6 6

67 6 6 7

2 2.83374

2 6.16986

arctan 2 0.012304 0 1

x x y

y y x

z z

σ

σ

σ σ

−

−

−

= − ⋅ ⋅ =

= + ⋅ ⋅ =

= − ⋅ = > ⇒ =

( )

78 7 7 7

78 7 7 7

78 7 7 8

2 2.78554

2 6.192

arctan 2 0.004492 0 1

x x y

y y x

z z

σ

σ

σ σ

−

−

−

= − ⋅ ⋅ =

= + ⋅ ⋅ =

= − ⋅ = > ⇒ =

Skalierungsfaktor nach Gl. (2.11): 8

0.607259n

K = = ⇒ 8 8

8 8

K 1.69154

3.76015

x x

K y y

⋅ → =⋅ → =

Für den exakten Ergebnisvektor [ , ]Tn nx y gilt nach Gl. (2.2):

1.67464

3.7677n

n

x

y

=

5

Beispiel 2.2: Mit dem CORDIC-Algorithmus sollen die trigonometrischen Funktionen sin(φ) und cos(φ) für φ im Bogenmaß mit |φ|<π/2 approximiert werden. Nach Gl. (2.2) lautet der Ausgangsvektor [x0,y0]

T:

0

0

1 cos( ) sin( ) 1 cos( )

0 sin( ) cos( ) 0 sin( )n

n

x x

y y

ϕ ϕ ϕϕ ϕ ϕ

− = ⇒ = ⋅ =

Der im Anhang A1 angegebene C-Quellcode führt n=15 Rotationen mit φ=0.6 durch und gibt anschließend den skalierten Ergebnisvektor aus. Man erhält folgende Werte: x = 0.82530795 // cos(0.6) y = 0.56468291 // sin(0.6) In Anlehnung an [WEB02] wird die Skalierung über g=K-1 (gain) vorgenommen (Genauigkeit: 11 Stellen nach dem Komma). Mit dem unter dieser Adresse angegebenen Dokument kann der CORDIC-Algorithmus interaktiv untersucht werden. Die von JavaScript zur Verfügung gestellten Bibliotheksfunktionen Math.sin und Math.cos dienen als Referenz (in grün). Zunächst wird im Feld “symbolic φ“ der Wert des Winkels eingegeben (hier 0.6). Anschließend wird dieser Wert durch Drücken der Schaltfläche “get φ and evaluate library functions“ übernommen. Mittels “step“ werden nun die einzelnen Iterationsschritte ausgeführt. Bild 2.2 zeigt die Ergebnisse für n=15.

Bild 2.2: Berechnung von cos(0.6) und sin(0.6) mittels interaktiven Dokuments [WEB02] Die Ergebnisse nach Bild 2.2 stimmen mit denen des C-Quellcodes (gerundet) überein. 2.2 Betriebsarten des CORDIC-Algorithmus Der bisherige Iterationsablauf des CORDIC-Algorithmus bezeichnet man als “Rotating“, da ein Vektor [x0,y0]

T um einen definierten Winkel φ gedreht wird [Pirs96]. Für die mit z0=φ initialisierte Hilfsvariable gilt hierbei z→0. In diesem Rotating-Mode lassen sich mit x0=1 und y0=0 die trigonometrischen Funktionen sin(φ) und cos(φ) berechnen, siehe Beispiel 2.2.

Arkustangensfunktion: Die zweite Betriebsart des Verfahrens ist das sog. “Vectoring“. Hier wird der Ausgangsvektor derartig gedreht, dass der Vektor nach Abschluss der Iteration auf der positiven x-Achse liegt (y→0). Dann repräsentieren die aufakkumulierten Teilwinkel den Winkel arctan(y0/x0) und die x-Koordinate den noch zu skalierenden Betrag des Vektors. Die einzelnen Teilschritte berechnen sich bis auf die Initialisierung von z0 und der Vorzeichenbestimmung σi gemäß [Andr98] analog zum Rotating-Mode nach Gl. (2.12):

6

xK e⋅ 1xe =nx

1r = argny =

xK e⋅ 1xe =nx

1r = argny =

( )

1

1

1

2

2

arctan 2

ii i i i

ii i i i

ii i i

x x y

y y x

z z

σ

σ

σ

−+

−+

−+

= − ⋅ ⋅

= + ⋅ ⋅

= − ⋅

mit: 0

1 00 ;

1 0i

ii

yz

yσ

+ <= = − >

(2.15)

Daraus ergibt sich schließlich:

�

1 2 20 0

0 00

0 00

0

arctan arctan

n

n

n

x K x y

y

y yz z

x x

−= ⋅ +=

= + =

(2.16)

Die Drehrichtung des Vektors hängt vom Quadranten ab. Mit der in Gl. (2.15) angegebenen Vorzeichensteuerung lässt sich arctan(y0/x0) für den 1. und 4. Quadranten bestimmen. Der C-Quellcode im Anhang A2 berechnet wieder n=15 Iterationen mit den Initialwerten x=1.8 und y=0.7. Gemäß Taschenrechner gilt hierfür arctan(0.7/1.8)=0.370891. Die mittels Programm berechneten Ergebnisse lauten:

x = 1.93132079 // Betrag des Vektors z = 0.37086993 // arctan(0.7/1.8)

Arkussinusfunktion: Der Vectoring-Mode eignet sich nach [WEB02] und [Andr98] auch zur Berechnung der inversen trigonometrischen Funktion φ=arcsin(arg). Hierzu muss der skalierte Einheitsvektor K·ex einer Folge von Rotationen unterzogen werden, bis die y-Komponente des gedrehten Vektors mit arg übereinstimmt. Durch die anfängliche Skalierung wird erreicht, dass nach der Rotation der Betrag des Vektors gleich Eins ist, siehe Bild 2.3:

2 2 2 2 2n1 arg x 1 argn n nr x y x= = + = + ⇒ = − (2.17)

Bild 2.3: Berechnung von arcsin(arg) Mit der Hilfsvariablen zi werden nach Bild 2.3 alle Teilwinkel zum Winkel arctan(yn/xn) aufakkumuliert. Nach Gl. (2.17) kann man hierfür auch schreiben:

2

argarctan arctan

1 arg

n

n

y

x = −

Nach [Papu01] gilt zwischen dem Arkussinus und dem Arkustangens folgender Zusammenhang:

2 2

argarcsin( ) arctan arcsin(arg) arctan

1 1 arg

xx

x

= ⇒ = − −

(2.18)

Somit entsprechen die aufakkumulierten Teilwinkel der inversen trigonometrischen Funktion φ=arcsin(arg). Es ergeben sich folgende Iterationsgleichungen:

7

( )

1

1

1

2

2

arctan 2

ii i i i

ii i i i

ii i i

x x y

y y x

z z

σ

σ

σ

−+

−+

−+

= − ⋅ ⋅

= + ⋅ ⋅

= + ⋅

mit: 0 0 0

1 (arg ) 0 ; 0 ;

1 (arg ) 0i

yx K y z

yσ

+ − >= = = = − − <

(2.19)

Die Vorzeichensteuerung berücksichtigt für (arg-y)>0 ein Anwachsen von yi+1, d.h. es gilt hier yi+1>yi. Der im Anhang A3 dargestellte C-Quellcode berechnet wieder n=15 Iterationen mit dem Initialwert arg=0.7. Hierfür gilt gemäß Taschenrechner arcsin(0.7)=0.775397. Das Programm liefert folgende Ergebnisse: x = 0.71415968 y = 0.69998283 // arg = 0.7 z = 0.77537345 // arcsin(0.7) Das Dokument [WEB02] ermöglicht auch hier wieder eine interaktive Untersuchung. Die Ergebnisse für arg=0.7 nach n=15 Iterationen zeigt Bild 2.4, sie stimmen erwartungsgemäß mit denen des C-Quellcodes (gerundet) überein. [WEB02] berechnet zusätzlich noch die arccos-Funktion, für die folgender Zusammenhang gilt:

arccos(arg) arcsin(arg)2

π= − (2.20)

Bild 2.4: Berechnung von arcsin(0.7) und arccos(0.7) mittels interaktiven Dokuments [WEB02]

2.3 Verallgemeinerung des CORDIC-Algorithmus Der CORDIC-Algorithmus in seiner bisher definierten Form erlaubt mit den beiden Betriebsarten Rotating und Vectoring die Bestimmung der trigonometrischen Funktionen sin und cos (damit auch tan) sowie der Arkusfunktionen arctan und arcsin (damit auch arccos). Durch die von J. S. Walther formulierte Verallgemeinerung des Verfahrens können zusätzliche Elementarfunktionen wie beispielsweise die hyperbolischen Funktionen und die Exponentialfunktion ermittelt werden [Walt71]. Um die Verallgemeinerung des CORDIC-Algorithmus zu verstehen, betrachten wir die hyperbolischen Additionstheoreme:

0 0

0 0

cosh( ) cosh( ) cosh( ) sinh( ) sinh( )

sinh( ) sinh( ) cosh( ) cosh( ) sinh( )

x y

y x

β ϕ β ϕ β ϕ

β ϕ β ϕ β ϕ

± = ⋅ ± ⋅

± = ⋅ ± ⋅

����� �����

����� �����

(2.21)

Diese Additionstheoreme lassen sich wieder zu einer Drehmatrix zusammenfassen:

8

0 0

0 0

cosh( ) cosh( ) sinh( ) 1 tanh( )cosh( )

sinh( ) sinh( ) cosh( ) tanh( ) 1

x x

y y

β ϕ ϕ ϕ ϕϕ

β ϕ ϕ ϕ ϕ± ± ±

= ⋅ = ⋅ ⋅ ± ± ± (2.22)

Diese Drehmatrix nach [WEB02] legt es nahe, die Funktionen cosh(φ) und sinh(φ) für ein gegebenes Argument φ durch eine Folge von hyperbolischen Rotationen zeitgleich zu berechnen. Hierzu muss φ nach [WEB02] wieder als Summe geeigneter Dreh- bzw. Teilwinkel αi und Drehrichtungen σi dargestellt werden:

( )tanh( ) 2 bzw. artanh 2i ii iα α− −= = (2.23)

( )1 1

artanh 2 { 1,1}n n

ii i i i

i i

ϕ σ α σ σ−

= == ⋅ = ⋅ ∈ −∑ ∑ (2.24)

Die Summation beginnt mit i=1, da artanh(0) nicht definiert ist. Somit lassen sich auch hier die einzelnen Rotationen – bis auf die Skalierung am Ende – ausschließlich mit SHIFT- und ADD-Operationen ausführen. Da der Kosinus hyperbolicus eine gerade Funktion ist, spielt das Vorzeichen σi bei der Bestimmung des Skalierungsfaktors keine Rolle:

( )1 1

cosh( ) cosh artanh 2n n

ii

i i

K α −

= =

= = ∏ ∏ (2.25)

Über den hyperbolischen Pythagoras kann nach [WEB02] folgende allgemeine Beziehung angegeben werden:

2 2

22

1 cosh ( ) sinh ( )1 tanh ( )

cosh( ) cosh ( )

x xx

x x

−= = −

Für den Skalierungsfaktor ergibt sich jetzt:

( ) ( )artanh 2 tanh( ) tanh artanh 2 2i i ix x− − − = ⇒ = =

2

1

1

1 2

n

ii

K− ⋅=

⇒ =−

∏ (2.26)

Für n → ∞ liefert K folgenden Wert:

� 21

1lim 1.20514

1 2

n

in i

K− ⋅→∞ =

= ≈−

∏ (2.27)

Mit β=0 erhält man: 0 00 0cosh( ) 1 ; sinh( ) 0x yβ ββ β= == = = =

Die Hilfsvariable zi enthält wieder die aufakkumulierten Teilwinkel und steuert (im Rotating-Mode) das Vorzeichen σi=sgn(zi). Damit gilt für die Folge von hyperbolischen Rotationen zur Berechnung von cosh(φ) und sinh(φ):

( )

21 1

1 0

1 2 1 2cosh( ) 1 11

sinh( ) 0 02 1 2 11 2

0 sgn( ) artanh 2

i in ni i

i iii ii i

ii i i

K

z z z z

σ σϕϕ σ σ

− −

− −− ⋅= =

−+

⋅ ⋅ = ⋅ ⋅ = ⋅ ⋅

⋅ ⋅ −

← = − ⋅

∏ ∏

ϕ=

(2.28)

9

Gl. (2.28) zeigt, dass sich auf Grundlage der hyperbolischen Additionstheoreme die hyperbolischen Funktionen cosh(φ) und sinh(φ) mittels CORDIC-Verfahren ermitteln lassen. Nun sollen die Berechnungsvorschriften auf Basis der trigonometrischen und hyperbolischen Additionstheoreme gemeinsam durch eine allgemeine Iterationsgleichung beschrieben werden. Betrachten wird zunächst die Iterationsgleichungen für beide Fälle getrennt.

trigonometrisch (Gl. (2.12)):

( )

1

1

1

2

2

arctan 2

ii i i i

ii i i i

ii i i

x x y

y y x

z z

σ

σ

σ

−+

−+

−+

= − ⋅ ⋅

= + ⋅ ⋅

= − ⋅

hyperbolisch (Gl. (2.28)):

( )

1

1

1

2

2

artanh 2

ii i i i

ii i i i

ii i i

x x y

y y x

z z

σ

σ

σ

−+

−+

−+

= + ⋅ ⋅

= + ⋅ ⋅

= − ⋅

Wir führen nun nach [Pirs96] folgende Größen ein: {1,0, 1}m∈ − Parameter für den Modus zirkular (m=1), linear (m=0) und

hyperbolisch (m=-1)

, 2 im iδ −= Schrittweite des i-ten Iterationsschrittes (0<δm,i≤1)

1mK − Skalierungsfaktor

,m iα Dreh- bzw. Teilwinkel im i-ten Iterationsschritt

Der lineare Modus (m=0) ermöglicht die Multiplikation und Division und wird in Abschnitt 2.3.2 ausführlich erläutert. Es lässt sich nun folgende allgemeine Form des CORDIC-Algorithmus angeben:

1 ,

1 ,

1 ,

1 0 0

1 0 0

0 0 11

ii i m i

ii i m i

ii i m i

xx m

yy

zz

σ δσ δ

σ α

+

+

+

−

= ⋅ −

(2.29)

Die sich für die verschiedenen Werte von m ergebenden Skalierungsfaktoren (in Gl. (2.29) nicht berücksichtigt) und Drehwinkel der einzelnen Iterationsschritte (auch Mikrorotationen genannt) sind in Tabelle 2.1 zusammengefasst.

m ,m iα 1

mK − i

1 1,arctan( )iδ 121,

0

1n

ii

δ−

=+∏

{0,1, , 1}n −…

0 0,iδ 1 {0,1, , 1}n −…

1− 1,artanh( )iδ− 21,

1

1n

ii

δ−=

−∏ (Ko

nve

{

rgenzp

1,2

rob

, ,

m

}

le )

n…

Tabelle 2.1: Teilwinkel und Skalierungsfaktoren in Abhängigkeit von m Der zirkulare Modus ist bereits in Abschnitt 2.1 und 2.2 für die beiden Betriebsarten Rotating und Vectoring ausführlich behandelt worden. Daher wird nachfolgend nur noch auf die durch die Verallgemeinerung hinzugekommenen Moden “hyperbolisch“ und “linear“ eingegangen, jeweils in beiden Betriebsarten.

10

2.3.1 Hyperbolischer Modus (m = -1) Rotating: Im letzten Abschnitt wurde durch Gl. (2.28) bereits gezeigt, dass mit dem Ausgangsvektor [x0,y0]

T=[1,0]T die hyperbolischen Funktionen cosh(φ) und sinh(φ) approximiert werden können. Der im Anhang A4 angegebene C-Quellcode führt dies mit dem Argument φ=0.8 für n=15 Iterationen aus. Gemäß Taschenrechner erhält man folgende Werte: cosh(0.8) 1.3374349463 ; sinh(0.8) 0.8881059822= = Die berechneten Ergebnisse nach Anhang A4 lauten: x = 1.33741188 // cosh(0.8) y = 0.88807125 // sinh(0.8) Gemäß [Papu01] weisen die hyperbolischen Funktionen cosh und sinh einen Definitionsbereich von −∞ bis +∞ auf. Werden alle Mikrorotationen artanh(δ-1,i) nach Tabelle 2.1 mit dem gleichen Drehsinn aufaddiert, dann erhält man allerdings folgenden Konvergenzbereich des C-Quellcode nach A4 für φ:

( )1

artanh 2 1.05547n

i

i

ϕ −

== ± ≈ ±∑ (2.30)

Das bedeutet, dass der in Gl. (2.30) angegebene Wert das maximal zulässige Argument für den Quellcode nach A4 darstellt. Im Zusammenhang mit der Exponentialfunktion wird in Abschnitt 2.3.3 gezeigt, wie mit unzulässigen Argumenten zu verfahren ist. Konvergenzproblem: Zur einfachen Realisierung des CORDIC-Verfahrens wurden bisher für die Iterationsschrittweiten δm,i Zweierpotenzen eingesetzt, wodurch sich die Matrixmultiplikation durch eine Shiftoperation implementieren lässt (siehe Tabelle 2.1):

( , ) , 2 ( , ) {0,1, , 1} mit: {1,0}S m i

m i S m i n mδ −= = − ∈… (2.31a)

( , ), 2 ( , ) {1,2, , } mit: { 1}S m i

m i S m i n mδ −= = ∈ −… (2.31b)

Die ganzzahlige Folge S(m,i) wird nach [Pirs96] als Shiftfolge des CORDIC-Algorithmus bezeichnet. Sie beeinflusst die Konvergenz, den Konvergenzbereich und den Skalierungsfaktor. Die aus ihr resultierenden Teilwinkel αm,i sollen zunächst einmal eine monoton fallende Folge bilden. Darüber hinaus muss jeder Teilwinkel einer Iteration durch eine nachfolgende Richtungsumkehr durch alle folgenden Teilwinkel bis auf einen Restfehler kompensiert werden können. Dieser ist durch den Teilwinkel der letzten Iteration definiert. Somit muss für die Winkel gelten:

1

, , , 11

{0,1, , 2} mit: {1,0}n

m i m j m nj i

i n mα α α−

−= +

− < ∈ − ∈∑ … (2.32a)

, , ,1

{1,2, , 1} mit: { 1}n

m i m j m nj i

i n mα α α= +

− < ∈ − ∈ −∑ … (2.32b)

Für eine Folge, in der der Nachfolgewinkel mindestens die Hälfte des vorhergehenden beträgt, ist Gl. (2.32) erfüllt:

11 2i iα α+ ≥ ⋅ (2.33)

Bei der zirkularen Drehung ist diese Bedingung erfüllt, denn für i=3 gilt beispielsweise:

11

( ) ( )4 312arctan 2 0.062419 arctan 2 0.062177− −= > ⋅ =

Die Verhältnisse im linearen Fall werden in Abschnitt 2.3.2 behandelt. Im hyperbolischen Modus wird Gl. (2.33) für eine ganzzahlige Shiftfolge nach Gl. (2.31b) nicht erfüllt. Mit dem C-Quellcode nach A4 soll nun gezeigt werden, wie sich diese Verletzung der Konvergenzbedingung auf die Bestimmung der Funktionen cosh(0.01) und sinh(0.01) auswirkt. Tabelle 2.2 zeigt die exakten Werte sowie die Ergebnisse für eine steigende Anzahl von Iterationsschritten (10 Stellen nach dem Komma). Man sieht, dass das Verfahren mit der gewählten Shiftfolge bei einem kleinen Argument völlig versagt.

n x (cosh(0.01)) y (sinh(0.01)) 8 1.0011095622 0.0470666476

16 1.0009314586 0.0431715729 24 1.0009308025 0.0431563596 32 1.0009307999 0.0431563001 40 1.0009307999 0.0431562999

exakt: 1.0000500004 0.0100001667 Tabelle 2.2: Ergebnisse für 0.01ϕ = bei Nichteinhaltung der

Konvergenzbedingung

Um die Konvergenzbedingung im hyperbolischen Fall zu erfüllen, muss eine Winkeldrehung an bestimmten Stellen wieder vollständig rückgängig gemacht werden können. Dies wird durch die Wiederholung der Iterationen i=4,13,40,k,3k+1 erreicht, siehe auch [Andr98, Pirs96, WEB01]. Man erhält somit die komplexere Shiftfolge nach Gl. (2.34). ( 1, ) {1,2,3,4,4,5,6, ,13,13,14,15, ,40,40,41, , }S i n− = … … … (2.34) Da die Iterationen i=4,13,40,k,3k+1 doppelt ausgeführt werden, ändert sich auch der Skalierungsfaktor. Für ihn gilt jetzt:

4 13 40 121

2 2 2 21 4 13 40

1 1 1 11.207497

1 2 1 2 1 2 1 2i i i ii i i i

K− ⋅ − ⋅ − ⋅ − ⋅= = = =

= ⋅ ⋅ ⋅ ⋅ ≈ − − − − ∏ ∏ ∏ ∏ … (2.35)

Der C-Quellcode im Anhang A5 berücksichtigt die Shiftfolge nach Gl. (2.34) und den neuen Skalierungsfaktor g=K-1 nach Gl. (2.35). Tabelle 2.3 zeigt die berechneten Ergebnisse für das Argument φ=0.01. Das CORDIC-Verfahren konvergiert jetzt gegen die exakten Werte.

n x (cosh(0.01)) y (sinh(0.01)) 8 1.0000337989 0.0079055564

16 1.0000499596 0.0099960772 24 1.0000500007 0.0100001902 32 1.0000500004 0.0100001667

exakt: 1.0000500004 0.0100001667 Tabelle 2.3: Ergebnisse für 0.01ϕ = bei Einhaltung der Konvergenzbedingung

Der Konvergenzbereich für φ erhöht sich nun gegenüber Gl. (2.30) etwas:

( ) ( ) ( ) ( )4 13 40

1

artanh 2 artanh 2 artanh 2 artanh 2 1.118173n

i

i

ϕ − − − −

=

= ± + + + + ≈ ±

∑ … (2.36)

Ist φ vom Betrag her größer als der in Gl. (2.36) angegebene Wert, so liefert der Quellcode nach A5 ebenfalls falsche Ergebnisse. Wie bereits erwähnt, wird in Abschnitt 2.3.3 gezeigt, wie mit solchen unzulässigen Argumenten im Zusammenhang mit der Exponentialfunktion zu verfahren ist.

12

Vectoring: Nach Abschluss der Iteration liegt der Ausgangsvektor wieder auf der positiven x-Achse (y→0). Die aufakkumulierten Teilwinkel repräsentieren jetzt den Winkel artanh(y0/x0) und die x-Koordinate den noch zu skalierenden hyperbolischen Betrag des Vektors. Die einzelnen Teilschritte entsprechen bis auf die Initialisierung von z0 und der Vorzeichenbestimmung σi gemäß [Andr98] dem Rotating-Mode nach Gl. (2.28):

( )

1

1

1

2

2

artanh 2

ii i i i

ii i i i

ii i i

x x y

y y x

z z

σ

σ

σ

−+

−+

−+

= + ⋅ ⋅

= + ⋅ ⋅

= − ⋅

mit: 0

1 00 ;

1 0i

ii

yz

yσ

+ <= = − >

(2.37)

Daraus ergibt sich:

�

1 2 2 1 2 20 0 0 0

0 00

0 00

0

artanh artanh

n

n

n

x K x m y K x y

y

y yz z

x x

− −= ⋅ + ⋅ = ⋅ −=

= + =

(2.38)

Die Vorzeichensteuerung nach Gl. (2.37) ermöglicht es, artanh(y0/x0) für den 1. und 4. Quadranten zu bestimmen. Gemäß [Papu01] hat die Funktion artanh einen Definitionsbereich von 1− bis 1+ und einen Wertebereich von −∞ bis +∞ . Beim CORDIC-Algorithmus ist dieser Wertebereich für φ auf den Konvergenzbereich nach Gl. (2.36) beschränkt. Der im Anhang A6 angegebene C-Quellcode berechnet n=15 Iterationen für die Werte x=3.9 und y=2.5. Der Taschenrechner liefert hierfür artanh(2.5/3.9)=0.759913. Die Ergebnisse nach A6 lauten: x = 2.99332591 // hyperbolischer Betrag z = 0.75993257 // artanh(2.5/3.9) Der maximal mögliche Quotient der Werte x0 und y0 für den Quellcode A6 beträgt:

0 0

0 0

artanh 1.118173 tanh( 1.118173) 0.806932y y

x xϕ = = ± ⇒ = ± = ±

(2.39)

Logarithmus naturalis: Über den Areatangens hyperbolicus lässt sich auch der natürliche Logarithmus mittels CORDIC-Verfahren berechnen:

( ) 12

ln ln12ln ln

1 1tanh ln ; ln ln ln

11

x x

x x

e e x x xx x x x

xx xe e

−

−

− − −= = = = = ⋅+++

1

ln 2 artanh1

xx

x

− ⇒ = ⋅ +

(2.40)

Quadratwurzel: Für r kann man auch schreiben:

( ) ( )2 2 2 21 1 1 1 1 1 1 14 4 2 16 2 16 2 2r r r r r r r r r r r+ − − = + + − + − = + = (2.41)

mit: 1 10 04 4 ; x r y r= + = −

13

2.3.2 Linearer Modus (m = 0) Der CORDIC-Algorithmus im linearen Modus multipliziert (Rotating) oder dividiert (Vectoring). Auch wenn diese Implementierung den “herkömmlichen“ Implementierungen für Multiplizierer bzw. Dividierer nicht ebenbürtig ist, so bietet sich hier doch eine attraktive Variante an, beispielsweise wenn es um die Einsparung zusätzlicher funktionaler Baugruppen geht. Im linearen Fall entspricht die Schrittweite δ0,i dem Teilwinkel α0,i im i-ten Iterationsschritt:

0, 0, 2 {0,1, , 1}ii i i nα δ −= = ∈ −… (2.42)

Der Skalierungsfaktor K0 ist 1. Aus Gl. (2.42) ist ersichtlich, dass sich im linearen Fall das System am Rande des Konvergenzbereiches befindet, da die Gleichheitsbeziehung in Gl. (2.33) gilt. Für Gl. (2.29) erhält man jetzt:

1

1 0,

1 0,

1 0 0 0

1 0 0

0 0 11

ii

ii i i

ii i i

xx

yy

zz

σ δσ δ

+

+

+

= ⋅ −

(2.43)

Rotating: In dieser Betriebsart gilt für die einzelnen Iterationsschritte:

1

1

1

2

0 2

i i

ii i i i

ii i i

x x

y y x

z z

σ

σ

+−

+−

+

=

= + ⋅ ⋅

← = − ⋅

mit: 0

1 z 00 ;

1 z 0i

ii

y σ+ >

= = − < (2.44)

Mit Gl. (2.5): 1 1

0, 00 0

2n n

ii i i

i i

zϕ σ δ σ− −

−

= == ⋅ = ⋅ =∑ ∑

Daraus resultiert schließlich unter Beachtung des Konvergenzbereiches:

�

0

0 0 0 0 0

0

0

n

n

n

x x

y y x z x z

z

== + ⋅ = ⋅

=

für: 1

00

lim 2 2n

i

n i

z−

−

→∞ == >∑ (2.45)

Der Betrag von z0 muss somit <2 sein. Im Anhang A7 ist der zugehörige C-Quellcode für die Initialwerte x=3.5 und z=1.9 angegeben. Nach n=24 Iterationen ergibt sich für das Produkt x·z: y = 6.64999992 // x*z

Vectoring: Die einzelnen Iterationsschritte lauten:

1

1

1

0 2

2

i i

ii i i i

ii i i

x x

y y x

z z

σ

σ

+−

+−

+

=

← = + ⋅ ⋅

= − ⋅

mit: 0

1 00 ;

1 0i

ii

yz

yσ

+ <= = − >

(2.46)

Wegen α0,i=δ0,i gilt für die aufakkumulierten Teilwinkel: 1 1

00,

00 0

2n n

ii i i

i i

y

xσ δ σ

− −−

= =⋅ = ⋅ =∑ ∑

14

Daraus ergibt sich:

�

0

0 00

0 00

0n

n

n

x x

y

y yz z

x x

==

= + =

für: 1

0

00

lim 2 2n

i

n i

y

x

−−

→∞ == >∑ (2.47)

Die Vorzeichensteuerung in Gl. (2.46) ermöglicht es, den Quotienten y0/x0 für den 1. und 4. Quadranten zu bestimmen. Der C-Quellcode im Anhang A8 führt die Division für die Initialwerte x=4.1 und y=-5.8 aus. Das Ergebnis nach n=24 Iterationen lautet: z = -1.41463411 // y/x 2.3.3 Berechnung der Exponentialfunktion Nach [Papu01] gilt für cosh(x) und sinh(x):

cosh( ) ; sinh( )2 2

x x x xe e e ex x

− −+ −= =

Durch Addition dieser beiden hyperbolischen Funktionen erhält man die Exponentialfunktion:

cosh( ) sinh( )2 2

x x x xx xe e e e

e e x x− −+ −+ = ⇒ = + (2.48)

Die Berechnung dieser hyperbolischen Funktionen mittels CORDIC-Algorithmus wurde in Abschnitt 2.3.1 behandelt. Als maximal zulässiges Argument für φ gilt nach Gl. (2.36):

1.118173ϕ = ± Hierfür wollen wir nachfolgend ein beliebiges Argument x ∈ℝ einsetzen, das sich als Summe eines zulässigen Argumentes z mit |z|<ln2≈0.693147 und einem Vielfachen von ln2 darstellen lässt: ln 2x z m= + ⋅ (2.49) Unter Anwendung der hyperbolischen Additionstheoreme gilt dann [WEB02]:

* **

cosh( ) cosh( ln 2) cosh( ) cosh( ln 2) sinh( ) sinh( ln 2)x z m z m z m= + ⋅ = ⋅ ⋅ + ⋅ ⋅��������� ���������

(2.50a)

( )ln 2 ln 2 2 mm me e⋅ = = ⇒

( ) ( )

( ) ( )

112

112

cosh( ln 2) 2 2 2 2 2

sinh( ln 2) 2 2 2 2 2

m m m m

m m m m

m

m

− − −

− − −

⋅ = ⋅ + = ⋅ +

⋅ = ⋅ − = ⋅ −

( )( )

1 1 1

1 1 1

cosh( ) 2 2 2 2 cosh( ) 2 cosh( )

sinh( ) 2 2 2 2 sinh( ) 2

*

sinh** ( )

m m m m

m m m m

z z z

z z z

− − − − −

− − − − −

⋅ ⋅ + = ⋅ + ⋅

⋅ ⋅ − = ⋅ − ⋅

( )1 2 cosh( ln 2) 2 cosh( ) sinh( ) 2 cosh( ) sinh( )m mz m z z z z− − ⇒ + ⋅ = ⋅ + + ⋅ − (2.50b)

15

* **

sinh( ) sinh( ln 2) sinh( ) cosh( ln 2) cosh( ) sinh( ln 2)x z m z m z m= + ⋅ = ⋅ ⋅ + ⋅ ⋅��������� ���������

(2.51a)

( )( )

1 1 1

1 1 1

sinh( ) 2 2 2 2 sinh( ) 2 sinh( )

cosh( ) 2 2 2 2 cosh( ) 2

*

cosh** ( )

m m m m

m m m m

z z z

z z z

− − − − −

− − − − −

⋅ ⋅ + = ⋅ + ⋅

⋅ ⋅ − = ⋅ − ⋅

( )1 2 sinh( ln 2) 2 cosh( ) sinh( ) 2 cosh( ) sinh( )m mz m z z z z− − ⇒ + ⋅ = ⋅ + − ⋅ − (2.51b)

Die Gln. (2.50b) und (2.51b) ermöglichen es, die hyperbolischen Funktionen cosh(x) und sinh(x) für ein beliebiges Argument x ∈ℝ zu berechnen. Nach Gl. (2.48) erhält man daher für exp(x):

[ ] [ ]1exp( ) exp( ln 2) 2 2 cosh( ) sinh( ) 2 cosh( ) sinh( )m mx z m z z z z−= + ⋅ = ⋅ ⋅ + = ⋅ + (2.52)

Der im Anhang A9 angegebene C-Quellcode berechnet die Exponentialfunktion exp(arg) mit arg∈ℝ nach Gl. (2.52). Zunächst wird die Variable arg mittels Cast-Operator (int) durch ln_2 dividiert, um den ganzzahligen Faktor m für die spätere Stellenverschiebung zu erhalten. Das Argument z, für welches die hyperbolischen Funktionen direkt bestimmt werden, geht aus Gl. (2.49) hervor: m = (int)(arg / ln_2); z = arg – m * ln_2; Bei einer Hardware-Realisierung ist z der echt-gebrochenrationale Divisionsrest von arg/ln_2. Für die Initialwerte arg=0.5, arg=5.7 und arg=-5.7 berechnet A9 nach jeweils n=32 Iterationen folgende Funktionswerte: exp_fkt = 1.6487212710 // für arg = 0.5 exp_fkt = 298.8674009805 // für arg = 5.7 exp_fkt = 0.0033459655 // für arg = -5.7 Vergleicht man die Ergebnisse für arg=0.5 und arg=5.7 mit denen des Zahlenbeispiels aus Abschnitt 1, dann erkennt man, dass das Konvergenzverhalten der Exponentialfunktion hier nicht vom Argument abhängig ist! 2.3.4 Erweiterung des zulässigen Argument-Bereic hes für artanh(), ln() und sqrt() artanh(): Das betragsmäßig größte Argument für den Areatangens hyperbolicus innerhalb des Konvergenz-Bereiches beträgt nach Gl. (2.39) ≈0.81, die Funktion weist aber einen Definitionsbereich von -1 bis +1 auf. Um auch den unzulässigen Bereich abzufangen, führen wir eine Beschreibung mittels Mantisse M und Exponent E ein, wobei gelten soll:

arg 1 2 bzw. 2 1 argE EM M− −= − ⋅ ⋅ = − (2.53)

mit: 0.5 1.0 ; 0M E≤ < ≥ Die Normalisierung der Mantisse entspricht nicht der Norm ANSI/IEEE 754, sondern ist der Original-Veröffentlichung von J. S. Walther entnommen [Walt71]. Da der ganzzahlige Exponent E nicht <0 werden kann, liegt der Definitionsbereich des Argumentes jetzt bei 0<arg<1. Für die Berechnung stellt dies kein Problem dar, weil zwischen einem positiven und einem negativen Argument folgende Beziehung besteht: artanh( arg) artanh(arg)− = − (2.54)

16

Gesucht ist somit folgende Lösung:

( )artanh 1 2 ?EM −− ⋅ = Diese Funktion wird nun so umgeformt, dass ein Term artanh() entsteht, dessen Argument sich im Konvergenzbereich befindet. Nach [Papu01] gilt für artanh(x):

12

1artanh( ) ln

1

xx

x

+ = ⋅ −

( ) 1 12 2

1 1 2 2 2 2 2 artanh 1 2 ln ln

1 1 2 2 2 2

E E EE

E E E

M MM

M M

− −−

− −

+ − ⋅ − ⋅ ⋅⇒ − ⋅ = ⋅ = ⋅

− + ⋅ ⋅ ⋅

( ) ( )1 12 2

4 2 2 4 2 2 2 artanh 1 2 ln ln

2 2

E E EE M M

MM M

−− ⋅ − − ⋅ ⋅

⇒ − ⋅ = ⋅ = ⋅

( ) ( )

( )( )

12

1 1 2 2

2 2 2 2 2artanh 1 2 ln

2 2 2 2 ln l

2 2 2 2

n 22 2 2 2

E E EE

E EE

E E

E E

M MM

M M

M M M M

M M

M

M

M

M M

M

− −−

− −

− −

− −

+ − ⋅ − ⋅ ⋅− ⋅ = ⋅ +

+ − ⋅ + − − ⋅= ⋅ + ⋅ + − ⋅ − −

+ −+ − + ⋅

⋅

−

− ⋅

( ) 12 2

2 21

2 2artanh 1 2 ln ln 22 2

12 2

E

EE EE

E

M M

M MMM M

M M

−

−−−

−

− − ⋅+ + − ⋅ − ⋅ = ⋅ + ⋅

− − ⋅− + − ⋅

( )2 artanh 1 2 artanh( ) ln 2E EM T−⇒ − ⋅ = + ⋅ (2.55)

mit: 2 2

; 0.5 1.0 ; 02 2

E

E

M MT M E

M M

−

−− − ⋅= ≤ < ≥+ − ⋅

Um den Wertebereich von T zu ermitteln, betrachten wir folgende Zahlenkombinationen von M und E: M → 1; E = 0 => T = 0 M = 0.5; E = 0 => T = 0.5 M → 1; E = 1 => T = 0.2 M = 0.5; E = 1 => T = 0.555556 M → 1; E = 100 => T = 0.3 M = 0.5; E = 100 => T = 0.6 0 0.6T⇒ < ≤ (2.56) Der Wertebereich von T liegt somit im Konvergenzbereich der artanh-Funktion nach Gl. (2.39).

Beispiel 2.3: Es ist der Areatangens hyperbolicus für die folgenden Argumente mittels Gl. (2.55) zu bestimmen:

a) artanh(0.99609375) = 3.11818480 (Taschenrechner) b) artanh(0.00390625) = 0.00390627 “ c) artanh(0.50390625) = 0.55452812 “

a)

7

72

arg 1 2 0.99609375 0.11111111

2 0.00390625 0.00000001 0.1 2

0.5 ; 7 0.59937402

artanh(0.59937402) ln 2 3.11818479

E

E

M b

M b b

M E T

−

− −

= − ⋅ = =

⋅ = = = ⋅= = ⇒ =

⇒ + ⋅ =

17

b)

02

arg 1 2 0.00390625 0.00000001

2 0.99609375 0.11111111

0.99609375 ; 0 0.00390625

artanh(0.00390625) ln 2 0.00390627

E

E

M b

M b

M E T

−

−

= − ⋅ = =

⋅ = == = ⇒ =

⇒ + ⋅ =

c)

1

12

arg 1 2 0.50390625 0.10000001

2 0.49609375 0.01111111 0.11111110 2

0.9921875 ; 1 0.20500782

artanh(0.20500782) ln 2 0.55452812

E

E

M b

M b b

M E T

−

− −

= − ⋅ = =

⋅ = = = ⋅= = ⇒ =

⇒ + ⋅ =

N

ln(): Die Logarithmus-Funktion lässt sich nach Gl. (2.40) berechnen:

1

ln 2 artanh1

xx

x

− = ⋅ +

Für das Argument der artanh-Funktion gilt wieder der Konvergenzbereich nach Gl. (2.39), d.h. für den Konvergenzbereich des ln ergibt sich zunächst:

maxmax max max

max

1 arg10.806932 arg ; 1 arg arg

1 1 arg

xx x x

x

±− = ± = ± − = ± ⋅ ± ⇒ =+ ∓

max maxmin max

max max

1 arg 1 arg0.11 ; 9.36 0.11 9.36

1 arg 1 argx x x

− += ≈ = ≈ ⇒ ≤ ≤

+ − (2.57)

Der Numerus lässt sich auch hier durch Mantisse und Exponent ausdrücken, wobei im Hinblick auf eine spätere Umsetzung in VHDL eine von [Walt71] abweichende, jedoch ANSI/IEEE-konforme Darstellung gewählt wird. Die Mantisse liegt im Konvergenzbereich von Gl. (2.57):

( )ln 2 ln( ) ln 2EM M E⋅ = + ⋅ (2.58)

mit: 1.0 2.0M≤ <

Beispiel 2.4: Es ist der natürliche Logarithmus für folgende Werte mittels Gl. (2.58) zu bestimmen:

a) ln(0.00390625) = -5.54517744 (Taschenrechner) b) ln(11.8750000) = 2.47443535 “

a)

82 0.00390625 0.00000001 1.0 2

1.0 ; 8 ln(1.0) 8 ln 2 5.54517744

EM b b

M E

−⋅ = = = ⋅= = − ⇒ − ⋅ = −

b)

32 11.875 1011.111 1.011111 2

1.484375 ; 3 ln(1.484375) 3 ln 2 2.47443535

EM b b

M E

⋅ = = = ⋅= = ⇒ + ⋅ =

N

sqrt(): Zunächst bestimmen wir den Konvergenzbereich der Wurzelfunktion nach Gl. (2.41):

1

0 41 10 0 max4 4 1

0 4

; 0.806932 argry

x r y rx r

−= + = − ⇒ = = ± = ±

+

18

( ) ( )1 1 1max max max max4 4 4arg arg 1 arg 1 argr r r− = ± ⋅ ± ⋅ ⇒ ⋅ = ± ⋅∓

max14

max

1 arg 0.03 2.34

1 argr r

±= ⋅ ⇒ < <

∓ (2.59)

Der Radikand wird auch hier ANSI/IEEE-konform durch Mantisse und Exponent beschrieben, wobei nach [Walt71] zwischen geraden und ungeraden Exponenten unterschieden werden muss. In beiden Fällen wird der Konvergenzbereich von Gl. (2.59) eingehalten (für E mod 2 = 1 wird die Erhöhung des Exponenten durch Division der Mantisse ausgeglichen):

( )1

122 22 2 2E E Er M M M= ⋅ = ⋅ = ⋅

2 2 für: 1.0 2.0

mod2 0 (gerade)

Er M M

E

⇒ = ⋅ ≤ <=

(2.60a)

( 1) 2 2 2 für: 0.5 2 1.0

mod2 1 (ungerade)

Er M M

E

+⇒ = ⋅ ≤ <=

(2.60b)

2.3.5 Zusammenfassung CORDIC-Algorithmus Nachfolgend ist der CORDIC-Algorithmus für die beiden Betriebsarten Rotating und Vectoring sowie für die drei Modi linear (m=0), zirkular (m=1) und hyperbolisch (m=-1) zusammengefasst.

Mode m Rotating

0 0 01 , 0 , 0 , sgn( )n i i

x y zz z

ϕσ

= = =→ =

Vectoring

0 0 , 0sgn( )

n

i i

z yyσ

= →= −

0

01

(0, ) {0,1,..., 1}

mKS i n

==

= −

1 0(0, )

1(0, )

1

mit: 0

2

0 2

i iS i

i i i iS i

i i i

x x y

y y x

z z

σσ

+−

+−

+

= == + ⋅ ⋅

← = − ⋅

⇒ 0

0 0

n

n

x xy x z

== ⋅

1(0, )

1(0, )

1

0 2

2

i iS i

i i i iS i

i i i

x x

y y x

z z

σσ

+−

+−

+

=← = + ⋅ ⋅

= − ⋅

⇒ 0

0 0

n

n

x xz y x

==

1

10.607252935

(1, ) {0,1,..., 1}

mKS i n

==

= −

( )

(1, )1

(1, )1

(1, )1

2

2

0 arctan 2

S ii i i i

S ii i i i

S ii i i

x x y

y y x

z z

σσ

σ

−+

−+

−+

= − ⋅ ⋅= + ⋅ ⋅

← = − ⋅

⇒1

11

1

cos( )

sin( )n

n

x K

y K

ϕϕ

−

−= ⋅= ⋅

( )

(1, )1

(1, )1

(1, )1

2

0 2

arctan 2

S ii i i i

S ii i i i

S ii i i

x x y

y y x

z z

σσ

σ

−+

−+

−+

= − ⋅ ⋅← = + ⋅ ⋅

= − ⋅

⇒1 2 2

1 0 0

0 0arctan( )n

n

x K x y

z y x

−= ⋅ +=

1

11.207497068

( 1, ) {1,2,3,4,4,5,..,13,13,14,..., }

mKS i n

−

= −=

− =

( )

( 1, )1

( 1, )1

( 1, )1

2

2

0 artanh 2

S ii i i i

S ii i i i

S ii i i

x x y

y y x

z z

σσ

σ

− −+

− −+

− −+

= + ⋅ ⋅= + ⋅ ⋅

← = − ⋅

⇒1111

cosh( )

sinh( )n

n

x K

y K

ϕϕ

−−−−

= ⋅= ⋅

( )

( 1, )1

( 1, )1

( 1, )1

2

0 2

artanh 2

S ii i i i

S ii i i i

S ii i i

x x y

y y x

z z

σσ

σ

− −+

− −+

− −+

= + ⋅ ⋅← = + ⋅ ⋅

= − ⋅

⇒1 2 21 0 0

0 0artanh( )n

n

x K x y

z y x

−−= ⋅ −

=

Tabelle 2.4: Zusammenfassung CORDIC-Algorithmus

19

3 Prinzipielle Möglichkeiten der Implementierung Die Implementierung des CORDIC-Algorithmus ist immer unter dem Gesichtspunkt bestimmter Entwurfsziele zu sehen. Diese sind beispielsweise:

• Geschwindigkeit • Hardwareaufwand (notwendige Chipfläche) • Genauigkeit • Verlustleistung

Von den zahlreichen Umsetzungsmöglichkeiten sollen 4 Verfahren erläutert werden. Aspekte wie die Umwandlung in ein anderes Datenformat (Gleit- in Festkommaformat bzw. umgekehrt) sowie die Skalierung bleiben in diesem Abschnitt unberücksichtigt. 3.1 Bitparallele iterative CORDIC-Struktur In der bitparallelen iterativen CORDIC-Struktur nach Bild 3.1 wird jede der drei Differenzengleichungen für x, y und z in Hardware realisiert. Das bisher σi genannte Vorzeichen wird dort di genannt und ist abhängig von der Betriebsart, dem Modus sowie dem Vorzeichen des y- oder z-Registers. Zu Beginn der Iteration werden die Initialwerte über Multiplexer in die x-, y- und z-Register geladen.

In jedem Iterationszyklus durchlaufen die xi und yi die Barrelshifter und Addierer/Subtrahierer und das Ergebnis wird im zugehörigen Register gespeichert. Die gewünschte Anzahl an Shiftoperationen wird, abhängig von dem Modus, in jedem Schritt angepasst. Um den durch z repräsentierten Gesamtwinkel bestimmen zu können, sind die erforderlichen Teilwinkel (arctan bzw. artanh) in jeweils einem ROM-Speicher abgelegt, dessen Adresse in jedem Iterationszyklus inkrementiert wird (in nebenstehendem Bild ist nur ein ROM-Speicher skizziert). In der letzten Iteration werden die Ergebnisse für xn, yn und zn direkt von den Addierern/Subtrahierern gelesen. Zur Ablaufsteuerung und ROM-Adressierung muss das in Bild 3.1 gezeigte Rechenwerk noch um einen Zustandsautomaten (Steuerwerk) erweitert werden. Die hier beschriebene CORDIC-Struktur führt die Shift- und arithmetischen Operationen mit Datenwortbreite innerhalb eines Taktes aus. Für die beiden Barrelshifter bedeutet diese variabel durchführbare parallele Schiebeoperation aller Bits, dass eine synthetisierte Hardware auf einem FPGA nach [Andr98] zu einer hohen Anzahl an hintereinander geschalteten Logikzellen und dementsprechend zu relativ hohen Signallaufzeiten führt.

Bild 3.1: Bitparallele iterative CORDIC- Struktur [Andr98] 3.2 Bitserielle iterative CORDIC-Struktur Der Hardwareaufwand für die Umsetzung des CORDIC-Verfahrens kann, ausgehend von der bitparallelen CORDIC-Struktur, durch Verwendung bitserieller Arithmetik weiter reduziert werden. Die einfachere Logik und Verdrahtung eines solchen bitseriellen Designs erlaubt eine wesentlich höhere Taktrate als ein vergleichbares bitparalleles Design, allerdings werden jetzt auch pro Iteration w Taktzyklen benötigt (w entspricht der Datenwortbreite). Die bitserielle iterative CORDIC-Struktur nach Bild 3.2 setzt sich zusammen aus drei bitseriellen Addierern/Subtrahierern, drei Shiftregistern und

20

einem (bzw. zwei) seriellen ROM-Speicher und benötigt somit ein Minimum an Hardware. Die Länge der Shiftregisters ist identisch mit der Wortbreite. Ihre Initialisierung kann entweder bitseriell über die linksseitig skizzierten Multiplexer in w Takten oder bitparallel in einem Takt erfolgen. Einmal initialisiert, werden die Daten nach rechts geshiftet, durchlaufen den seriellen Addierer/Subtrahierer und werden linksseitig des zugehörigen Shiftregisters wieder eingeschoben.

Bild 3.2: Bitserielle iterative CORDIC-Struktur [Andr98] Jede Iteration erfordert w Taktzyklen, um das Zwischenergebnis in das Register zurückzuschreiben. Die Funktionalität der Barrelshifter wird zum einen ersetzt durch den Rechtsshift der Register x und y in jedem Takt, zum anderen über jeweils einen Multiplexer, der in jeder Iteration entsprechend der Anzahl der Schiebeoperationen das höchstwertige Bit selektiert. Ist beispielsweise i=3, dann ist durch die beiden Multiplexer bei der Berechnung von x3 und y3 das 4. Bit von rechts das höchstwertige Bit der gekreuzten Zweige in Bild 3.2 und geht in die Addition/Subtraktion ein. Die niederwertigsten drei Bits bleiben in diesem Zweig unberücksichtigt, für die letzten drei Schiebeoperationen muss jeweils eine ’0’ verwendet werden. Diese Tatsache sowie die komplette Ablaufsteuerung werden wieder von einem (nicht skizzierten) Automaten übernommen. Während der n-ten Iteration lassen sich die Ergebnisse an den Ausgängen der seriellen Addierer/Subtrahierer lesen, während die nächsten Initialisierungsdaten in die Register geshiftet werden. 3.3 Bitparallele CORDIC-Struktur, rein kombinato risch Die beiden bisherigen Hardware-Realisierungen des CORDIC-Algorithmus sind iterativ, d.h. für die Berechnung eines Ergebnisses werden n Iterationen benötigt. Es ist jedoch auch möglich, die für jede Iterationsstufe erforderliche Hardware zu implementieren. Dieser Ansatz führt auf eine parallele, rein kombinatorische CORDIC-Struktur, die nach [Andr98] als “unrolled CORDIC processor“ bezeichnet wird und in Bild 3.3 dargestellt ist. Zwei wesentliche Punkte sind hier zu nennen: Zum einen führt jeder Shifter eine konstante Shiftoperation durch, d.h. er ist hard-wired. Zum anderen sind die Look-Up-Werte der Teilwinkel als Konstanten für jeden Addierer/Subtrahierer im z-Zweig aufgeführt. Da die Logik rein kombinatorisch arbeitet, sind keine Register erforderlich. Die Verzögerungszeit bis zum Anliegen der gültigen Ergebnisse am Ausgang ist entsprechend hoch, jedoch fällt hier im Vergleich zu Bild 3.1 die Einhaltung der Setup- und Hold-Zeit der Register zur Zwischenspeicherung weg.

21

Bild 3.3: Bitparallele CORDIC-Struktur, rein kombinatorisch [Andr98] 3.4 Bitparallele CORDIC-Struktur mit Pipeline-Re gistern Beim Einsatz von FPGAs als Zielhardware macht es keinen Sinn, eine umfangreiche kombinatorische Logik nach Bild 3.3 zu synthetisieren. Besser ist es, zwischen den Addierern/Subtrahierern Pipeline-Register vorzusehen. Da in einem FPGA jede Logikzelle ohnehin ausgangsseitig über ein Flipflop verfügt, entsteht durch zusätzliche Pipeline-Register kein Mehraufwand an Hardware. Bei einem konstanten Strom an Eingangsdaten weist diese Architektur die höchste Durchsatzrate auf. Nach einer gewissen Latenz, die abhängig ist von der Stufenzahl der Pipeline und damit der Genauigkeit, liegt in jedem Takt ein Ergebnis am Ausgang vor! Ein Beispiel für eine solche Pipeline-Architektur ist der in [Timm93] vorgestellte und für Anwendungen in der digitalen Signalverarbeitung entwickelte Chip IMSCOR24 (die im nächsten Kapitel beschriebene VHDL-Implementierung orientiert sich teilweise hieran). Extern liegen die Daten im 32-Bit-Gleitkommaformat vor, intern arbeitet der Chip mit einem 32-Bit-Festkommaformat und einer 44-stufigen Pipeline (Iterationen einschließlich Skalierung und Datenformatumwandlung). Mit diesem Chip wurden erstmals alle möglichen CORDIC-Berechnungen durchgeführt.

22

XFP[31:0] YFP[31:0] ZFP[31:0]

VorstufeX_INT[31:0] Y_INT[31:0] Z_INT[31:0]

24 CORDIC-Iterationen

QX_INT[31:0] QY_INT[31:0] QZ_INT[31:0]

SkalierungQX_MULT[32:0] QY_MULT[32:0]

Endstufe

XFP_OUT[31:0] YZFP_OUT[31:0]

Operationen fürln-Berechnung

XFP[31:0] YFP[31:0] ZFP[31:0]

VorstufeX_INT[31:0] Y_INT[31:0] Z_INT[31:0]

24 CORDIC-Iterationen

QX_INT[31:0] QY_INT[31:0] QZ_INT[31:0]

SkalierungQX_MULT[32:0] QY_MULT[32:0]

Endstufe

XFP_OUT[31:0] YZFP_OUT[31:0]

Operationen fürln-Berechnung

4 Implementierung in VHDL Der in diesem Abschnitt beschriebene strukturelle VHDL-Entwurf berücksichtigt nicht den vollständigen CORDIC-Algorithmus nach Tabelle 2.4, sondern lediglich die Berechnung der trigonometrischen Funktionen cos und sin sowie Exponential-, Logarithmus- und Wurzelfunktion. Alle implementierten Funktionen sowie der zugehörige Wertebereich der Eingabedaten sind in Tabelle 4.1 zusammengefasst. OP_CODE[7:0]

Eingabe

XFP YFP ZFP Ergebnis

XFP_OUT YZFP_OUT Wertebereich der

Eingabedaten X“01“ 1 0 z cos(z) sin(z) 0 2z π≤ ≤ X“03“ 1 0 z exp(z) 0 16z < X“04“ r r 0 0 ln(r) 1.0842022E-19 1.8446743E19r≤ ≤ X“05“ r r 0 sqrt(r) 0 9.536743E-7 2.0971519E6r≤ ≤

Tabelle 4.1: Implementierte Funktionen des CORDIC-Verfahrens 4.1 Funktionaler Überblick Bild 4.1 zeigt das Blockschaltbild für die Hardware-Realisierung. Die Eingangsdaten liegen extern im 32-Bit-Gleitkommaformat vor (1 Vorzeichenbit, 8-Bit-Charakteristik und 23-Bit-Mantisse) und werden mit XFP, YFP und ZFP bezeichnet. In der Vorstufe werden diese Eingangsdaten in das interne 32-Bit-Festkommaformat umgewandelt, wobei zur möglichst einfachen Umsetzung der Addition/Subtraktion die Zweikomplement-Darstellung gewählt wird. Für XFP und YFP wird ein gemeinsamer 10-Bit-Referenzexponent QREF_EXP_XY gespeichert. Das in Bild 4.2 gezeigte interne Festkommaformat berücksichtigt neben der um das Hidden-Bit auf 24 Bit erweiterten Mantisse noch 2 Overflow- und 5 Addition-Bits.

Bild 4.1: Blockdiagramm der Hardware-Umsetzung Die eigentliche Berechnung des Cordic-Algorithmus erfolgt in 24 Iterationen, in denen entweder die einfache Shiftfolge nach Gl. (2.31a) für den zirkularen Modus oder die komplexere Shiftfolge nach Gl. (2.34) für den hyperbolischen Modus zugrunde gelegt wird. Die Funktionswerte für arctan(2-i) bzw. artanh(2-i) sind intern in einer Tabelle (ROM-Speicher) mit 28 Nachkommastellen abgelegt.

23

S OV H 23-Bit-Mantisse Addition-Bits

31 30 28 27 05 4

S OV H 23-Bit-Mantisse Addition-Bits

31 30 28 27 05 4

ManX[23:0]

24-Bit-Barrelshifter

ManXshft[23:0] OV_BITS_X[1:0]

ManY[23:0]

24-Bit-Barrelshifter

ManYshft[23:0] OV_BITS_Y[1:0]

CONV_FP_TC

XFP[31] YFP[31]

X_INT[31:0]

CONV_FP_TCY_INT[31:0]

QX_INT[31:0] QY_INT[31:0]

XFP[30:23] YFP[30:23]

Referenzexponent

QREF_EXP_XY[9:0]

CountY[7:0] CountX[7:0]

ManZ[23:0]

24-Bit-Barrelshifter

ManZshft[23:0] OV_BITS_Z[1:0]

CONV_FP_TCZ_INT[31:0]

QZ_INT[31:0]

ZFP[30:23]

Argument-Erweiterung(exp-Funktion)

QUOTIENT[4:0]REST[23:0]

Fallunterscheidung(cos/sin-Funktion):

z <= pi/2z > pi/2 and z <= pi

z > pi and z <= 3pi/2z > 3pi/2 and z <= 2pi

QQUOTIENT[9:0]

MUX

ZFP[31]

QX_INT <= X_INT QY_INT <= Y_INT

QX_INT <= X_INT + 1QY_INT <= Y_INT - 1

QX_INT <= X_INT / 2 + 0.25QY_INT <= Y_INT / 2 - 0.25

QX_INT <= X_INT + 0.25QY_INT <= Y_INT - 0.25

OP_CODE[7:0]

X“01“ or X“03“:

X“04“:

Initialisierung oi

oi

X“05“:when VE_MOD_2=‘1‘

else

Anpassungan 32 Bit

ZFP[31]

ManX[23:0]

24-Bit-Barrelshifter

ManXshft[23:0] OV_BITS_X[1:0]

ManY[23:0]

24-Bit-Barrelshifter

ManYshft[23:0] OV_BITS_Y[1:0]

CONV_FP_TC

XFP[31] YFP[31]

X_INT[31:0]

CONV_FP_TCY_INT[31:0]

QX_INT[31:0] QY_INT[31:0]

XFP[30:23] YFP[30:23]

Referenzexponent

QREF_EXP_XY[9:0]

CountY[7:0] CountX[7:0]

ManZ[23:0]

24-Bit-Barrelshifter

ManZshft[23:0] OV_BITS_Z[1:0]

CONV_FP_TCZ_INT[31:0]

QZ_INT[31:0]

ZFP[30:23]

Argument-Erweiterung(exp-Funktion)

QUOTIENT[4:0]REST[23:0]

Fallunterscheidung(cos/sin-Funktion):

z <= pi/2z > pi/2 and z <= pi

z > pi and z <= 3pi/2z > 3pi/2 and z <= 2pi

QQUOTIENT[9:0]

MUX

ZFP[31]

QX_INT <= X_INT QY_INT <= Y_INT

QX_INT <= X_INT + 1QY_INT <= Y_INT - 1

QX_INT <= X_INT / 2 + 0.25QY_INT <= Y_INT / 2 - 0.25

QX_INT <= X_INT + 0.25QY_INT <= Y_INT - 0.25

OP_CODE[7:0]

X“01“ or X“03“:

X“04“:

Initialisierung oi

oi

X“05“:when VE_MOD_2=‘1‘

else

Anpassungan 32 Bit

ZFP[31]

Nach den 24 CORDIC-Iterationen erfolgt die Skalierung des x-Wertes (exp- und sqrt-Funktion) bzw. der x- und y-Werte (cos- und sin-Funktion) über einen Konstantenmultiplizierer in 4 Taktzyklen, wobei die Konstante abhängig ist von den beiden möglichen Shiftfolgen und das Zweierkomplement zunächst in die Vorzeichen-Betrag-Darstellung umgewandelt werden muss. Für den Logarithmus ist keine Skalierung notwendig, jedoch muss der hier berechnete artanh gemäß Gl. (2.40) durch eine Shiftoperation mit 2 multipliziert und gemäß Gl. (2.58) das Produkt E·ln2 addiert/subtrahiert werden. Anschließend wird das Zweierkomplement in die Vorzeichen-Betrag-Darstellung umgewandelt. All diese Operationen werden auf die 4 für die Konstantenmultiplikation notwendigen Taktzyklen aufgeteilt. Das Produkt E·ln2 ist intern in einem ROM-Speicher abgelegt, wobei die Adressierung über die 6 niederwertigsten Bits des Referenzexponenten QREF_EXP_XY erfolgt.

Bild 4.2: Internes 32-Bit-Festkommaformat für QX_INT, QY_INT und QZ_INT 4.2 Vorstufe Der interne Aufbau der Vorstufe ist in Bild 4.3 dargestellt. Da die Addition/Subtraktion von Gleitkommazahlen nur für gleiche Exponenten direkt möglich ist, wird zunächst für XFP und YFP ein Referenzexponent QREF_EXP_XY ermittelt und abgespeichert. Hierbei handelt es sich um den größeren der beiden Exponenten, d.h. der kleinere Operand wird auf den Exponenten des größeren Operanden umgeformt. Zu diesem Zweck ist für jeden der beiden Operanden ein 24-Bit-Barrelshifter als kombinatorische Logik vorgesehen, siehe hierzu auch [Beul08]. Das zugehörige VHDL-File Barrelshifter.vhd steht im Anhang A10. Die Differenz zwischen dem größeren und dem kleineren Exponenten, also die Anzahl der Rechtsverschiebungen, wird hier mit dem Signalvektor COUNT[7:0] vorgegeben. Das zusätzlich noch vorhandene Bit SHIFT_DIRECTION (siehe weiter unten) muss für einen Rechtsshift auf 0 gesetzt sein. Die geshifteten Mantissen werden mit ManXshft[23:0] bzw. ManYshft[23:0] bezeichnet.

Bild 4.3: Vorstufe

24

Für ZFP sind die Verhältnisse etwas komplizierter, weil die Funktionswerte für arctan()bzw. artanh() intern in einer Tabelle im 32-Bit-Festkommaformat nach Bild 4.2 abgelegt sind, d.h. es handelt sich um denormalisierte Mantissen und der Exponent (ohne Bias) ist 0. Ist der Exponent von ZFP, d.h. ZFP(30 downto 23), größer als 127d, dann erfolgt ein Linksshift um die Differenzstellenzahl “ZFP(30 downto 23) – 127“, andernfalls ein Rechtsshift um die Differenzstellenzahl “127 – ZFP(30 downto 23)“. Um zwischen Links- und Rechtsshift unterscheiden zu können, besitzt der 24-Bit-Barrelshifter das Bit SHIFT_DIRECTION. Bei SHIFT_DIRECTION = ’0’ erfolgt ein Rechtsshift, bei ’1’ entsprechend ein Linksshift um die über COUNT angegebene Stellenzahl. Die geshiftete Mantisse wird mit ManZshft[23:0] bezeichnet. Ihr zugehöriger Referenzexponent ist immer 0, so dass dieser nicht explizit abgespeichert werden muss. Die 24-Bit-Barrelshifter liefern zusätzlich noch die notwendigen Overflow-Bits OV_BITS_X[1:0] OV_BITS_Y[1:0] und OV_BITS_Z[1:0]. Bei einem Rechtsshift sind diese OV-Bits “00“ (für XFP und YFP immer der Fall), bei einem Linksshift um eine bzw. zwei Stellen werden die geshifteten Bits in die OV-Bits übernommen. Argument-Erweiterung der exp-Funktion: Im hyperbolischen Modus erfolgt zur Erweiterung des Argument-Bereiches bei der Berechnung der Exponentialfunktion eine Division von ZFP durch ln(2), bei der ein ganzzahliger Quotient, QUOTIENT[4:0], und ein echt-gebrochenrationaler Rest, REST[23:0], entsteht. Eine ausführliche Erläuterung dieser Zusammenhänge erfolgt im nächsten Abschnitt. Die eigentliche Umwandlung der geshifteten Mantissen ManXshft[23:0], ManYshft[23:0] und ManZshft[23:0] in das interne 32-Bit-Festkommaformat in Zweikomplement-Darstellung erfolgt in den Blöcken CONV_FP_TC, siehe Anhang A10. Ist beispielsweise die Gleitkommazahl XFP positiv, d.h. XFP(31) = ’0’, dann werden XFP(31), OV_BITS_X[1:0] sowie ManXshft[23:0] direkt in das Festkommaformat X_INT bzw. QX_INT übernommen (der Index Q deutet an, dass es sich um ein Flipflop bzw. Register handelt). Die noch fehlenden Addition-Bits werden alle auf 0 gesetzt. Ist die Gleitkommazahl hingegen negativ, dann wird nur XFP(31) direkt übernommen. Alle Stellen der OV-Bits und der geshifteten Mantisse werden invertiert und die Addition-Bits alle auf 1 gesetzt. Um das Komplement zu erhalten, wird abschließend noch eine 1 addiert. Fallunterscheidung für cos und sin: Der Argument-Bereich für die trigonometrischen Funktionen cos und sin soll zwischen 0 und 2π liegen. Damit das Verfahren konvergiert, darf φ nach Gl. (2.14) den Wert 1.74 nicht überschreiten. Somit muss bei Auswahl des OP-Codes X“01“ das Argument (arg=Z_INT) im Automatenzustand Z1 gemäß nachfolgender Tabelle modifiziert sowie der 2-Bit-Signalvektor FALLUNTERSCHEIDUNG entsprechend gesetzt werden:

Argument FALLUNTER-SCHEIDUNG

modifiziertes Argument

cos-Zweig (x)

sin-Zweig (y)

arg 2π≤ 00 QZ_INT Z_INT<= cos sin

arg 2 and argπ π> ≤ 01 QZ_INT Z_INT 2π<= − sin− cos

arg and arg 3 2π π> ≤ 10 QZ_INT Z_INT π<= − cos− sin−

arg 3 2 and arg 2π π> ≤ 11 QZ_INT Z_INT 3 2π<= − sin cos−

Tabelle 4.2: Argument-Modifizierung der trigonometrischen Funktionen Die CORDIC-Berechnung erfolgt mit dem modifizierten Argument wie bisher. Im Zustand Z8 werden die berechneten Ergebnisse über das Signal FALLUNTERSCHEIDUNG an den cos- bzw. sin-Zweig durch Umkopieren und/oder Vorzeicheninvertierung (falls notwendig) angepasst, siehe auch [Walt71]. Logarithmus und Quadratwurzel: Soll mittels OP-Code X“04“ der Logarithmus berechnet werden, dann steht nach Tabelle 4.1 der Numerus sowohl in XFP als auch in YFP. Nach Gl. (2.40) wird im x-Zweig eine 1 addiert und im y-Zweig eine 1 subtrahiert. Diese Operationen erfolgen im Zustand Z1. Für die Berechnung der Quadratwurzel mittels OP-Code X“05“ muss der Radikand ebenfalls in XFP und YFP stehen (ein negativer Radikand wird nicht abgefangen). Nach Gl. (2.60) muss zunächst geprüft werden, ob der Exponent (egal ob von XFP oder YFP, da beide Werte ja identisch) gerade oder ungerade ist. Hierzu wird im Prozess A_BEST im Zustand Z1 die Variable VE_MOD_2 mit der jeweiligen Bitstelle des ROM-Speichers E_MOD_2 gesetzt. In der derzeitigen Implementierung berücksichtigt diese interne Tabelle nicht alle möglichen 8-Bit-Exponenten, sondern nur Werte

25

zwischen 107d und 147d. Ausgehend von der Dezimalzahl 1 kann das Komma damit um 20 Stellen nach links bzw. rechts verschoben werden. Ist der Exponent ungerade, dann wird die Mantisse durch 2 geteilt, andernfalls bleibt sie unverändert. Der Teilung durch 2 entspricht ein Rechtsshift, bei welchem nach dem Vorzeichenbit von links eine ’0’ eingeschoben wird. 4.3 Erweiterung des Argument-Bereiches der exp-F unktion Nach Gl. (2.48) setzt sich die Exponentialfunktion additiv aus den hyperbolischen Funktionen cosh und sinh zusammen, für deren maximal zulässiges Argument Gl. (2.36) gilt. Für ein beliebiges Argument x ∈ℝ erhält man nach Gl. (2.52):

exp( ) exp( ln 2) 2 [cosh( ) sinh( )]mx z m z z= + ⋅ = ⋅ + Um die Exponentialfunktion berechnen zu können, müssen also zunächst die Größen z und m aus x bestimmt werden. Dies führt auf eine Division mit einem ganzzahligen Quotienten (m) und einem echt-gebrochenrationalen Rest (z): ln 2z x m= − ⋅ (4.1) Um die Vorgehensweise zu verdeutlichen, gehen wir zunächst von einer reduzierten Stellenanzahl mit einer 10-Bit-Mantisse (einschließlich Hidden-Bit) aus und betrachten allgemein folgendes Divisionsbeispiel: x = 9.171875000d = 1001.001011 = 1.001001011 * 2 ^3 ln2 = 0.693359375d = 0.101100011 Der Dezimalpunkt des Dividenden wird nun um eine Stelle nach links verschoben, wodurch sich seine Stellenzahl um ein Bit erhöht: x = 0.1001001011 * 2^4 Nun erfolgt die Division mit Rückstellen des Zwischenrestes (siehe z.B. [Beul08]): Stufe 0: 01001001011 : 0101100011 = 0.1101 * 2^4 = 1101 = 13d -0101100011 0 11110000101 => 0

Stufe 1: 10010010110 -0101100011 0 00111010000 => 1

Stufe 2: 01110100000 -0101100011 0 00011011010 => 1

Stufe 3: 00110110100 -0101100011 0 11011101110 => 0

Stufe 4: 01101101000 -0101100011 0 00010100010 => 1

� x = 1101 = 13d ; z = 0.001010001 = 0.158203125d

Die hier beschriebene Umsetzung des Divisionsverfahrens über 5 (stellenverschobene) Subtraktionen erlaubt einen maximalen Exponenten (ohne Bias) von E=3. Damit wird das maximal mögliche Argument für die Exponentialfunktion beschränkt auf:

26

3max 1111.1111111111 1.1111111111111 2 16x = = ⋅ <… …

16 6 16 7 8.886 10 ; 1.125 10e e− −⇒ = ⋅ = ⋅ Trotz dieser Beschränkung können schon relativ große bzw. kleine Ergebnisse bestimmt werden. Mit der Verwendung einer weiteren Subtraktionsstufe ließe sich dieser Zahlenbereich noch erheblich erweitern auf dann exp(32) bzw. exp(-32). In dem gezeigten Rechenbeispiel ist der Exponent E=3, so dass die 4. Stufe den Divisionsrest repräsentiert. Ist dieses Ergebnis negativ (höchstwertige Stelle eine ’1’), dann erfolgt eine Rückstellung des Zwischenrestes ohne anschließende Stellenverschiebung. Der Wertebereich des Exponenten ist nach oben auf E=3 und nach unten zunächst auf E=-1 beschränkt. Eine in diesem Wertebereich durchgeführte Division hat immer einen echt-gebrochenrationalen Rest zur Folge, welcher vom Exponenten E abhängig ist, siehe Tabelle 4.3.

E Rest -1 Stufe 0 0 Stufe 1 1 Stufe 2 2 Stufe 3 3 Stufe 4

Tabelle 4.3: Abhängigkeit des Restes von E Für -1≤E≤3 kann dieser Rest unmittelbar als Argument z verwendet werden, da er einen Referenzexponenten von 0 aufweist. Ist E>3, so führt dies zu falschen Ergebnissen! Der im Anhang angegebene VHDL-Quellcode fängt diese Wertebereichsüberschreitung nicht ab. Ist E<-1 (|x|<0.5) dann muss als Divisionsrest die geshiftete Mantisse ManZshft[23:0] verwendet werden, da für sie der Referenzexponent 0 ist. Hardwaremäßig werden alle 5 Subtraktionsstufen kombinatorisch realisiert. Die korrekten Bits für den Quotienten sowie der zugehörige Divisionsrest werden über den Exponenten E ermittelt. Sind hier nicht alle 5 Stellen zur Darstellung des Quotienten erforderlich, dann werden die höherwertigen Bits auf ’0’ gesetzt. Der 5-stellige ganzzahlige Quotient wird ausgehend von dem Vorzeichenbit ZFP(31) im Automatenzustand z1 in das 10-stellige Signal QQUOTIENT[9:0] in Zweikomplement-Darstellung umgerechnet und fließt später im Zustand z8 bei der Bestimmung des Ergebnis-Exponenten VREF_EXP_X der Exponentialfunktion mit ein. Abschließend sei noch ein Rechenbeispiel für E=2 angegeben: x = 5.156250000d = 101.0010100 = 1.010010100 * 2 ^2 = 0.1010010100 * 2^3 Stufe 0: 01010010100 : 0101100011 = 0.111 * 2^3 = 111 = 7d -0101100011 0 11111001110 => 0

Stufe 1: 10100101000 -0101100011 0 01001100010 => 1

Stufe 2: 10011000100 -0101100011 0 00111111110 => 1

Stufe 3: 01111111100 -0101100011 0 00100110110 => 1

� x = 111 = 7d ; z = 0.010011011 = 0.302734375d

27

4.4 CORDIC-Iterationen Nach der Vorstufe werden, ausgehend von dem 32-Bit-Festkommaformat nach Bild 4.2, 24 Iterationen ausgeführt. Da für dieses Format das Zweierkomplement gewählt wurde, sind die arithmetischen Operationen “Addition“ und “Subtraktion“ besonders einfach. Die Vorgehensweise soll der Vollständigkeit halber für beide Operationen anhand von typischen Zahlenbeispielen demonstriert werden, auch wenn die Abfrage der hierbei auftretenden Bereichsüberschreitungen im VHDL-Code nicht implementiert ist. Addition: Nachfolgend aufgeführte Beispiele für die Addition von 4-stelligen Zweikomplementzahlen sind [Hoff93] entnommen. C OV C OV C OV a1) +2 0010 a2) +7 0111 b 1) -1 1111 +2 0010 +7 0111 -1 1111 +4 0 0100 0 0 +14 01110 0 1 -2 1 1110 1 0

b2) -7 1001 c1) +2 0010 c 2) +1 0001 -7 1001 -1 1111 -7 1001 -14 10010 1 1 +1 1 0001 1 0 -6 0 1010 0 0 Tabelle 4.4 zeigt die sich aus den Zahlenbeispielen ergebenden Fälle der Addition im Zweierkomplement.

Fall A n-1 B n-1 S n S n-1 Summe a1 0 0 0 0 S[n] = A + B a2 0 0 0 1 S[n+1] = A + B “Überlauf“ b1 1 1 1 1 S[n] = A + B – 2^n b2 1 1 1 0 S[n+1] = A + B “Überlauf“ c1 0 1 1 0 S[n] = A + B – 2^n c2 0 1 0 1 S[n] = A + B

Tabelle 4.4: Die verschiedenen Fälle der Addition im Zweikomplement In den Fällen a2 und b2 findet eine Bereichsüberschreitung (Überlauf, Overflow) statt, da das Ergebnis nicht mehr mit n=4 Stellen darstellbar ist. Die hinzukommende Summenstelle Sn kann dann als Vorzeichenstelle interpretiert werden. Die Stelle Sn entspricht dem Übertrag (carry out). Der Übertrag, der in den Fällen b1 und c1 auftritt, gehört nicht zur Darstellung des n -stelligen Ergebnisses und wird deshalb abgeschnitten (A+B-2n). Alle geschilderten Fälle können nach [Hoff93] in zwei Schritten zusammengefasst werden:

1. S[n]=A+B bilden und den Übertrag Sn merken. Die Überlaufbedingung berechnen. 2. Wenn kein Überlauf aufgetreten ist (OV=0), dann ist die Summe gleich S[n], ansonsten

S[n+1]. Subtraktion: Auch die Subtraktion von Zweikomplementzahlen soll anhand von Zahlenbeispielen für n=4 veranschaulicht werden. C OV C OV C OV a1) (+7) 0111 a2) (+5) 0101 b1) (+2) 0010 -(+3) 0011 -(+7) 0111 -(-5) 1011 +4 0 0100 0 0 -2 1 1110 1 0 +7 1 0111 1 0

b2) (+3) 0011 c1) (-5) 1011 c2) (-5) 1011 -(-6) 1010 -(+3) 0011 -(+6) 0110 +9 11001 1 1 -8 0 1000 0 0 -11 00101 0 1

28

C OV C OV d1) (-6) 1010 d2) (-3) 1101 -(-4) 1100 -(-7) 1001 -2 1 1110 1 0 +4 0 0100 0 0 Die sich ergebenden Fälle sind wieder in Tabelle 4.5 zusammengefasst.

Fall A n-1 B n-1 D n D n-1 Summe a1 0 0 0 0 D[n] = A - B a2 0 0 1 1 D[n] = A – B – 2^n b1 0 1 1 0 D[n] = A - B – 2^n b2 0 1 1 1 “Überlauf“ → Wert zu groß c1 1 0 0 1 D[n] = A - B c2 1 0 0 0 “Überlauf“ → Wert zu klein d1 1 1 1 1 D[n] = A - B – 2^n d2 1 1 0 0 D[n] = A - B

Tabelle 4.5: Die verschiedenen Fälle der Subtraktion im Zweikomplement Hier findet in den Fällen b2 und c2 eine Bereichsüberschreitung aufgrund einer zu geringen Stellenzahl statt. Invertiert man die hinzukommende Differenzstelle Dn, dann kann dieses Bit als Vorzeichenstelle interpretiert werden. Der Übertrag in den Fällen a2, b1 und d1 kann wieder abgeschnitten werden, da er nicht zum Ergebnis gehört (A-B-2n). Auch die Subtraktion im Zweikomplement kann in zwei Schritten zusammengefasst werden:

1. D[n]=A-B bilden und den Übertrag Dn merken. Die Überlaufbedingung berechnen. 2. Wenn kein Überlauf aufgetreten ist (OV=0), dann ist die Differenz gleich D[n], ansonsten wird

das Bit Dn invertiert und die Differenz ist D[n+1]. Wie bereits zu Beginn dieses Abschnitts erwähnt, fängt der im Anhang A10 aufgeführte VHDL-Quellcode die hier angegebenen Bereichsüberschreitungen sowohl bei Addition als auch bei Subtraktion nicht ab. Bei allen eingegebenen Testdaten ist ein solcher Fall (fehlerhaftes Ergebnis) bisher nicht aufgetreten! CORDIC-Architektur: Als Architektur kommt das Prinzip der in Bild 3.1 gezeigten bitparallelen iterativen CORDIC-Struktur zur Anwendung. Es werden im Prozess A_BEST im Zustand Z3 24 CORDIC-Iterationen ausgeführt. Nach der 24. Iteration wechselt der Automat in den Zustand Z4. Als Laufindex für den aktuellen CORDIC-Iterationszyklus dient der 5-Bit-Vektor QS[4:0]. Der Index i der Shiftfolgen wird im Quellcode über die Variable N ausgedrückt, die in jeder Iteration die Anzahl der Rechtsverschiebungen für die beiden bitparallelen Barrelshifter im x- und y-Zweig angibt. Sie ist vom Typ integer, da mit ihrem Wert innerhalb einer for-loop-Anweisung eine Umsortierung der Signalvektoren und damit die eigentliche Shiftoperation erfolgt. Zwischen QS und N besteht in Abhängigkeit von m folgender Zusammenhang:

QS (dezimal) 0 1 2 3 4 5 6 7 8 9 10 11 N (m=1) 0 1 2 3 4 5 6 7 8 9 10 11 N (m=-1) 1 2 3 4 4 5 6 7 8 9 10 11

QS (dezimal) 12 13 14 15 16 17 18 19 20 21 22 23

N (m=1) 12 13 14 15 16 17 18 19 20 21 22 23 N (m=-1) 12 13 13 14 15 16 17 18 19 20 21 22

Tabelle 4.6: Zusammenhang zwischen QS und N in Abhängigkeit von m Nachfolgender Codeausschnitt zeigt die Bestimmung von N aus QS gemäß Tabelle 4.6. when Z3 => if m = 1 then -- einfache Shiftfolge für zirkular en Modus N:=Bit5Int(QS); else -- komplexere Shiftfolge für hyperb olischen Modus if QS < "00100" then

29

N:=Bit5Int(QS) + 1; elsif QS < "01110" then N:=Bit5Int(QS); else N:=Bit5Int(QS) - 1; end if ; end if ; Die Ergebnisse der jeweiligen Shiftoperationen für den x- bzw. y-Zweig stehen in den Variablen VX_SHIFT[31:0] und VY_SHIFT[31:0]. Um die CORDIC-Berechnungen durchzuführen, ist für jeden Datenpfad eine zusätzliche 32-Bit-Variable vorgesehen, namentlich VX_INT, VY_INT und VZ_INT. Auf diese Weise kann durch Zugriff auf die Bitstelle VX_INT(31) bzw. VY_INT(31) das Vorzeichen oi für die nachfolgende CORDIC-Iteration ermittelt werden. Die Teilwinkel arctan(2-i) und artanh(2-i) sind jeweils als ROM-Speicher (ROM_arctan bzw. ROM_artanh) mit der Organisation 24*32 in der Quellcodedatei PACK_CORDIC.vhd abgelegt. Jedes 32-Bit-Datenwort setzt sich zusammen aus einem Vorzeichenbit (immer positiv, d.h. ’0’) sowie 3 Vor- und 28 Nachkommastellen. Die Umwandlung der Dezimalzahlen in die entsprechenden Dualzahlen erfolgt sowohl hier als auch bei dem Konstantenmultiplizierer im nächsten Abschnitt über das Programm dez_bin.exe der Website http://delphi.zsg-rottenburg.de/download.html. Im hyperbolischen Modus startet die Shiftfolge mit dem Wert 1, die Adressierung des Speichers ROM_artanh beginnt aber ab der Adresse 0. Deshalb wird hier von der Variablen N der Wert 1 subtrahiert. Die Iterationsgleichungen in Abhängigkeit von m und oi können nun direkt mittels Tabelle 2.4 angegeben werden: when Z3 => if m = 1 then -- zirkularer Modus if oi = 1 then VX_INT := QX_INT - VY_SHIFT; VY_INT := QY_INT + VX_SHIFT; VZ_INT := QZ_INT - ROM_arctan(N) ; else VX_INT := QX_INT + VY_SHIFT; -- oi ist -1 VY_INT := QY_INT - VX_SHIFT; VZ_INT := QZ_INT + ROM_arctan(N) ; end if ; else -- hyperbolischer Modus if oi = 1 then VX_INT := QX_INT + VY_SHIFT; VY_INT := QY_INT + VX_SHIFT; VZ_INT := QZ_INT - ROM_artanh(N- 1); else VX_INT := QX_INT - VY_SHIFT; -- oi ist -1 VY_INT := QY_INT - VX_SHIFT; VZ_INT := QZ_INT + ROM_artanh(N- 1); end if ; end if ;

Abschließend findet die Vorzeichenbestimmung oi aus VZ_INT(31) bzw. VY_INT(31) für die nächste Iteration statt und die Variablenwerte werden in die Register QX_INT, QY_INT und QZ_INT übernommen. Der Abbruch des CORDIC-Verfahrens nach 24 Iterationen und damit der Wechsel in den Zustand Z4 erfolgt im Prozess F_BEST: when Z3 => if QS = "10111" then FOLGE_Z <= Z4; -- Abbruch des COR DIC-Verfahrens nach 24 Iterationen end if ;

Über den Prozess Z_SPEICHER wird QS im Zustand Z0 zurückgesetzt und im Zustand Z3 mit jeder steigenden Flanke um 1 inkrementiert.

30

4.5 Skalierung – Konstantenmultiplizierer Eine dem CORDIC-Verfahren innewohnende Eigenschaft ist die Skalierung der x- und y-Werte in Abhängigkeit von der verwendeten Shiftfolge, siehe hierzu die Gln. (2.11) und (2.35). Da hier ein Faktor konstant ist, kann die Hardware-Realisierung der Multiplikation vereinfacht werden. Bevor eine solche Konstantenmultiplikation ausgeführt wird, erfolgt die Umwandlung des nicht konstanten Faktors (QX_INT bzw. QY_INT; im Folgenden Dualzahl LD genannt) von der Zweikomplement- in die Vorzeichen-Betrag-Darstellung. Da die Konstanten zur Skalierung positiv sind, ist das Produktvorzeichen identisch mit dem Vorzeichen für LD und kann am Ende der Multiplikation einfach übernommen werden. Als einführendes Beispiel betrachten wir die Multiplikation einer 12-Bit-Dualzahl V[11:0] mit einer 8-Bit-Konstanten K[7:0]. Das Ergebnis ist eine 20-Bit-Dualzahl P[19:0]:

( )( )( )( )

11 10 9 8 2 1 011 10 9 8 2 1 0

11 10 9 8 11 10 9 8

7 6 5 47 6 5 4

3 2 1 03 2 1 0

2 2 2 2 2 2 2

2 2 2 2

2 2 2 2

2 2 2 2

P K v v v v v v v

K v v v v

K v v v v

K v v v v

= ⋅ ⋅ + ⋅ + ⋅ + ⋅ + + ⋅ + ⋅ + ⋅

= ⋅ ⋅ + ⋅ + ⋅ + ⋅

+ ⋅ ⋅ + ⋅ + ⋅ + ⋅

+ ⋅ ⋅ + ⋅ + ⋅ + ⋅

…

Die folgende Struktur aus [Jork04] zeigt die Bildung der Teilprodukte und ihre Zusammenfassung zum Gesamtprodukt: k7 k6 k5 k4 k3 k2 k1 k0 * v0 k7 k6 k5 k4 k3 k2 k1 k0 * v1 k7 k6 k5 k4 k3 k2 k1 k0 * v2 k7 k6 k5 k4 k3 k2 k1 k0 * v3 --------------- ------------------ u11u10 u9 u8 u7 u6 u5 u4 u3 u2 u1 u0

k7 k6 k5 k4 k3 k2 k1 k0 * v4 k7 k6 k5 k4 k3 k2 k1 k0 * v5 k7 k6 k5 k4 k3 k2 k1 k0 * v6 k7 k6 k5 k4 k3 k2 k1 k0 * v7 --------------------------- ------ r11r10 r9 r8 r7 r6 r5 r4 r3 r2 r1 r0

k7 k6 k5 k4 k3 k2 k1 k0 * v8 k7 k6 k5 k4 k3 k2 k1 k0 * v9 k7 k6 k5 k4 k3 k2 k1 k0 * v10 k7 k6 k5 k4 k3 k2 k1 k0 * v11 --------------------------------- s11s10 s9 s8 s7 s6 s5 s4 s3 s2 s1 s0

Die Zwischensummen von je vier Partialprodukten werden zu einem Modul zusammengefasst. Jede Teilsumme kann in Abhängigkeit von den vier Bitstellen der Variablen 16 verschiedene Werte mit einer Wortlänge von 12 Bit annehmen. Ein solches Modul lässt sich als ROM-Speicherblock mit der Größe 16*12 Bit bilden. Dieses Prinzip soll nun zur Skalierung der Ergebnisse des CORDIC-Verfahrens angewandt werden. Zunächst können wir einen allgemein gültigen Zusammenhang zwischen den Längen der Konstanten und der Dualzahl, der Anzahl der ROM-Speicher und deren jeweiliger Kapazität (Organisation) sowie der Anzahl der Adressbits angeben: Länge der Konstanten: KL

Länge der Dualzahl: DL

Anzahl der ROM-Speicher: ROMN

Anzahl der Adressbits: ROMA

Anzahl der Datenbits: ROMD

Kapazität des Speichers: 2 ROMAROM ROMK D= ∗

31

31 025 24

0

0

0

0

0

025

25

25

25

25

25

∑

Vorkommastellen Nachkommastellen

32-Bit-ROM-Datenwort31 025 24

0

0

0

0

0

025

25

25

25

25

25

∑

Vorkommastellen Nachkommastellen

32-Bit-ROM-Datenwort

D ROM ROML N A⇒ = ⋅ (4.2)

( ) 2 2ROM ROMA AROM ROM K ROMK D L A= ∗ = ∗ + (4.3)

Die Dualzahl LD weist nach Bild 4.2 zunächst 31 Stellen (ohne Vorzeichenbit) auf. Um eine möglichst sinnvolle Einteilung zu erhalten, wird LD auf 30 Stellen reduziert, d.h. das niederwertigste Addition-Bit fällt weg. Mit 5 ROM-Speichern ergeben sich auf diese Weise nach Gl. (4.2) 6 Adressbits: 30 5 6D ROM ROML N A= = ⋅ = ⋅ Weiterhin wird als Datenwortbreite DROM eine Zweierpotenz angestrebt, was auf eine Konstantenlänge von LK=26 Bit führt, siehe Bild 4.4. Für die Darstellung eines Skalierungsfaktors von ≈1.2 ist eine Vorkommastelle ausreichend. Ein ROM-Speicher muss somit nach Gl. (4.3) folgende Kapazität aufweisen:

( ) 6 2 2 (26 6) 64 32ROMAROM K ROMK L A= ∗ + = ∗ + = ∗

Bild 4.4: Datenwortbreite des ROM-Speichers In Tabelle 4.7 ist ein kleiner Ausschnitt eines ROM-Speichers für die Skalierung mit K≈0.607252 (grau hinterlegte Felder) aufgeführt. Die vollständigen Tabellen für beide Skalierungsfaktoren (zirkular und hyperbolisch) sind im Package-File PACK_CORDIC.vhd im Anhang A10 zu finden. Bild 4.5 zeigt schließlich den sich ergebenden Konstantenmultiplizierer, bestehend aus 5 ROM-Speichern und 4 Addierern. Das Produkt (ohne Vorzeichen) umfasst LK+LD=56 Bits, wobei die höchstwertigen 4 Bits die Vorkommastellen sind. Die 6 niederwertigsten Bits des 1. ROMs sind bereits die 6 niederwertigsten Ergebnisbits. Jede der ersten drei Addierstufen liefert weitere 6 Bits zum Endergebnis. Die Summe der letzten Addierstufe stellt die 32 höchstwertigen Bits des Produktes dar. Eingang SUM SUM (dezimal) SUM (binär) 000000 0K 00.000000000000000 0000000.0000000000000000000000000 000001 1K 00.607252925634384 0000000.1001101101110100111011011 000010 2K 01.214505851268768 0000001.0011011011101001110110110 000011 3K 01.821758776903152 0000001.1101001001011110110010001 000100 4K 02.429011702537536 0000010.0110110111010011101101100 000101 5K 03.036264628171920 0000011.0000100101001000101000111 000110 6K 03.643517553806304 0000011.1010010010111101100100010 000111 7K 04.250770479440688 0000100.0100000000110010011111101

… … … … … … … …

111111 63K 38.256934314966192 0100110.0100000111000110011100101

Tabelle 4.7: Ausschnitt eines ROM-Speichers für 0.607252K ≈

32

ROM_1[31:0]ROM 64*32E[6:1]

ROM 64*32E[12:7]

ROM 64*32E[18:13]

ROM 64*32E[24:19]

ROM 64*32E[30:25]

ROM_2[31:0]

ROM_3[31:0]

ROM_4[31:0]

ROM_5[31:0]

ROM_1[31:6]

ADD_1[31:6]

“000000“

“000000“

ADD_2[31:6]

“000000“

“000000“

ADD_3[31:6]

ADD

ADD

ADD

ADD

E[31:0]

ADD_1[5:0]

ROM_1[5:0]

ADD_2[5:0]

ADD_4[31:0]

ADD_3[5:0]

A[55:0]

ROM_1[31:0]ROM 64*32E[6:1]

ROM 64*32E[12:7]

ROM 64*32E[18:13]

ROM 64*32E[24:19]

ROM 64*32E[30:25]

ROM_2[31:0]

ROM_3[31:0]

ROM_4[31:0]

ROM_5[31:0]

ROM_1[31:6]

ADD_1[31:6]

“000000“

“000000“

ADD_2[31:6]

“000000“

“000000“

ADD_3[31:6]

ADD

ADD

ADD

ADD

E[31:0]

ADD_1[5:0]

ROM_1[5:0]

ADD_2[5:0]

ADD_4[31:0]

ADD_3[5:0]

A[55:0]

ROM 64*32E[6:1]

ROM 64*32

A[31:0]E[12:7]

E[18:13]

E[24:19]

E[30:25]

MUX

“000000“

ADDADD[31:0]

MUX

REGREG[31:0]

REG[31:6]E[31:0]

X[31:6]

Y[31:0]

ROM 64*32E[6:1]

ROM 64*32

A[31:0]E[12:7]

E[18:13]

E[24:19]

E[30:25]

MUX

“000000“

ADDADD[31:0]

MUX

REGREG[31:0]

REG[31:6]E[31:0]

X[31:6]

Y[31:0]

Bild 4.5: Kombinatorische Realisierung der Konstantenmultiplikation

Der Aufwand an ROM-Speichern und Addierern lässt sich reduzieren, wenn die Operationen nicht kombinatorisch, sondern sequentiell ausgeführt werden. Der Konstantenmultiplizierer nach Bild 4.6 kommt neben zwei zusätzlichen Multiplexern mit nur 2 ROM-Speichern und einem Addierer aus. Weiterhin ist noch eine nicht skizzierte Ablaufsteuerung erforderlich. Diese Lösung erfordert 4 Taktzyklen für eine Multiplikation und ist im VHDL-Quellcode im Anhang A10 implementiert.

Bild 4.6: Sequentielle Realisierung der Konstantenmultiplikation Zur Simulation mit ModelSim XE werden die ROM-Speicher im Package-File als Array deklariert: type MMX64_32 is array (63 downto 0) of std_logic_vector(31 downto 0);

Im zirkularen Modus müssen für die Kosinus- und Sinusfunktion sowohl der x- als auch der y-Wert skaliert werden, d.h. in 4 Taktzyklen erfolgen parallel zwei Konstantenmultiplikationen. Somit sind für den Skalierungsfaktor nach Gl. (2.11) insgesamt vier ROM-Speicher notwendig: constant K1_ROM_X1,K1_ROM_X2,K1_ROM_Y1,K1_ROM_Y2: MMX64_32:= -- zirkularer Modus ("01001100100000111000110011100101", -- 111111 "01001011010011001010001100001010", -- 111110 ... ... "00000000000000000000000000000000"); -- 000000

33

Im hyperbolischen Modus muss für die Berechnung der Exponential- und Quadratwurzelfunktion lediglich der x-Wert skaliert werden, d.h. hier reichen für den Skalierungsfaktor nach Gl. (2.35) zwei ROM-Speicher aus: constant K2_ROM_X1,K2_ROM_X2: MMX64_32:= -- hyperbolischer Modus ("10011000001001010000011001110010", -- 111111 "10010101101110101100100101100100", -- 111110 ... ... "00000000000000000000000000000000"); -- 000000

In Bild 4.6 fällt auf, dass nicht alle 56 Ergebnisstellen verwendet werden, sondern die niederwertigsten 26 Stellen wegfallen. Die nicht skizzierte Ablaufsteuerung wird im Prozess A_BEST des Top-Level-Moduls durch die Automatenzustände Z4 bis Z7 übernommen. In Z4 wird zunächst das Zweierkomplement in die Vorzeichen-Betrag-Darstellung umgewandelt. Lediglich bei der Berechnung der exp-Funktion müssen vor dieser Umwandlung die Register QX_INT (cosh) und QY_INT (sinh) addiert werden. Um in Z4 direkt auf die Ergebnisse der Umwandlung zugreifen zu können, werden diese den 32-Bit-Variablen VX_INT bzw. VY_INT zugewiesen. Zusätzlich erfolgt am Prozessende die Speicherung in QX_INT bzw. QY_INT. Das in Bild 4.6 aufgeführte Register wird im VHDL-Entwurf um ein Bit auf 33 Bit zur Aufnahme des Vorzeichens erweitert und mit QX_MULT bzw. QY_MULT bezeichnet. Die serielle Konstantenmultiplikation sieht dann für den zirkularen Modus (OP_CODE = X“01“) prinzipiell wie folgt aus: when Z4 => -- OP_CODE = X“01“ -- Umwandlung des Zweierkomplements in die Vorzeichen-Betrag-Darstellung if QX_INT(31) = '1' then VX_INT(31 downto 0) := QX_INT(31) & ( not (QX_INT(30 downto 0)) + 1); else VX_INT(31 downto 0) := QX_INT; end if ; if QY_INT(31) = '1' then VY_INT(31 downto 0) := QY_INT(31) & ( not (QY_INT(30 downto 0)) + 1); else VY_INT(31 downto 0) := QY_INT; end if ; QX_INT <= VX_INT(31 downto 0); -- unskalierte Vorzeichen-Betrag-Darstellung QY_INT <= VY_INT(31 downto 0); -- unskalierte Vorzeichen-Betrag-Darstellung --Beginn der seriellen Konstantenmulti plikation V_ROM_X1_OUT := K1_ROM_X1(Bit6Int(VX_I NT(6 downto 1))); V_ROM_Y1_OUT := K1_ROM_Y1(Bit6Int(VY_I NT(6 downto 1))); QX_MULT <= QX_INT(31) & (("000000" & V _ROM_X1_OUT(31 downto 6)) + K1_ROM_X2(Bit6Int(VX_INT(12 downto 7)))); QY_MULT <= QY_INT(31) & (("000000" & V _ROM_Y1_OUT(31 downto 6)) + K1_ROM_Y2(Bit6Int(VY_INT(12 downto 7)))); when Z5 => QX_MULT(31 downto 0) <= ("000000" & QX_MULT(31 downto 6)) + K1_ROM_X2(Bit6Int(QX_INT(18 downto 13))); QY_MULT(31 downto 0) <= ("000000" & QY_MULT(31 downto 6)) + K1_ROM_Y2(Bit6Int(QY_INT(18 downto 13))); when Z6 => QX_MULT(31 downto 0) <= ("000000" & QX_MULT(31 downto 6)) + K1_ROM_X2(Bit6Int(QX_INT(24 downto 19))); QY_MULT(31 downto 0) <= ("000000" & QY_MULT(31 downto 6)) + K1_ROM_Y2(Bit6Int(QY_INT(24 downto 19))); when Z7 => QX_MULT(31 downto 0) <= ("000000" & QX_MULT(31 downto 6)) + K1_ROM_X2(Bit6Int(QX_INT(30 downto 25))); QY_MULT(31 downto 0) <= ("000000" & QY_MULT(31 downto 6)) + K1_ROM_Y2(Bit6Int(QY_INT(30 downto 25)));

Aus dem Codeausschnitt erkennt man, dass ab Z5 auf die Register QX_INT und QY_INT zur Adressbildung zugegriffen wird. Die zusätzlichen 32-Bit-Variablen V_ROM_X1_OUT und V_ROM_Y1_OUT werden benötigt, um die jeweils höchstwertigen 24 Bits des Datenwortes des oberen, nur im Zustand Z4 adressierten ROMs auf den Addierer zu legen. Am Addierer stellen diese 24 Bits die niederwertigsten Bits dar und werden linksseitig um 6 Nullen ergänzt. Weiterhin wird in Z4 das Vorzeichen von QX_INT bzw. QY_INT in die höchstwertige Stelle von QX_MULT bzw. QY_MULT übernommen.

34

S OV H 23-Bit-Mantisse Addition-Bits

31 30 28 27 05 4

S Vorkommastellen Nachkommastellen

31 30 2425 0

S OV H 23-Bit-Mantisse AB

31 30 25 24 02 1

a)

b)

c) AOV

drei zusätzliche Overflow-Bits

0

32

S OV H 23-Bit-Mantisse Addition-Bits

31 30 28 27 05 4

S Vorkommastellen Nachkommastellen

31 30 2425 0

S OV H 23-Bit-Mantisse AB

31 30 25 24 02 1

a)

b)

c) AOV

drei zusätzliche Overflow-Bits

0

32

4.6 Zusätzliche Operationen für ln-Berechnung Da der Logarithmus über den Areatangens hyperbolicus bestimmt wird, ist hier keine Skalierung der CORDIC-Berechnung notwendig. Allerdings ist der artanh nach Gl. (2.40) noch mit dem Faktor 2 zu multiplizieren und hierzu nach Gl. (2.58) noch der Faktor E·ln2 zu addieren. All diese Operationen werden in den für die Konstantenmultiplikation vorgesehenen Zuständen Z4 bis Z7 ausgeführt. Das Produkt E·ln2 ist als ROM-Speicher (K3_ROM) der Organisation 64*32 im Package-File deklariert. In Z4 wird die Multiplikation mit 2 durch einen Linksshift realisiert. Zusätzlich wird der Referenzexponent QREF_EXP_XY in die Vorzeichen-Betrag-Darstellung umgewandelt. In Z5 werden die 6 niederwertigsten Bits des Referenzexponenten zur Adressierung von K3_ROM verwendet und das 32-Bit-Datenwort im 33-Bit-Register QE_LN_2[32:0] gespeichert. Die Breite dieses Registers orientiert sich an QX_MULT bzw. QY_MULT und spielt eine Rolle bei der Normalisierung der Mantisse innerhalb der Endstufe im Zustand Z8. In Z5 wird das höchstwertige Bit QE_LN_2(32) auf ’0’ gesetzt. Die übrigen Bits repräsentieren das adressierte ROM-Datenwort. Ist der Referenzexponent negativ, dann erfolgt noch eine entsprechende Komplementbildung, wie der nachfolgende Codeausschnitt zeigt: when Z5 => -- OP_CODE = X’04’ if QREF_EXP_XY(9) = '1' then QE_LN_2 <= ('0' & not (K3_ROM(Bit6Int(QREF_EXP_XY(5 downto 0))))) + 1; else QE_LN_2 <= '0' & K3_ROM(Bit6Int(QRE F_EXP_XY(5 downto 0))); end if ;

Um nun die Größen artanh (QZ_INT) und E·ln2 (QE_LN_2) addieren zu können, betrachten wir uns zunächst deren Datenformat nach Bild 4.7:

Bild 4.7: Datenformate für: a) QE_LN_2; b) QZ_INT; c) QZ_INT mit drei zusätzlichen OV-Bits Das Datenformat für QE_LN_2 weist 6 Vor- und 25 Nachkommastellen auf. Gemäß Bild 4.2 sind für QZ_INT jedoch nur drei Vorkommastellen vorgesehen, d.h. das Komma von QZ_INT muss um drei Stellen nach rechts verschoben werden. Da die Datenwortbreite von 32 Bit beibehalten werden soll, gehen somit die drei niederwertigsten Addition-Bits verloren. Die drei zusätzlichen Overflow-Bits (AOV in Bild 4.7) werden auf “000“ gesetzt, wenn QZ_INT(31) = ’0’ ist, andernfalls auf “111“. Diese Verschiebung des Dezimalpunktes erfolgt ebenfalls im Zustand Z5: when Z5 => -- OP_CODE = X’04’ -- Für die Addition von QE_LN_2 werden in QZ_INT drei zusätzliche Overflow-Bits -- eingefügt, wodurch die drei niederwe rtigsten Addition-Bits verloren gehen; if QZ_INT(31) = '1' then QZ_INT <= QZ_INT(31) & "111" & QZ_IN T(30 downto 3); else QZ_INT <= QZ_INT(31) & "000" & QZ_IN T(30 downto 3); end if ;

Im Zustand Z6 werden die beiden Register QE_LN_2 und QZ_INT addiert. Die Summe steht im Hinblick auf die Normalisierung der Mantisse (siehe Abschnitt 4.8) im 33-Bit-Register QE_LN_2:

35

when Z6 => -- OP_CODE = X’04’ QE_LN_2(31 downto 0) <= QZ_INT + QE_LN_2(31 downto 0);

Im Zustand Z7 wird schließlich das Zweierkomplement in die Vorzeichen-Betrag-Darstellung umgewandelt, wobei die beiden höchstwertigen Bits vertauscht werden (in Bild 4.7a durch einen Pfeil angedeutet). Damit steht das Vorzeichen jetzt, analog zu QX_MULT und QY_MULT, an Bitposition 32 und QE_LN_2(31) ist gelöscht (es werden ja nur 6 Vorkommastellen benötigt). Somit sind die 33-Bit-Datenformate für QX_MULT, QY_MULT und QE_LN_2 identisch: when Z7 => -- OP_CODE = X’04’ if QE_LN_2(31) = '1' then QE_LN_2 <= (QE_LN_2(31) & '0' & ( not (QE_LN_2(30 downto 0)))) + 1; else QE_LN_2 <= QE_LN_2(31) & '0' & QE_L N_2(30 downto 0); end if ;

4.7 Zusätzliche Operationen für sqrt-Berechnung Neben der Skalierung des in QX_INT stehenden Wertes sind bei der Berechnung der Quadratwurzel aufgrund der Erweiterung des Argument-Bereiches noch weitere Operationen notwendig, die parallel zur Konstantenmultiplikation ausgeführt werden. Zu Beginn wurde in der Vorstufe bereits geprüft, ob der Exponent gerade ist oder nicht (siehe Abschnitt 4.2) und das Resultat im Flipflop QE_MOD_2 gespeichert. Im Falle eines ungeraden Exponenten ist eine Teilung der Mantisse nach Gl. (2.60b) ebenfalls in der Vorstufe vorgesehen (siehe Zustand Z1). Nun wird in den Zuständen Z5 und Z6 der Referenzexponent angepasst. Ist QE_MOD_2 = ’1’, dann wird in Z5 gemäß Gl. (2.60b) QREF_EXP_XY um 1 inkrementiert: when Z5 => -- OP_CODE = X“05“ if QE_MOD_2 = '1' then QREF_EXP_XY <= QREF_EXP_XY + 1; end if ;

Im Zustand Z6 erfolgt schließlich die Teilung des Referenzexponenten mittels Shiftoperation. Da der 10-stellige Referenzexponent im Zweierkomplement vorliegt, muss zwischen einem positiven und negativen Vorzeichen unterschieden werden: when Z6 => -- OP_CODE = X“05“ if QREF_EXP_XY(9) = '1' then -- Exponent ist negativ => von link s eine ’1’ nachziehen QREF_EXP_XY <= QREF_EXP_XY(9) & '1' & QREF_EXP_XY(8 downto 1); else QREF_EXP_XY <= QREF_EXP_XY(9) & '0' & QREF_EXP_XY(8 downto 1); end if ;

4.8 Endstufe Die Endstufe, d.h. die Umwandlung des internen 33-Bit-Festkommaformats (siehe Bild 4.8) in das externe 32-Bit-Gleitkommaformat, entspricht dem Automatenzustand Z8. Hier sind für die Korrektur des Referenzexponenten in Abhängigkeit vom jeweiligen Datenpfad die drei 10-Bit-Variablen VREF_EXP_X, VREF_EXP_Y und VREF_EXP_Z vorgesehen. Die Notwendigkeit einer Korrektur ergibt sich zum einen aus der vorzeichenlosen Darstellung als Charakteristik, zum anderen aus der Kommaverschiebung bei der Normalisierung der Mantisse. Die Anzahl der Schiebeoperationen zur Mantissen-Normalisierung wird für jeden Datenpfad kombinatorisch über den Block COUNT_SH_OP nach Bild 4.8 ermittelt. Das 33-Bit-Datenformat orientiert sich hierbei an den Registern QX_MULT und QY_MULT zur Konstantenmultiplikation und setzt sich entsprechend zusammen aus einem VZ-Bit, 3 OV-Bits, einem Hidden-Bit, einer 23-Bit-Mantisse und 5 Addition-Bits. Eine Verschiebung des Kommas nach links ist im Exponenten positiv, eine nach rechts entsprechend negativ (10-Bit-Zweierkomplement) zu bewerten, siehe Bild 4.8. Der VHDL-Quellcode für diese kombinatorische Logik findet sich in dem File COUNT_SH_OP.vhd.

36

S OV H 23-Bit-Mantisse Addition-Bits

32 31 28 27 05 4

1

1

X

1

1

X X X X X X X X X X X X X X X X X X X X X X X X X

1

1

0

0

1

1

0

0 0 0

XX

X

X

X

X

X

X

X

X

X

0

0 0

0

0 00

0 0 0

X X X X X X X X X X X X X X X X X X X X X

X X X X X X X X X X X X X X X X X X X X X

X X X X X X X X X X X X X X X X X X X X X

X X X X X X X X X X X X X X X X X X X X X

X X X X X X X X X X X X X X X X X X X X X

0

0 0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0

0 0 0

0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0

-3

-2

-1

0

1

2

22

23

S OV H 23-Bit-Mantisse Addition-Bits

32 31 28 27 05 4

1

1

X

1

1

X X X X X X X X X X X X X X X X X X X X X X X X X

1

1

0

0

1

1

0

0 0 0

XX

X

X

X

X

X

X

X

X

X

0

0 0

0

0 00

0 0 0

X X X X X X X X X X X X X X X X X X X X X

X X X X X X X X X X X X X X X X X X X X X

X X X X X X X X X X X X X X X X X X X X X

X X X X X X X X X X X X X X X X X X X X X

X X X X X X X X X X X X X X X X X X X X X

0

0 0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0

0 0 0

0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0

-3

-2

-1

0

1

2

22

23

QX_MULT[32:0]

COUNT_SH_OPCountX_mult[9:0]

VREF_EXP_X[9:0]

QY_MULT[32:0]

COUNT_SH_OPCountY_mult[9:0]

QE_LN_2[32:0]

COUNT_SH_OPCountZ_add[9:0]

VREF_EXP_Y[9:0] VREF_EXP_Z[9:0]

Korrektur des Referenzexponenten

SH

FT

32LR

Man

Zou

t[23:

0]

QE

_LN

_2[3

2:0]

SH

FT

32LR

Man

You

t[23:

0]

QY

_MU

LT[3

2:0]

SH

FT

32LR

Man

Xou

t[23:

0]

QX

_MU

LT[3

2:0]

Man

Xou

t[22

:0]

Man

You

t[22:

0]M

anZ

out[2

2:0]

VREF_EXP_X[7:0] VREF_EXP_Y[7:0] VREF_EXP_Z[7:0]

QX_MULT[32]

QY_MULT[32]

QE_LN_2[32]

QQUOTIENT[9:0]

OP_CODE[7:0]

FALLUNTERSCHEIDUNG[1:0]

X“03“: exp-Funktion

QXFP_OUT[31:0] speichern;QYZFP_OUT[31:0] löschen;

X“04“: ln-Funktion

QXFP_OUT[31:0] löschen;QYZFP_OUT[31:0] speichern;

X“05“: sqrt-Funktion

QXFP_OUT[31:0] speichern;QYZFP_OUT[31:0] löschen;

X“01“: cos und sin-Funktion

(Fallunterscheidung nach Tabelle 4.2)QXFP_OUT[31:0] speichern;QYZFP_OUT[31:0] speichern;

QXFP_OUT[31:0] QYZFP_OUT[31:0]

QX_MULT[32:0]

COUNT_SH_OPCountX_mult[9:0]

VREF_EXP_X[9:0]

QY_MULT[32:0]

COUNT_SH_OPCountY_mult[9:0]

QE_LN_2[32:0]

COUNT_SH_OPCountZ_add[9:0]

VREF_EXP_Y[9:0] VREF_EXP_Z[9:0]

Korrektur des Referenzexponenten

SH

FT

32LR

Man

Zou

t[23:

0]

QE

_LN

_2[3

2:0]

SH

FT

32LR

Man

You

t[23:

0]

QY

_MU

LT[3

2:0]

SH

FT

32LR

Man

Xou

t[23:

0]

QX

_MU

LT[3

2:0]

Man

Xou

t[22

:0]

Man

You

t[22:

0]M

anZ

out[2

2:0]

VREF_EXP_X[7:0] VREF_EXP_Y[7:0] VREF_EXP_Z[7:0]

QX_MULT[32]

QY_MULT[32]

QE_LN_2[32]

QQUOTIENT[9:0]

OP_CODE[7:0]

FALLUNTERSCHEIDUNG[1:0]

X“03“: exp-Funktion

QXFP_OUT[31:0] speichern;QYZFP_OUT[31:0] löschen;

X“04“: ln-Funktion

QXFP_OUT[31:0] löschen;QYZFP_OUT[31:0] speichern;

X“05“: sqrt-Funktion

QXFP_OUT[31:0] speichern;QYZFP_OUT[31:0] löschen;

X“01“: cos und sin-Funktion

(Fallunterscheidung nach Tabelle 4.2)QXFP_OUT[31:0] speichern;QYZFP_OUT[31:0] speichern;

QXFP_OUT[31:0] QYZFP_OUT[31:0]

Bild 4.8: Anzahl der Schiebeoperationen zur Mantissen-Normalisierung Durch nachfolgend angegebene Komponenteninstanziierung im Top-Level-Modul steht die (nur in Z8 in Abhängigkeit vom OP-Code gültige) Anzahl der Schiebeoperationen in den 10-Bit-Signalvektoren CountX_mult, CountY_mult sowie CountZ_add. Mit diesen Werten können dann durch die Blöcke SHFT32LR die normalisierten 24-Bit-Mantissen ManXout, ManYout und ManZout bestimmt werden: count_sh_op_x: COUNT_SH_OP port map (QX_MULT,CountX_mult); count_sh_op_y: COUNT_SH_OP port map (QY_MULT,CountY_mult); count_sh_op_z: COUNT_SH_OP port map (QE_LN_2,CountZ_add); shftX_out: SHFT32LR port map (QX_MULT,CountX_mult,ManXout); shftY_out: SHFT32LR port map (QY_MULT,CountY_mult,ManYout); shftZ_out: SHFT32LR port map (QE_LN_2,CountZ_add,ManZout);

Bild 4.9: Endstufe

37

Jetzt können, ausgehend vom vorliegenden OP-Code, im Zustand Z8 die in den Registern QXFP_OUT[31:0] bzw. QYZFP_OUT[31:0] zu speichernden 32-Bit-Gleitkommazahlen aus den Vorzeichenbits, den korrigierten Referenzexponenten sowie den normalisierten Mantissen zusammengesetzt werden, siehe Bild 4.9 bzw. Anhang A10. Man erkennt, dass das Hidden-Bit der normalisierten Mantisse sowie die beiden höchstwertigen Bits des Referenzexponenten (Darstellung als Charakteristik) jeweils unberücksichtigt bleiben. Bei der Berechnung der cos- und sin-Funktion ist noch die Fallunterscheidung nach Tabelle 4.2 vorzunehmen. 4.9 Simulation Zur rein logischen Verifikation des VHDL-Quellcodes wird der Simulator ModelSim XE Starter der Fa. Mentor Graphics verwendet. Die Umrechnung zwischen Dezimalzahlen und 32-Bit-Gleitkommazahlen erfolgt über ein Java-Applet unter der URL http://www.h-schmidt.net/FloatApplet/IEEE754de.html. Die Stimulidaten stehen in einer do-Datei, die für die Bilder 4.10 bzw. 4.11 wie folgt lautet: # Kommandodatei für Bild 4.10; Berechnung von exp(4 .85) restart -f radix hex force reset 0 0, 1 70ns force start_cordic 0 0, 1 175ns, 0 320ns force clk 0 0, 1 50ns -repeat 100ns force op_code x"03" force xfp x"3F800000" force yfp x"00000000" force zfp x"409B3333" run 3.6us # Kommandodatei für Bild 4.11; Berechnung von cos(3 .5) und sin(3.5) sowie ln(112.26) restart -f radix hex force reset 0 0, 1 70ns force start_cordic 0 0, 1 175ns, 0 320ns, 1 3.56us, 0 3.72us force clk 0 0, 1 50ns -repeat 100ns force op_code x"01" 0, x"04" 3.57us force xfp x"3F800000" 0, x"42E0851F" 3.53us force yfp x"00000000" 0, x"42E0851F" 3.53us force zfp x"40600000" 0, x"00000000" 3.53us run 7.2us In der Kommandodatei für Bild 4.10 werden XFP auf 1.0d, YFP auf 0.0 und ZFP auf 4.85d gesetzt. Mittels OP_CODE = X“03“ wird folglich die Exponentialfunktion berechnet. Das Ergebnis steht im Register QXFP_OUT und entspricht umgewandelt folgendem dezimalem Wert (zum Vergleich ist der mit einem Taschenrechner ermittelte Wert ebenfalls angegeben): QXFP_OUT[31:0] = 42FF7B0Fh → 127.74035d

exp(4.85) = 127.7403898… (Tasch enrechner)

In Bild 4.11 sind zwei CORDIC-Berechnungen zu sehen, zunächst für OP_CODE = X“01“ (cos- und sin-Funktion) und anschließend für OP_CODE = X“04“ (ln-Funktion). Das Argument für die cos- und sin-Funktion ist 3.5d, der Numerus für die ln-Funktion ist 112.26. Mit diesen Eingabedaten erhält man folgende Ergebnisse: 1. CORDIC-Berechnung: QXFP_OUT[31:0] = BF6FBBA0h → -0.9364567d cos(3. 5) = -0.936456687… (Taschenrechner)

QYZFP_OUT[31:0] = BEB399D 8h → -0.3507831d sin(3. 5) = -0.350783227… (Taschenrechner) 2. CORDIC-Berechnung: QYZFP_OUT[31:0] = 409710EFh → 4.720817d ln(112.2 6) = 4.72081760… (Taschenrechner)

In Tabelle 4.8 sind für jeden OP-Code einige Testdaten mit den zugehörigen Ergebnissen sowie einem Vergleich zur jeweils exakten Lösung zusammengefasst. Man erkennt, dass die 24 CORDIC-Iterationen auf ein etwa 6 Stellen genaues Ergebnis führen.

38

Bild 4.10: Berechnung von exp(4.85)

Bild 4.11: Berechnung von cos(3.5) und sin(3.5) sowie ln(112.26)

39

OP_CODE z bzw. r XFP YFP ZFP QXFP_OUT dezimal Taschenrechner QYZFP_OUT dezimal Taschenrechner

01 0.000 3F800000 00000000 00000000 3F7FFFFF 0.99999994 1.000000000 B3800000 -5.9604645E-8 0.000000000 01 0.080 3F800000 00000000 3DA3D70A 3F7F2E65 0.9968017 0.996801706 3DA3AA47 0.07991462 0.079914693 01 0.670 3F800000 00000000 3F2B851F 3F48A888 0.7838216 0.783821665 3F1EF8F0 0.620986 0.620985987 01 1.520 3F800000 00000000 3FC28F5C 3D4FF906 0.050774597 0.050774484 3F7FAB77 0.9987101 0.998710143 01 2.270 3F800000 00000000 401147AE BF24C385 -0.6436084 -0.643608418 3F43EE4C 0.7653549 0.765354952 01 3.055 3F800000 00000000 4043851F BF7F0A72 -0.99625313 -0.996253198 3DB11EC4 0.08648446 0.086484478 01 4.937 3F800000 00000000 409DFBE7 3E641281 0.22272684 0.222727169 BF7991CB -0.9748809 -0.974880817 01 5.377 3F800000 00000000 40AC1062 3F1DE384 0.61675286 0.616752990 BF49831C -0.7871568 -0.787156749 01 5.900 3F800000 00000000 40BCCCCD 3F6D6F3A 0.92747843 0.927478430 BEBF6CC1 -0.3738766 -0.373876664 01 6.180 3F800000 00000000 40C5C28F 3F7EA36B 0.99468106 0.994681117 BDD2F2FE -0.10300253 -0.103002298 03 0.000 3F800000 00000000 00000000 3F800001 1.0000001 1.000000000 00000000 - - 03 0.430 3F800000 00000000 3EDC28F6 3FC4C4D9 1.5372573 1.537257523 00000000 - - 03 1.550 3F800000 00000000 3FC66666 4096C45E 4.7114706 4.711470182 00000000 - - 03 3.905 3F800000 00000000 4079EB85 424699AF 49.65008 49.65007991 00000000 - - 03 8.330 3F800000 00000000 410547AE 45819356 4146.417 4146.417552 00000000 - - 03 15.950 3F800000 00000000 417F3333 4B00FA79 8452729.0 8452729.796 00000000 - - 03 -0.005 3F800000 00000000 BBA3D70A 3F7EB925 0.9950126 0.995012479 00000000 - - 03 -0.520 3F800000 00000000 BF051EB8 3F183280 0.59452057 0.594520547 00000000 - - 03 -1.107 3F800000 00000000 BF8DB22D 3EA93DBA 0.33054906 0.330549122 00000000 - - 03 -8.336 3F800000 00000000 C1056042 397B5FD3 2.3972923E-4 2.39729345E-4 00000000 - -

04 0.00005 3851B717 3851B717 00000000 00000000 - - C11E74AE -9.903486 -9.903487552 04 0.008 3C03126F 3C03126F 00000000 00000000 - - C09A818A -4.828313 -4.828313737 04 1.825 3FE9999A 3FE9999A 00000000 00000000 - - 3F1A011F 0.6015796 0.601579987 04 2.71828 402DF84D 402DF84D 00000000 00000000 - - 3F7FFFF3 0.9999992 0.999999327 04 22.6 41B4CCCD 41B4CCCD 00000000 00000000 - - 40478C7D 3.1179497 3.117949906 04 500.8 43FA6666 43FA6666 00000000 00000000 - - 40C6EB29 6.216206 6.216206819 04 4200.5 45834400 45834400 00000000 00000000 - - 41057CC1 8.3429575 8.342958844 04 22000 46ABE000 46ABE000 00000000 00000000 - - 411FFB12 9.998796 9.998797732 04 522560 48FF2800 48FF2800 00000000 00000000 - - 4152A9F5 13.166493 13.16649508 04 1.52E+6 49B98C00 49B98C00 00000000 00000000 - - 4163BF5C 14.234219 14.23422089 05 0.00005 3851B717 3851B717 00000000 3BE7B469 0.007071067 0.007071067 00000000 - - 05 0.008 3C03126F 3C03126F 00000000 3DB72DBE 0.089442715 0.089442719 00000000 - - 05 1.825 3FE9999A 3FE9999A 00000000 3FACEB21 1.3509256 1.350925608 00000000 - - 05 2.71828 402DF84D 402DF84D 00000000 3FD30947 1.6487206 1.648720716 00000000 - - 05 22.6 41B4CCCD 41B4CCCD 00000000 40982052 4.7539454 4.753945729 00000000 - - 05 500.8 43FA6666 43FA6666 00000000 41B3074A 22.37856 22.37856116 00000000 - - 05 4200.5 45834400 45834400 00000000 42819F5D 64.81126 64.81126445 00000000 - - 05 22000 46ABE000 46ABE000 00000000 431452EF 148.32396 148.3239697 00000000 - - 05 522560 48FF2800 48FF2800 00000000 4434B884 722.88306 722.8831164 00000000 - - 05 1.52E+6 49B98C00 49B98C00 00000000 449A1C3F 1232.8827 1232.882800 00000000 - -

Tabelle 4.8: Zusammenfassung einiger Testdaten

40

Anhang A1 C-Quellcode zur Berechnung von sin() und cos( ) // Berechnung von sin und cos #include <stdio.h> #include <conio.h> #include <math.h> int main (void) { int n = 1; double g = 1.64676025812, delta = 1; double tmpx, tmpy, VZ; double x = 1, y = 0, z = 0.6; // Initialwert für z entspricht phi do { if (z < 0) VZ = -1; else VZ = 1; tmpx = x - VZ*delta*y; tmpy = y + VZ*delta*x; x = tmpx; y = tmpy; z = z - VZ*atan(delta); delta = delta / 2; n++; } while (n <= 15); x = x / g; y = y / g; printf("\n\n x = %1.8lf", x); printf("\n\n y = %1.8lf\n\n", y); getch(); return 0; } A2 C-Quellcode zur Berechnung von arctan() // Berechnung von arctan #include <stdio.h> #include <conio.h> #include <math.h> int main (void) { int n = 1; double g = 1.64676025812, delta = 1; double tmpx, tmpy, VZ; double x = 1.8, y = 0.7, z = 0; do { if (y < 0) VZ = 1; else VZ = -1; tmpx = x - VZ*delta*y; tmpy = y + VZ*delta*x; x = tmpx; y = tmpy; z = z - VZ*atan(delta); delta = delta / 2; n++; } while (n <= 15); x = x / g; printf("\n\n\n x = %1.8lf", x); printf("\n z = %1.8lf", z); getch(); return 0; }

41

A3 C-Quellcode zur Berechnung von arcsin() // Berechnung von arcsin #include <stdio.h> #include <conio.h> #include <math.h> int main (void) { int n = 1; double g = 1.64676025812, delta = 1; double tmpx, tmpy, VZ; double x = 1 / g, y = 0, z = 0, arg = 0.7; do { if (arg - y < 0) VZ = -1; else VZ = 1; tmpx = x - VZ*delta*y; tmpy = y + VZ*delta*x; x = tmpx; y = tmpy; z = z + VZ*atan(delta); delta = delta / 2; n++; } while (n <= 15); printf("\n\n\n x = %1.8lf", x); printf("\n y = %1.8lf", y); printf("\n\n z = %1.8lf\n\n", z); getch(); return 0; } A4 C-Quellcode zur Berechnung von sinh() und cos h() // Berechnung von sinh und cosh // einfache Shiftfolge => Konvergenzkriterium nicht erfüllt #include <stdio.h> #include <conio.h> #include <math.h> int main (void) { int n = 1; double g = 0.82978162013, delta = 0.5; double tmpx, tmpy, VZ; double x = 1, y = 0, z = 0.8; do { if (z < 0) VZ = -1; else VZ = 1; tmpx = x + VZ*delta*y; tmpy = y + VZ*delta*x; x = tmpx; y = tmpy; z = z - VZ*0.5*log((1+delta)/(1-delta)); delta = delta / 2; n++; } while (n <= 15); x = x / g; y = y / g; printf("\n\n x = %1.8lf", x); printf("\n\n y = %1.8lf\n\n", y); getch(); return 0; }

42

A5 C-Quellcode zur Berechnung von sinh() und cos h() // Berechnung von sinh und cosh // komplexere Shiftfolge => Konvergenzkriterium erf üllt #include <stdio.h> #include <conio.h> #include <math.h> int main (void) { int n = 1, k = 4; double g = 0.82815936096, delta = 0.5; double tmpx, tmpy, VZ; double x = 1, y = 0, z = 0.01; do { if (z < 0) VZ = -1; else VZ = 1; tmpx = x + VZ*delta*y; tmpy = y + VZ*delta*x; x = tmpx; y = tmpy; z = z - VZ*0.5*log((1+delta)/(1-delta)); if (n == k) { k = 3*k + 1; } else { delta = delta / 2; n++; } } while (n <= 32); x = x / g; y = y / g; printf("\n\n x = %1.10lf", x); printf("\n\n y = % 1.10lf\n\n", y); getch(); return 0; } A6 C-Quellcode zur Berechnung von artanh() // Berechnung von artanh(); Konvergenzkriterium erf üllt #include <stdio.h> #include <conio.h> #include <math.h> int main (void) { int n = 1, k = 4; double g = 0.82815936096, delta = 0.5; double tmpx, tmpy, VZ; double x = 3.9, y = 2.5, z = 0; do { if (y < 0) VZ = 1; else VZ = -1; tmpx = x + VZ*delta*y; tmpy = y + VZ*delta*x; x = tmpx; y = tmpy; z = z - VZ*0.5*log((1+delta)/(1-delta)); if (n == k) { k = 3*k + 1; } else { delta = delta / 2; n++; } } while (n <= 15); x = x / g; printf("\n\n x = %1.8lf", x); printf("\n\n z = %1 .8lf\n\n", z); getch(); return 0; }

43

A7 C-Quellcode zur Multiplikation im linearen Mo dus // Multiplikation im linearen Modus (Rotating) #include <stdio.h> #include <conio.h> #include <math.h> int main (void) { int n = 1; double delta = 1, VZ; double x = 3.5, y = 0, z = 1.9; do { if (z < 0) VZ = -1; else VZ = 1; y = y + VZ*delta*x; z = z - VZ*delta; delta = delta / 2; n++; } while (n <= 24); printf("\n\n y = %1.8lf\n\n", y); getch(); return 0; }

A8 C-Quellcode zur Division im linearen Modus // Division im linearen Modus (Vectoring) #include <stdio.h> #include <conio.h> #include <math.h> int main (void) { int n = 1; double delta = 1, VZ; double x = 4.1, y = -5.8, z = 0; do { if (y < 0) VZ = 1; else VZ = -1; y = y + VZ*delta*x; z = z - VZ*delta; delta = delta / 2; n++; } while (n <= 24); printf("\n\n z = %1.8lf\n\n", z); getch(); return 0; }

44

A9 C-Quellcode zur Berechnung von exp() // Berechnung von exp() über cosh und sinh // komplexere Shiftfolge => Konvergenzkriterium für cosh und sinh erfüllt // Berücksichtigung unzulässiger Argumente #include <stdio.h> #include <conio.h> #include <math.h> int main (void) { int n = 1, k = 4, m; double g = 0.82815936096, de lta = 0.5; double tmpx, tmpy, VZ; double z, x = 1, y = 0, arg = -5.7; double ln_2 = 0.69314718056; double exp_fkt; m = (int) (arg / ln_2); z = arg - m*ln_2; do { if (z < 0) VZ = -1; else VZ = 1; tmpx = x + VZ*delta*y; tmpy = y + VZ*delta*x; x = tmpx; y = tmpy; z = z - VZ*0.5*log((1+delta)/(1-delta)); if (n == k) { k = 3*k + 1; } else { delta = delta / 2; n++; } } while (n <= 32); x = x / g; y = y / g; // Berechnung der Exponentialfunktion exp_fkt = pow(2,m)*(x + y); printf("\n\n exp_fkt = %1.10lf", exp_fkt); getch( ); return 0; }

45

A10 VHDL-Quellcode für CORDIC-Algorithmus (Stand : 19.06.2008) -- Top_Level_Modul.vhd library IEEE; use IEEE.STD_LOGIC_1164. ALL; use IEEE.STD_LOGIC_ARITH. ALL; use IEEE.STD_LOGIC_UNSIGNED.ALL; use work.PACK_CORDIC. ALL; entity CORDIC_VERSION_1 is port (CLK,RESET,START_CORDIC: in std_logic; XFP,YFP,ZFP: in std_logic_vector(31 downto 0); OP_CODE: in std_logic_vector(7 downto 0); READY_CORDIC: out std_logic; XFP_OUT,YZFP_OUT: out std_logic_vector(31 downto 0)); end CORDIC_VERSION_1; architecture CORDIC_VERSION_1_A of CORDIC_VERSION_1 is -- Komponenten-Deklarationen component SHFT24LR port (E: in std_logic_vector(23 downto 0); COUNT: in std_logic_vector(7 downto 0); SHIFT_DIRECTION: in std_logic; OV_BITS: out std_logic_vector(1 downto 0); A: out std_logic_vector(23 downto 0)); end component SHFT24LR; component CONV_FP_TC port (MANTISSE: in std_logic_vector(23 downto 0); VZ: in std_logic; OV_BITS: in std_logic_vector(1 downto 0); S_INT: out std_logic_vector(31 downto 0)); end component CONV_FP_TC; component COUNT_SH_OP port (E: in std_logic_vector(32 downto 0); COUNT: out std_logic_vector(9 downto 0)); end component COUNT_SH_OP; component SHFT32LR port (E: in std_logic_vector(32 downto 0); COUNT: in std_logic_vector(9 downto 0); A: out std_logic_vector (23 downto 0)); end component SHFT32LR; component EXTENDING_DOMAIN port (EXPONENT: in std_logic_vector(7 downto 0); MANTISSE: in std_logic_vector(23 downto 0); MANTISSE_SHIFT: in std_logic_vector(23 downto 0); QUOTIENT: out std_logic_vector(4 downto 0); REST: out std_logic_vector (23 downto 0)); end component EXTENDING_DOMAIN;

46

-- Signal-Deklarationen type ZUSTAENDE is (Z0,Z1,Z2,Z3,Z4,Z5,Z6,Z7,Z8); --Aufzählungstyp signal ZUSTAND,FOLGE_Z: ZUSTAENDE; --Inte rprozesskommunikation signal ManX,ManY,ManZ,ManXshft,ManYshft,ManZshft: std_log ic_vector(23 downto 0); signal ManXout,ManYout,ManZout: std_logic_vector(23 downto 0); signal XgtY, YgtX,Zgteq127: std_logic; signal CountX,CountY,CountZ: std_logic_vector(7 downto 0); signal OV_BITS_X,OV_BITS_Y,OV_BITS_Z: std_logic_vector(1 downto 0); signal X_INT,Y_INT,Z_INT: std_logic_vector(31 downto 0); signal QX_INT,QY_INT,QZ_INT: std_logic_vector(31 downto 0); signal QREF_EXP_XY: std_logic_vector(9 downto 0); signal QX_MULT,QY_MULT,QE_LN_2: std_logic_vector(32 downto 0); signal CountX_mult,CountY_mult,CountZ_add: std_logic_vect or(9 downto 0); signal QS: std_logic_vector(4 downto 0); signal oi,m: integer range -1 to +1; signal QUOTIENT: std_logic_vector(4 downto 0); signal QQUOTIENT: std_logic_vector(9 downto 0); signal REST: std_logic_vector(23 downto 0); signal FALLUNTERSCHEIDUNG: std_logic_vector(1 downto 0); signal QXFP_OUT,QYZFP_OUT: std_logic_vector(31 downto 0); signal QE_MOD_2: std_logic; begin -- Kombinatorik ManX <= '1' & XFP(22 downto 0); ManY <= '1' & YFP(22 downto 0); ManZ <= '1' & ZFP(22 downto 0); XgtY <= '1' when XFP(30 downto 23) > YFP(30 downto 23) else '0'; YgtX <= '1' when YFP(30 downto 23) > XFP(30 downto 23) else '0'; Zgteq127 <= '1' when ZFP(30 downto 23) >= 127 else '0'; CountX <= YFP(30 downto 23) - XFP(30 downto 23) when YgtX = '1' else "00000000"; CountY <= XFP(30 downto 23) - YFP(30 downto 23) when XgtY = '1' else "00000000"; CountZ <= ZFP(30 downto 23) - 127 when Zgteq127 = '1' else (127 - ZFP(30 downto 23)); shftX: SHFT24LR port map (ManX,CountX,'0',OV_BITS_X,ManXshft); shftY: SHFT24LR port map (ManY,CountY,'0',OV_BITS_Y,ManYshft); shftZ: SHFT24LR port map (ManZ,CountZ,Zgteq127,OV_BITS_Z,ManZshft); convx_int: CONV_FP_TC port map(ManXshft,XFP(31),OV_BITS_X,X_INT); convy_int: CONV_FP_TC port map(ManYshft,YFP(31),OV_BITS_Y,Y_INT); convz_int: CONV_FP_TC port map(ManZshft,ZFP(31),OV_BITS_Z,Z_INT); count_sh_op_x: COUNT_SH_OP port map(QX_MULT,CountX_mult); count_sh_op_y: COUNT_SH_OP port map(QY_MULT,CountY_mult); count_sh_op_z: COUNT_SH_OP port map(QE_LN_2,CountZ_add); shftX_out: SHFT32LR port map(QX_MULT,CountX_mult,ManXout); shftY_out: SHFT32LR port map(QY_MULT,CountY_mult,ManYout); shftZ_out: SHFT32LR port map(QE_LN_2,CountZ_add,ManZout); ext_domain: EXTENDING_DOMAIN port map(ZFP(30 downto 23),ManZ,ManZshft,QUOTIENT,REST);

47

Z_SPEICHER: process (CLK,RESET) begin if RESET = '0' then ZUSTAND <= Z0; QS <= "00000"; elsif CLK = '1' and CLK'event then ZUSTAND <= FOLGE_Z; case ZUSTAND is when Z0 => QS <= "00000"; when Z3 => QS <= QS + 1; when others => end case ; end if ; end process Z_SPEICHER; F_BEST: process (START_CORDIC,ZUSTAND,QS) begin case ZUSTAND is when Z0 => if START_CORDIC = '1' then FOLGE_Z <= Z1; else FOLGE_Z <= Z0; end if ; when Z1 => FOLGE_Z <= Z2; when Z2 => FOLGE_Z <= Z3; when Z3 => if QS = "10111" then FOLGE_Z <= Z4; -- Abbruch d es CORDIC-Verfahrens nach 24 Iterationen end if ; when Z4 => FOLGE_Z <= Z5; when Z5 => FOLGE_Z <= Z6; when Z6 => FOLGE_Z <= Z7; when Z7 => FOLGE_Z <= Z8; when Z8 => FOLGE_Z <= Z0; end case ; end process F_BEST; A_BEST: process (ZUSTAND,QS) variable N,I,J: integer; variable VX_SHIFT,VY_SHIFT: std_logic_vector(31 downto 0); variable VX_INT,VY_INT,VZ_INT: std_logic_vector(31 downto 0); variable VREF_EXP_X,VREF_EXP_Y,VREF_EXP_Z: std_logic_vector (9 downto 0); variable V_ROM_X1_OUT,V_ROM_Y1_OUT: std_logic_vector(31 downto 0); variable VE_MOD_2: std_logic; begin case ZUSTAND is when Z0 => READY_CORDIC <= '1'; when Z1 => READY_CORDIC <= '0';

48

QQUOTIENT <= "0000000000"; -- Initialisierung QX_INT <= X_INT; QY_INT <= Y_INT; QZ_INT <= Z_INT; case OP_CODE is when X"01" => m <= 1; -- Falluntersc heidung für trigonometrische Funktionen cos und sin -- Argument-Be reich von 0 bis 2pi if Z_INT <= "00011001001000011111101101010100" then FALLUNTERSC HEIDUNG <= "00"; elsif (Z_INT > "00011001001000011111101101010100" and Z_INT <= "00110010010000111111011010101000") then FALLUNTERSC HEIDUNG <= "01"; QZ_INT <= Z _INT - "00011001001000011111101101010100"; elsif (Z_INT > "00110010010000111111011010101000" and Z_INT <= "01001011011001011111000111111100") then FALLUNTERSC HEIDUNG <= "10"; QZ_INT <= Z _INT - "00110010010000111111011010101000"; else -- Z_INT li egt zwischen 3pi/2 und 2pi FALLUNTERSC HEIDUNG <= "11"; QZ_INT <= Z _INT - "01001011011001011111000111111100"; end if ; when X"03" => m <= -1; -- Erweiterung des zulässigen Argument-Bereiches für exp-Funktion if ZFP(31) = '0' then QQUOTIENT < = "00000" & QUOTIENT; QZ_INT <= Z FP(31) & "00" & REST & "00000"; else QQUOTIENT < = ("11111" & not (QUOTIENT)) + 1; QZ_INT <= ( ZFP(31) & "11" & not (REST) & "11111") + 1; end if ; when X"04" => m <= -1; QX_INT <= X_IN T + X"10000000"; -- Addition von 1d QY_INT <= Y_IN T - X"10000000"; -- Subtraktion von 1d when others => m <= -1; VE_MOD_2 := E_ MOD_2(Bit8Int(XFP(30 downto 23))); QE_MOD_2 <= VE _MOD_2; -- falls Expon ent ungerade, dann wird Mantisse M geteilt => Besti mmung des Radikanden r -- es wird dav on ausgegangen, dass der Radikand positiv ist if VE_MOD_2 = '1' then QX_INT <= X _INT(31) & '0' & X_INT(30 downto 1); QY_INT <= Y _INT(31) & '0' & Y_INT(30 downto 1); end if ; end case ; when Z2 => -- Bestimmung des Referenzexponenten für X /Y if (XgtY = '0' and YgtX = '0') then QREF_EXP_XY <= ("00" & XFP(30 downto 23)) - 127; elsif XgtY = '1' then QREF_EXP_XY <= ("00" & XFP(30 downto 23)) - 127; else QREF_EXP_XY <= ("00" & YFP(30 downto 23)) - 127; end if ;

49

-- Initialisierung von oi und Be stimmung von M für Wurzelberechnung case OP_CODE is when X"01" => -- zirkularer Modus, Rotating if QZ_INT(31) = '0' then oi <= 1; - - z0 >= 0 => oi ist positiv else oi <= -1; - - z0 < 0 => oi ist negativ end if ; when X"03" => -- hyperbolischer Modus, Rotating if QZ_INT(31) = '0' then oi <= 1; - - z0 >= 0 => oi ist positiv else oi <= -1; - - z0 < 0 => oi ist negativ end if ; when X"04" => -- hyperbolischer Modus, Vectoring, Ber echnung ln if QY_INT(31) = '0' then oi <= -1; - - y0 >= 0 => oi ist negativ else oi <= 1; - - y0 < 0 => oi ist positiv end if ; when others => -- hyperbolischer Modus, Vectoring, Berechnung sqrt if QY_INT(31) = '0' then oi <= -1; - - y0 >= 0 => oi ist negativ else oi <= 1; - - y0 < 0 => oi ist positiv end if ; QX_INT <= QX_I NT + X"04000000"; -- Addition von 0.25d QY_INT <= QY_I NT - X"04000000"; -- Subtraktion von 0.25d end case ; when Z3 => -- Beginn der Iteration if m = 1 then -- einfache Shiftfolge für zi rkularen Modus N:=Bit5Int(QS); -- Rechtsshi ft um QS Stellen else -- komplexere Shiftfolge für hyperbolischen Modus if QS < "00100" then N:=Bit5Int(QS) + 1; elsif QS < "01110" then N:=Bit5Int(QS); else N:=Bit5Int(QS) - 1; end if ; end if ; for I in 0 to 31 loop if I >= N then VX_SHIFT(I-N):=QX_INT(I); VY_SHIFT(I-N):=QY_INT(I); else if QX_INT(31) = '0' then -- Zahl ist positiv => Null einschieben VX_SHIFT(32-N+I):='0';

50

else -- Zahl ist negativ => Eins einschieben VX_SHIFT(32-N+I):='1'; end if ; if QY_INT(31) = '0' then -- Zahl ist positiv => Null einschieben VY_SHIFT(32-N+I):='0'; else -- Zahl ist negativ => Eins einschieben VY_SHIFT(32-N+I):='1'; end if ; end if ; end loop ; if m = 1 then if oi = 1 then VX_INT := QX_INT - VY_SHIF T; VY_INT := QY_INT + VX_SHIF T; VZ_INT := QZ_INT - ROM_arc tan(N); else VX_INT := QX_INT + VY_SHIF T; -- oi is t -1 VY_INT := QY_INT - VX_SHIF T; VZ_INT := QZ_INT + ROM_arc tan(N); end if ; else if oi = 1 then VX_INT := QX_INT + VY_SHIF T; VY_INT := QY_INT + VX_SHIF T; VZ_INT := QZ_INT - ROM_art anh(N-1); else VX_INT := QX_INT - VY_SHIF T; VY_INT := QY_INT - VX_SHIF T; VZ_INT := QZ_INT + ROM_art anh(N-1); end if ; end if ; -- Bestimmung von oi case OP_CODE is when X"01" => if VZ_INT(31) = '0' then oi <= 1; -- z >= 0 => oi ist positiv else oi <= -1; -- z < 0 => oi ist negativ end if ; when X"03" => if VZ_INT(31) = '0' then oi <= 1; -- z >= 0 => oi ist positiv else oi <= -1; -- z < 0 => oi ist negativ end if ; when others => if VY_INT(31) = '0' then oi <= -1; -- y >= 0 => oi ist negativ else oi <= 1; -- y < 0 => oi ist positiv end if ; end case ;

51

QX_INT <= VX_INT(31 downto 0); QY_INT <= VY_INT(31 downto 0); QZ_INT <= VZ_INT; when Z4 => case OP_CODE is when X"01" => -- Umwandlung des Zweierkomplements in die Vorzeichen-Betrag-Darstellung if QX_INT(31) = '1' then VX_INT(31 downto 0) := QX_INT(31) & ( not (QX_INT(30 downto 0)) + 1); else VX_INT(31 downto 0) := QX_INT; end if ; if QY_INT(31) = '1' then VY_INT(31 downto 0) := QY_INT(31) & ( not (QY_INT(30 downto 0)) + 1); else VY_INT(31 downto 0) := QY_INT; end if ; QX_INT <= VX_I NT(31 downto 0); -- unskalierte Vorzeichen-Betrag-Darstellung QY_INT <= VY_I NT(31 downto 0); -- unskalierte Vorzeichen-Betrag-Darstellung --Beginn der s eriellen Konstantenmultiplikation V_ROM_X1_OUT : = K1_ROM_X1(Bit6Int(VX_INT(6 downto 1))); V_ROM_Y1_OUT : = K1_ROM_Y1(Bit6Int(VY_INT(6 downto 1))); QX_MULT <= QX_ INT(31) & (("000000" & V_ROM_X1_OUT(31 downto 6)) + K1_ROM_X2(Bit6Int(VX_INT(12 downto 7)))); QY_MULT <= QY_ INT(31) & (("000000" & V_ROM_Y1_OUT(31 downto 6)) + K1_ROM_Y2(Bit6Int(VY_INT(12 downto 7)))); when X"03" => --Addition von QX_INT (cosh) und QY_INT (sinh) zur Exponentialfunktion VX_INT(31 downto 0) := QX_INT + QY_INT; --Umwandlung d es Zweierkomplements in die Vorzeichen-Betrag-Darst ellung if VX_INT(31) = '1' then VX_INT(31 downto 0) := VX_INT(31) & ( not (VX_INT(30 downto 0)) + 1); -- else -- VX_INT(31 downto 0) := QX_INT; end if ; QX_INT <= VX_I NT(31 downto 0); -- unskalierte Vorzeichen-Betrag-Darstellung QY_INT <= "000 00000000000000000000000000000"; -- löschen --Beginn der s eriellen Konstantenmultiplikation V_ROM_X1_OUT : = K2_ROM_X1(Bit6Int(VX_INT(6 downto 1))); QX_MULT <= QX_ INT(31) & (("000000" & V_ROM_X1_OUT(31 downto 6)) + K2_ROM_X2(Bit6Int(VX_INT(12 downto 7)))); QY_MULT <= "00 0000000000000000000000000000000"; when X"04" => -- in QZ_INT steht der Areatangens hype rbolicus QZ_INT <= QZ_I NT(30 downto 0) & '0'; -- Linksshift => 2*artanh() if QREF_EXP_XY(9) = '1' then -- Umwandlung von QREF_EXP_XY in Vorzeich en-Betrag-Darstellung QREF_EXP_XY <= (QREF_EXP_XY(9) & not (QREF_EXP_XY(8 downto 0))) + 1; end if ; QX_MULT <= "00 0000000000000000000000000000000"; QY_MULT <= "00 0000000000000000000000000000000"; when others => -- in QX_INT steht sqrt(r) (unskaliert) -- Umwandlung des Zweierkomplements in die Vorzeichen-Betrag-Dars tellung

52

if QX_INT(31) = '1' then VX_INT(31 downto 0) := QX_INT(31) & ( not (QX_INT(30 downto 0)) + 1); else VX_INT(31 downto 0) := QX_INT; end if ; --Beginn der s eriellen Konstantenmultiplikation V_ROM_X1_OUT : = K2_ROM_X1(Bit6Int(VX_INT(6 downto 1))); QX_MULT <= QX_ INT(31) & (("000000" & V_ROM_X1_OUT(31 downto 6)) + K2_ROM_X2(Bit6Int(VX_INT(12 downto 7)))); QY_MULT <= "00 0000000000000000000000000000000"; end case ; when Z5 => case OP_CODE is when X"01" => QX_MULT(31 downto 0) <= ("000000" & QX_MULT(31 downto 6)) + K1_ROM_X2(Bit6Int(QX_INT(18 downto 13))); QY_MULT(31 downto 0) <= ("000000" & QY_MULT(31 downto 6)) + K1_ROM_Y2(Bit6Int(QY_INT(18 downto 13))); when X"03" => QX_MULT(31 downto 0) <= ("000000" & QX_MULT(31 downto 6)) + K2_ROM_X2(Bit6Int(QX_INT(18 downto 13))); when X"04" => if QREF_EXP_XY(9) = '1' then QE_LN_2 <= ('0' & not (K3_ROM(Bit6Int(QREF_EXP_XY(5 downto 0))))) + 1; else QE_LN_2 <= '0' & K3_ROM(Bit6Int(QREF_EXP_XY(5 downto 0))); end if ; -- Für die Add ition von QE_LN_2 werden in QZ_INT drei zusätzliche Overflow-Bits -- eingefügt, wodurch die drei niederwertigsten Addition-Bits ver loren gehen; if QZ_INT(31) = '1' then QZ_INT <= Q Z_INT(31) & "111" & QZ_INT(30 downto 3); else QZ_INT <= Q Z_INT(31) & "000" & QZ_INT(30 downto 3); end if ; when others => QX_MULT(31 downto 0) <= ("000000" & QX_MULT(31 downto 6)) + K2_ROM_X2(Bit6Int(QX_INT(18 downto 13))); -- Anpassung d es Referenzexponenten (E mod 2 = 1 => Addition von 1) if QE_MOD_2 = '1' then QREF_EXP_XY <= QREF_EXP_XY + 1; end if ; end case ; when Z6 => case OP_CODE is when X"01" => QX_MULT(31 downto 0) <= ("000000" & QX_MULT(31 downto 6)) + K1_ROM_X2(Bit6Int(QX_INT(24 downto 19))); QY_MULT(31 downto 0) <= ("000000" & QY_MULT(31 downto 6)) + K1_ROM_Y2(Bit6Int(QY_INT(24 downto 19))); when X"03" => QX_MULT(31 downto 0) <= ("000000" & QX_MULT(31 downto 6)) + K2_ROM_X2(Bit6Int(QX_INT(24 downto 19))); when X"04" => -- Addition von QZ_INT und QE_LN_2 QE_LN_2(31 downto 0) <= QZ_INT + QE_LN_2(31 downto 0); when others => QX_MULT(31 downto 0) <= ("000000" & QX_MULT(31 downto 6)) + K2_ROM_X2(Bit6Int(QX_INT(24 downto 19))); -- Anpassung d es Referenzexponenten (Shiftoperation) if QREF_EXP_XY(9) = '1' then -- Exponent ist negativ => von links eine '1' nachziehen QREF_EXP_XY <= QREF_EXP_XY(9) & '1' & QREF_EXP_XY(8 downto 1); else QREF_EXP_XY <= QREF_EXP_XY(9) & '0' & QREF_EXP_XY(8 downto 1); end if ; end case ;

53

when Z7 => case OP_CODE is when X"01" => QX_MULT(31 downto 0) <= ("000000" & QX_MULT(31 downto 6)) + K1_ROM_X2(Bit6Int(QX_INT(30 downto 25))); QY_MULT(31 downto 0) <= ("000000" & QY_MULT(31 downto 6)) + K1_ROM_Y2(Bit6Int(QY_INT(30 downto 25))); when X"03" => QX_MULT(31 downto 0) <= ("000000" & QX_MULT(31 downto 6)) + K2_ROM_X2(Bit6Int(QX_INT(30 downto 25))); -- in QX_MULT steht jetzt das skalierte 33-Bit-Ergebnis der CORDI C-Iteration when X"04" => -- Umwandlung des Zweierkomplements in die Vorzeichen-Betrag-Darstellung if QE_LN_2(31) = '1' then QE_LN_2 <= (QE_LN_2(31) & '0' & ( not (QE_LN_2(30 downto 0)))) + 1; else QE_LN_2 <= QE_LN_2(31) & '0' & QE_LN_2(30 downto 0); end if ; when others => QX_MULT(31 downto 0) <= ("000000" & QX_MULT(31 downto 6)) + K2_ROM_X2(Bit6Int(QX_INT(30 downto 25))); end case ; when Z8 => -- Umwandlung in das externe 32-Bit-Gleitk ommaformat -- Korrektur des Referenzexponen ten VREF_EXP_X := QREF_EXP_XY - Coun tX_mult + QQUOTIENT + 127; VREF_EXP_Y := QREF_EXP_XY - Coun tY_mult + QQUOTIENT + 127; VREF_EXP_Z := 127 - CountZ_add + 3; case OP_CODE is when X"01" => case FALLUNTERSCHEIDUNG is when "00" => QXFP_OUT <= QX_MULT(32) & VREF_EXP_X(7 downto 0) & ManXout(22 downto 0); QYZFP_OUT <= QY_MULT(32) & VREF_EXP_Y(7 downto 0) & ManYout(22 downto 0); when "01" => QXFP_OUT <= not (QY_MULT(32)) & VREF_EXP_Y(7 downto 0) & ManYout(22 downto 0); QYZFP_OUT <= QX_MULT(32) & VREF_EXP_X(7 downto 0) & ManXout(22 downto 0); when "10" => QXFP_OUT <= not (QX_MULT(32)) & VREF_EXP_X(7 downto 0) & ManXout(22 downto 0); QYZFP_OUT <= not (QY_MULT(32)) & VREF_EXP_Y(7 downto 0) & ManYout(22 downto 0); when "11" => QXFP_OUT <= QY_MULT(32) & VREF_EXP_Y(7 downto 0) & ManYout(22 downto 0); QYZFP_OUT <= not (QX_MULT(32)) & VREF_EXP_X(7 downto 0) & ManXout(22 downto 0); when others => end case ; when X"03" => QXFP_OUT <= QX_MULT(32) & VREF_EXP_X(7 downto 0) & ManXout(22 downto 0); QYZFP_OUT <= " 00000000000000000000000000000000"; when X"04" => QXFP_OUT <= "00000000000000000000000000 000000"; QYZFP_OUT <= Q E_LN_2(32) & VREF_EXP_Z(7 downto 0) & ManZout(22 downto 0); when others => QXFP_OUT <= QX_MULT(32) & VREF_EXP_X(7 downto 0) & ManXout(22 downto 0); QYZFP_OUT <= " 00000000000000000000000000000000"; end case ; when others => end case ;

54

end process A_BEST; XFP_OUT <= QXFP_OUT; YZFP_OUT <= QYZFP_OUT; end CORDIC_VERSION_1_A;

-- 24_Bit_Barrelshifter.vhd library IEEE; use IEEE.STD_LOGIC_1164. ALL; use IEEE.STD_LOGIC_ARITH. ALL; use IEEE.STD_LOGIC_UNSIGNED.ALL; entity SHFT24LR is port (E: in std_logic_vector(23 downto 0); COUNT: in std_logic_vector(7 downto 0); SHIFT_DIRECTION: in std_logic; OV_BITS: out std_logic_vector(1 downto 0); A: out std_logic_vector (23 downto 0)); end SHFT24LR; architecture SHFT24LR_A of SHFT24LR is signal SELECTOR: std_logic_vector(8 downto 0); begin SELECTOR <= SHIFT_DIRECTION & COUNT; with SELECTOR select A <= E when "000000000", -- Rechtsshift '0' & E(23 downto 1) when "000000001", "00" & E(23 downto 2) when "000000010", "000" & E(23 downto 3) when "000000011", "0000" & E(23 downto 4) when "000000100", "00000" & E(23 downto 5) when "000000101", "000000" & E(23 downto 6) when "000000110", "0000000" & E(23 downto 7) when "000000111", "00000000" & E(23 downto 8) when "000001000", "000000000" & E(23 downto 9) when "000001001", "0000000000" & E(23 downto 10) when "000001010", "00000000000" & E(23 downto 11) when "000001011", "000000000000" & E(23 downto 12) when "000001100", "0000000000000" & E(23 downto 13) when "000001101", "00000000000000" & E(23 downto 14) when "000001110", "000000000000000" & E(23 downto 15) when "000001111", "0000000000000000" & E(23 downto 16) when "000010000", "00000000000000000" & E(23 downto 17) when "000010001", "000000000000000000" & E(23 downto 18) when "000010010", "0000000000000000000" & E(23 downto 19) when "000010011", "00000000000000000000" & E(23 downto 20) when "000010100", "000000000000000000000" & E(23 downto 21) when "000010101", "0000000000000000000000" & E(23 downto 22) when "000010110", "00000000000000000000000" & E(23) when "000010111",

55

E when "100000000", -- Linksshift E(22 downto 0) & '0' when "100000001", -- um eine Stelle E(21 downto 0) & "00" when "100000010", -- um zwei Stellen "000000000000000000000000" when others ; OV_BITS <= '0' & E(23) when SELECTOR = "100000001" else E(23 downto 22) when SELECTOR = "100000010" else "00"; end SHFT24LR_A;

-- CONV_FP_TC.vhd -- Konvertierung einer Gleitkommazahl in das intern e 32-Bit-Festkommaformat (Zweikomplement-Darstellun g) library IEEE; use IEEE.STD_LOGIC_1164. ALL; use IEEE.STD_LOGIC_ARITH. ALL; use IEEE.STD_LOGIC_UNSIGNED.ALL; entity CONV_FP_TC is port (MANTISSE: in std_logic_vector(23 downto 0); VZ: in std_logic; OV_BITS: in std_logic_vector(1 downto 0); S_INT: out std_logic_vector(31 downto 0)); end CONV_FP_TC; architecture CONV_FP_TC_A of CONV_FP_TC is begin with VZ select S_INT <= VZ & OV_BITS & MANTISSE(23 downto 0) & "00000" when '0', (VZ & not (OV_BITS) & not (MANTISSE(23 downto 0)) & "11111") + 1 when '1', "00000000000000000000000000000000" when others ; end CONV_FP_TC_A;

56

-- EXTENDING_DOMAIN_EXP.vhd -- Erweiterung des Argumente-Bereichs der exp-Funkt ion auf +- 16 library IEEE; use IEEE.STD_LOGIC_1164. ALL; use IEEE.STD_LOGIC_ARITH. ALL; use IEEE.STD_LOGIC_UNSIGNED.ALL; entity EXTENDING_DOMAIN is port (EXPONENT: in std_logic_vector(7 downto 0); MANTISSE: in std_logic_vector(23 downto 0); MANTISSE_SHIFT: in std_logic_vector(23 downto 0); QUOTIENT: out std_logic_vector(4 downto 0); REST: out std_logic_vector (23 downto 0)); end EXTENDING_DOMAIN; architecture EXTENDING_DOMAIN_A of EXTENDING_DOMAIN is signal Man_1,Man_2,Man_3,Man_4,Man_5: std_logic_vector(24 downto 0); signal Diff_1,Diff_2,Diff_3,Diff_4,Diff_5: std_logic_vect or(24 downto 0); signal Rest_1,Rest_2,Rest_3,Rest_4,Rest_5: std_logic_vect or(23 downto 0); constant LN_2: std_logic_vector(23 downto 0) := "010110001011100100001100"; begin Man_1 <= '0' & MANTISSE; Diff_1 <= Man_1 - (LN_2 & '0'); Rest_1 <= Diff_1(24 downto 1) when Diff_1(24) = '0' else Man_1(24 downto 1); Man_2 <= Man_1(23 downto 0) & '0' when Diff_1(24) = '1' else Diff_1(23 downto 0) & '0'; Diff_2 <= Man_2 - (LN_2 & '0'); Rest_2 <= Diff_2(24 downto 1) when Diff_2(24) = '0' else Man_2(24 downto 1); Man_3 <= Man_2(23 downto 0) & '0' when Diff_2(24) = '1' else Diff_2(23 downto 0) & '0'; Diff_3 <= Man_3 - (LN_2 & '0'); Rest_3 <= Diff_3(24 downto 1) when Diff_3(24) = '0' else Man_3(24 downto 1); Man_4 <= Man_3(23 downto 0) & '0' when Diff_3(24) = '1' else Diff_3(23 downto 0) & '0'; Diff_4 <= Man_4 - (LN_2 & '0'); Rest_4 <= Diff_4(24 downto 1) when Diff_4(24) = '0' else Man_4(24 downto 1); Man_5 <= Man_4(23 downto 0) & '0' when Diff_4(24) = '1' else Diff_4(23 downto 0) & '0'; Diff_5 <= Man_5 - (LN_2 & '0'); Rest_5 <= Diff_5(24 downto 1) when Diff_5(24) = '0' else Man_5(24 downto 1); with EXPONENT select REST <= Rest_1 when "01111110", Rest_2 when "01111111", Rest_3 when "10000000", Rest_4 when "10000001", Rest_5 when "10000010", --"000000000000000000000000" when others ; MANTISSE_SHIFT when others ; with EXPONENT select QUOTIENT <= "0000" & ( not (Diff_1(24))) when "01111110", "000" & ( not (Diff_1(24) & Diff_2(24))) when "01111111", "00" & ( not (Diff_1(24) & Diff_2(24) & Diff_3(24))) when "10000000",

57

'0' & ( not (Diff_1(24) & Diff_2(24) & Diff_3(24) & Diff_4(24)) ) when "10000001", not (Diff_1(24) & Diff_2(24) & Diff_3(24) & Diff_4(24) & Diff_5(24)) when "10000010", "00000" when others ; end EXTENDING_DOMAIN_A;

-- COUNT_SH_OP.vhd -- Bestimmung der Anzahl der Schiebeoperationen zur Mantissen-Normalisierung library IEEE; use IEEE.STD_LOGIC_1164. ALL; use IEEE.STD_LOGIC_ARITH. ALL; use IEEE.STD_LOGIC_UNSIGNED.ALL; entity COUNT_SH_OP is port (E: in std_logic_vector(32 downto 0); COUNT: out std_logic_vector(9 downto 0)); end COUNT_SH_OP; architecture COUNT_SH_OP_A of COUNT_SH_OP is begin COUNT <= "1111111101" when E(31)='1' else "1111111110" when E(31)='0' and E(30)='1' else "1111111111" when E(31 downto 30)="00" and E(29)='1' else "0000000000" when E(31 downto 29)="000" and E(28)='1' else "0000000001" when E(31 downto 28)="0000" and E(27)='1' else "0000000010" when E(31 downto 27)="00000" and E(26)='1' else "0000000011" when E(31 downto 26)="000000" and E(25)='1' else "0000000100" when E(31 downto 25)="0000000" and E(24)='1' else "0000000101" when E(31 downto 24)="00000000" and E(23)='1' else "0000000110" when E(31 downto 23)="000000000" and E(22)='1' else "0000000111" when E(31 downto 22)="0000000000" and E(21)='1' else "0000001000" when E(31 downto 21)="00000000000" and E(20)='1' else "0000001001" when E(31 downto 20)="000000000000" and E(19)='1' else "0000001010" when E(31 downto 19)="0000000000000" and E(18)='1' else "0000001011" when E(31 downto 18)="00000000000000" and E(17)='1' else "0000001100" when E(31 downto 17)="000000000000000" and E(16)='1' else "0000001101" when E(31 downto 16)="0000000000000000" and E(15)='1' else "0000001110" when E(31 downto 15)="00000000000000000" and E(14)='1' else "0000001111" when E(31 downto 14)="000000000000000000" and E(13)='1' else "0000010000" when E(31 downto 13)="0000000000000000000" and E(12)='1' else "0000010001" when E(31 downto 12)="00000000000000000000" and E(11)='1' else "0000010010" when E(31 downto 11)="000000000000000000000" and E(10)='1' else "0000010011" when E(31 downto 10)="0000000000000000000000" and E(9)='1' else "0000010100" when E(31 downto 9) ="00000000000000000000000" and E(8)='1' else "0000010101" when E(31 downto 8) ="000000000000000000000000" and E(7)='1' else "0000010110" when E(31 downto 7) ="0000000000000000000000000" and E(6)='1' else "0000010111" when E(31 downto 6) ="00000000000000000000000000" and E(5)='1' else "0000011000"; end COUNT_SH_OP_A;

58

-- 32_Bit_Barrelshifter.vhd library IEEE; use IEEE.STD_LOGIC_1164. ALL; use IEEE.STD_LOGIC_ARITH. ALL; use IEEE.STD_LOGIC_UNSIGNED.ALL; entity SHFT32LR is port (E: in std_logic_vector(32 downto 0); COUNT: in std_logic_vector(9 downto 0); A: out std_logic_vector (23 downto 0)); end SHFT32LR; architecture SHFT32LR_A of SHFT32LR is begin with COUNT select A <= E(31 downto 8) when "1111111101", E(30 downto 7) when "1111111110", E(29 downto 6) when "1111111111", E(28 downto 5) when "0000000000", E(27 downto 4) when "0000000001", E(26 downto 3) when "0000000010", E(25 downto 2) when "0000000011", E(24 downto 1) when "0000000100", E(23 downto 0) when "0000000101", E(22 downto 0) & '0' when "0000000110", E(21 downto 0) & "00" when "0000000111", E(20 downto 0) & "000" when "0000001000", E(19 downto 0) & "0000" when "0000001001", E(18 downto 0) & "00000" when "0000001010", E(17 downto 0) & "000000" when "0000001011", E(16 downto 0) & "0000000" when "0000001100", E(15 downto 0) & "00000000" when "0000001101", E(14 downto 0) & "000000000" when "0000001110", E(13 downto 0) & "0000000000" when "0000001111", E(12 downto 0) & "00000000000" when "0000010000", E(11 downto 0) & "000000000000" when "0000010001", E(10 downto 0) & "0000000000000" when "0000010010", E(9 downto 0) & "00000000000000" when "0000010011", E(8 downto 0) & "000000000000000" when "0000010100", E(7 downto 0) & "0000000000000000" when "0000010101", E(6 downto 0) & "00000000000000000" when "0000010110", E(5 downto 0) & "000000000000000000" when "0000010111", "000000000000000000000000" when others ; end SHFT32LR_A;

59

-- PACK_CORDIC.vhd library IEEE; use IEEE.STD_LOGIC_1164. ALL; use IEEE.STD_LOGIC_ARITH. ALL; use IEEE.STD_LOGIC_UNSIGNED.ALL; -- Package Deklaration package PACK_CORDIC is type MMX24 is array (23 downto 0) of std_logic_vector(31 downto 0); type MMX64_32 is array (63 downto 0) of std_logic_vector(31 downto 0); type MMX41_1 is array (147 downto 107) of std_logic; constant ROM_arctan: MMX24:= ("00000000000000000000000000011111", "00000000000000000000000000111111", "00000000000000000000000001111111", "00000000000000000000000100000000", "00000000000000000000000111111111", "00000000000000000000010000000000", "00000000000000000000100000000000", "00000000000000000001000000000000", "00000000000000000010000000000000", "00000000000000000100000000000000", "00000000000000000111111111111111", "00000000000000001111111111111111", "00000000000000011111111111111111", "00000000000000111111111111111111", "00000000000001111111111111111111", "00000000000011111111111111111010", "00000000000111111111111111010101", "00000000001111111111111010101010", "00000000011111111111010101010110", "00000000111111111010101011011101", "00000001111111010101101110101001", "00000011111010110110111010111111", "00000111011010110001100111000001", "00001100100100001111110110101010"); constant ROM_artanh: MMX24:= ("00000000000000000000000000000000", "00000000000000000000000000000000", "00000000000000000000000000111111", "00000000000000000000000001111111", "00000000000000000000000100000000", "00000000000000000000000111111111", "00000000000000000000010000000000", "00000000000000000000100000000000", "00000000000000000001000000000000", "00000000000000000010000000000000", "00000000000000000100000000000000", "00000000000000000111111111111111", "00000000000000001111111111111111", "00000000000000100000000000000000", "00000000000001000000000000000000", "00000000000010000000000000000000",

60

"00000000000100000000000000000101", "00000000001000000000000000101010", "00000000010000000000000101010101", "00000000100000000000101010101100", "00000001000000000101010110001000", "00000010000000101011000100100011", "00000100000101100010101110111110", "00001000110010011111010100111101"); constant K1_ROM_X1,K1_ROM_X2,K1_ROM_Y1,K1_ROM_Y2: MMX64_32:= -- Skalierung, zirkularer Modus ("01001100100000111000110011100101", -- 111111 "01001011010011001010001100001010", -- 111110 "01001010000101011011100100101111", -- 111101 "01001000110111101100111101010100", -- 111100 "01000111101001111110010101111001", -- 111011 "01000110011100001111101110011110", -- 111010 "01000101001110100001000111000011", -- 111001 "01000100000000110010011111101000", -- 111000 "01000010110011000011111000001101", -- 110111 "01000001100101010101010000110010", -- 110110 "01000000010111100110101001010111", -- 110101 "00111111001001111000000001111100", -- 110100 "00111101111100001001011010100001", -- 110011 "00111100101110011010110011000110", -- 110010 "00111011100000101100001011101011", -- 110001 "00111010010010111101100100010000", -- 110000 "00111001000101001110111100110101", -- 101111 "00110111110111100000010101011010", -- 101110 "00110110101001110001101101111111", -- 101101 "00110101011100000011000110100100", -- 101100 "00110100001110010100011111001001", -- 101011 "00110011000000100101110111101110", -- 101010 "00110001110010110111010000010011", -- 101001 "00110000100101001000101000111000", -- 101000 "00101111010111011010000001011101", -- 100111 "00101110001001101011011010000010", -- 100110 "00101100111011111100110010100111", -- 100101 "00101011101110001110001011001100", -- 100100 "00101010100000011111100011110001", -- 100011 "00101001010010110000111100010110", -- 100010 "00101000000101000010010100111011", -- 100001 "00100110110111010011101101100000", -- 100000 "00100101101001100101000110000101", -- 011111 "00100100011011110110011110101010", -- 011110 "00100011001110000111110111001111", -- 011101 "00100010000000011001001111110100", -- 011100 "00100000110010101010101000011001", -- 011011 "00011111100100111100000000111110", -- 011010 "00011110010111001101011001100011", -- 011001 "00011101001001011110110010001000", -- 011000 "00011011111011110000001010101101", -- 010111 "00011010101110000001100011010010", -- 010110 "00011001100000010010111011110111", -- 010101 "00011000010010100100010100011100", -- 010100 "00010111000100110101101101000001", -- 010011 "00010101110111000111000101100110", -- 010010 "00010100101001011000011110001011", -- 010001

61

"00010011011011101001110110110000", -- 010000 "00010010001101111011001111010101", -- 001111 "00010001000000001100100111111010", -- 001110 "00001111110010011110000000011111", -- 001101 "00001110100100101111011001000100", -- 001100 "00001101010111000000110001101001", -- 001011 "00001100001001010010001010001110", -- 001010 "00001010111011100011100010110011", -- 001001 "00001001101101110100111011011000", -- 001000 "00001000100000000110010011111101", -- 000111 "00000111010010010111101100100010", -- 000110 "00000110000100101001000101000111", -- 000101 "00000100110110111010011101101100", -- 000100 "00000011101001001011110110010001", -- 000011 "00000010011011011101001110110110", -- 000010 "00000001001101101110100111011011", -- 000001 "00000000000000000000000000000000"); -- 000000 constant K2_ROM_X1,K2_ROM_X2: MMX64_32:= -- Skalierung, hyperbolischer Modus ("10011000001001010000011001110010", -- 111111 "10010101101110101100100101100100", -- 111110 "10010011010100001000110001010110", -- 111101 "10010000111001100100111101001000", -- 111100 "10001110011111000001001000111010", -- 111011 "10001100000100011101010100101100", -- 111010 "10001001101001111001100000011110", -- 111001 "10000111001111010101101100010000", -- 111000 "10000100110100110001111000000010", -- 110111 "10000010011010001110000011110100", -- 110110 "01111111111111101010001111100110", -- 110101 "01111101100101000110011011011000", -- 110100 "01111011001010100010100111001010", -- 110011 "01111000101111111110110010111100", -- 110010 "01110110010101011010111110101110", -- 110001 "01110011111010110111001010100000", -- 110000 "01110001100000010011010110010010", -- 101111 "01101111000101101111100010000100", -- 101110 "01101100101011001011101101110110", -- 101101 "01101010010000100111111001101000", -- 101100 "01100111110110000100000101011010", -- 101011 "01100101011011100000010001001100", -- 101010 "01100011000000111100011100111110", -- 101001 "01100000100110011000101000110000", -- 101000 "01011110001011110100110100100010", -- 100111 "01011011110001010001000000010100", -- 100110 "01011001010110101101001100000110", -- 100101 "01010110111100001001010111111000", -- 100100 "01010100100001100101100011101010", -- 100011 "01010010000111000001101111011100", -- 100010 "01001111101100011101111011001110", -- 100001 "01001101010001111010000111000000", -- 100000 "01001010110111010110010010110010", -- 011111 "01001000011100110010011110100100", -- 011110 "01000110000010001110101010010110", -- 011101 "01000011100111101010110110001000", -- 011100 "01000001001101000111000001111010", -- 011011 "00111110110010100011001101101100", -- 011010

62

"00111100010111111111011001011110", -- 011001 "00111001111101011011100101010000", -- 011000 "00110111100010110111110001000010", -- 010111 "00110101001000010011111100110100", -- 010110 "00110010101101110000001000100110", -- 010101 "00110000010011001100010100011000", -- 010100 "00101101111000101000100000001010", -- 010011 "00101011011110000100101011111100", -- 010010 "00101001000011100000110111101110", -- 010001 "00100110101000111101000011100000", -- 010000 "00100100001110011001001111010010", -- 001111 "00100001110011110101011011000100", -- 001110 "00011111011001010001100110110110", -- 001101 "00011100111110101101110010101000", -- 001100 "00011010100100001001111110011010", -- 001011 "00011000001001100110001010001100", -- 001010 "00010101101111000010010101111110", -- 001001 "00010011010100011110100001110000", -- 001000 "00010000111001111010101101100010", -- 000111 "00001110011111010110111001010100", -- 000110 "00001100000100110011000101000110", -- 000101 "00001001101010001111010000111000", -- 000100 "00000111001111101011011100101010", -- 000011 "00000100110101000111101000011100", -- 000010 "00000010011010100011110100001110", -- 000001 "00000000000000000000000000000000"); -- 000000 constant K3_ROM: MMX64_32:= -- Produkt E*ln2 für die Berechnung des na türlichen Logarithmus ("01010111010101100010011101010010", -- 111111 "01010101111100110100001100100100", -- 111110 "01010100100100000101111011110110", -- 111101 "01010011001011010111101011001000", -- 111100 "01010001110010101001011010011010", -- 111011 "01010000011001111011001001101100", -- 111010 "01001111000001001100111000111110", -- 111001 "01001101101000011110101000010000", -- 111000 "01001100001111110000010111100010", -- 110111 "01001010110111000010000110110100", -- 110110 "01001001011110010011110110000110", -- 110101 "01001000000101100101100101011000", -- 110100 "01000110101100110111010100101010", -- 110011 "01000101010100001001000011111100", -- 110010 "01000011111011011010110011001110", -- 110001 "01000010100010101100100010100000", -- 110000 "01000001001001111110010001110010", -- 101111 "00111111110001010000000001000100", -- 101110 "00111110011000100001110000010110", -- 101101 "00111100111111110011011111101000", -- 101100 "00111011100111000101001110111010", -- 101011 "00111010001110010110111110001100", -- 101010 "00111000110101101000101101011110", -- 101001 "00110111011100111010011100110000", -- 101000 "00110110000100001100001100000010", -- 100111 "00110100101011011101111011010100", -- 100110 "00110011010010101111101010100110", -- 100101 "00110001111010000001011001111000", -- 100100 "00110000100001010011001001001010", -- 100011

63

"00101111001000100100111000011100", -- 100010 "00101101101111110110100111101110", -- 100001 "00101100010111001000010111000000", -- 100000 "00101010111110011010000110010010", -- 011111 "00101001100101101011110101100100", -- 011110 "00101000001100111101100100110110", -- 011101 "00100110110100001111010100001000", -- 011100 "00100101011011100001000011011010", -- 011011 "00100100000010110010110010101100", -- 011010 "00100010101010000100100001111110", -- 011001 "00100001010001010110010001010000", -- 011000 "00011111111000101000000000100010", -- 010111 "00011110011111111001101111110100", -- 010110 "00011101000111001011011111000110", -- 010101 "00011011101110011101001110011000", -- 010100 "00011010010101101110111101101010", -- 010011 "00011000111101000000101100111100", -- 010010 "00010111100100010010011100001110", -- 010001 "00010110001011100100001011100000", -- 010000 "00010100110010110101111010110010", -- 001111 "00010011011010000111101010000100", -- 001110 "00010010000001011001011001010110", -- 001101 "00010000101000101011001000101000", -- 001100 "00001111001111111100110111111010", -- 001011 "00001101110111001110100111001100", -- 001010 "00001100011110100000010110011110", -- 001001 "00001011000101110010000101110000", -- 001000 "00001001101101000011110101000010", -- 000111 "00001000010100010101100100010100", -- 000110 "00000110111011100111010011100110", -- 000101 "00000101100010111001000010111000", -- 000100 "00000100001010001010110010001010", -- 000011 "00000010110001011100100001011100", -- 000010 "00000001011000101110010000101110", -- 000001 "00000000000000000000000000000000"); -- 000000 constant E_MOD_2: MMX41_1:= -- Operation E mod 2 ('0', -- 10010011 ( 20) (entspricht 147d) '1', -- 10010010 ( 19) '0', -- 10010001 ( 18) '1', -- 10010000 ( 17) '0', -- 10001111 ( 16) '1', -- 10001110 ( 15) '0', -- 10001101 ( 14) '1', -- 10001100 ( 13) '0', -- 10001011 ( 12) '1', -- 10001010 ( 11) '0', -- 10001001 ( 10) '1', -- 10001000 ( 9) '0', -- 10000111 ( 8) '1', -- 10000110 ( 7) '0', -- 10000101 ( 6) '1', -- 10000100 ( 5) '0', -- 10000011 ( 4) '1', -- 10000010 ( 3) '0', -- 10000001 ( 2) '1', -- 10000000 ( 1)

64

'0', -- 01111111 ( 0) '1', -- 01111110 ( -1) '0', -- 01111101 ( -2) '1', -- 01111100 ( -3) '0', -- 01111011 ( -4) '1', -- 01111010 ( -5) '0', -- 01111001 ( -6) '1', -- 01111000 ( -7) '0', -- 01110111 ( -8) '1', -- 01110110 ( -9) '0', -- 01110101 (-10) '1', -- 01110100 (-11) '0', -- 01110011 (-12) '1', -- 01110010 (-13) '0', -- 01110001 (-14) '1', -- 01110000 (-15) '0', -- 01101111 (-16) '1', -- 01101110 (-17) '0', -- 01101101 (-18) '1', -- 01101100 (-19) '0'); -- 01101011 (-20) (entspricht 107d) function Bit5Int (b: std_logic_vector(4 downto 0)) return INTEGER; function Bit6Int (b: std_logic_vector(5 downto 0)) return INTEGER; function Bit8Int (b: std_logic_vector(7 downto 0)) return INTEGER; end PACK_CORDIC; package body PACK_CORDIC is function Bit5Int (b: std_logic_vector(4 downto 0)) return INTEGER is variable value: INTEGER; begin value:=0; for I in 0 to 4 loop if (b(I)='1') then value:=value + (2**I); end if ; end loop ; return value; end Bit5Int; function Bit6Int (b: std_logic_vector(5 downto 0)) return INTEGER is variable value: INTEGER; begin value:=0; for I in 0 to 5 loop if (b(I)='1') then value:=value + (2**I); end if ; end loop ; return value; end Bit6Int; function Bit8Int (b: std_logic_vector(7 downto 0)) return INTEGER is

65

variable value: INTEGER; begin value:=0; for I in 0 to 7 loop if (b(I)='1') then value:=value + (2**I); end if ; end loop ; return value; end Bit8Int; end PACK_CORDIC;

66

Literaturverzeichnis [Andr98] Andraka, R.: A survey of CORDIC algorithms for FPGA based comput ers ; http://www.andraka.com/files/crdcsrvy.pdf, 1998. [Beul08] Beuler, M.: Realisierung von Arithmetik-Baugruppen für das 32-B it-Gleit-

kommaformat der Norm ANSI/IEEE 754 mittels VHDL ; Report, Fachhochschule Gießen-Friedberg, http://dok.bib.fh-giessen.de/opus/volltexte/2009/4146/, April 2008.

[Hoff93] Hoffmann, R.: Rechnerentwurf: Rechenwerke, Mikroprogrammierung, R ISC; 3. Auflage, Oldenbourg Verlag, 1993. [Jork04] Jorke, G.: Rechnergestützter Entwurf digitaler Schaltungen – S chaltungs-

synthese mit VHDL ; 1. Auflage, Carl Hanser Verlag, 2004. [Papu01] Papula, L.: Mathematische Formelsammlung ; 7. Auflage, Vieweg Verlag, 2001. [Pirs96] Pirsch, P.: Architekturen der digitalen Signalverarbeitung ; 1. Auflage, Teubner Verlag, 1996. [Riss04] Risse, T.: CORDIC-Algorithmen Verbinden Mathematik, Computer-A rchitektur

und Anwendungen ; Global J. of Engng. Educ., Vol.8, No.3; 2004. [Timm93] Timmermann, D.; Rix, B.; Hahn, H.; Hosticka, B.J.: 24 Bit CMOS Gleitkomma

Vektorarithmetik Chip ; Fraunhofer-Institut für Mikroelektronische Schaltungen und Systeme, 1993.

[Walt71] Walther, J. S.: A unified algorithm for elementary functions ; Proc. Spring Joint

Computer Conf., pp 379-385, 1971. [WEB01] http://de.wikipedia.org/wiki/CORDIC [WEB02] Risse, T.: http://www.weblearn.hs-bremen.de/risse/MAI/docs/numerik.pdf Numerik – Interaktiv; 2004.