Dipl.-lng. Anton Obermaier

Kleine Standalone-Systeme (wie zum Beispiel der Z80-EMUF) werden normalerweise nicht auf Fließkommarechnungen ausgelegt. Es gibt aber Methoden, nur mit Integerarithmetik Funktionen wie sin und cos im Rahmen solcher Minimalkonzepte mit hinreichender Genauigkeit zu implementieren.

Wenn man zum Beispiel den Z80 in einem kleinen Einplatinensystem einsetzt, dann muß man, will man nicht ein viel Speicherplatz fressendes Arithmetikpaket implementieren, schon für die gewöhnliche Multiplikation und Division von ganzen Zahlen eigene Routinen entwerfen. In der Praxis sollen aber oft auch höhere Funktionen von solchen Systemen berechnet werden. Bei geschickter Programmierung gelingt das auch - und solche Spezial-Routinen liefern ihr Ergebnis auch schneller ab, als allgemeine Gleitkommapakete das können. Die Schnelligkeit der Routinen und ihr geringer Speicherbedarf werden allerdings mit begrenzter Genauigkeit (hier vier Stellen) erkauft. In der Praxis ist aber oft sehr viel weniger Genauigkeit gefragt, als Gleitkommapakete sie liefern.
Als Beispiel soll angenommen werden, daß eine Positionierungsaufgabe im zweidimensionalen Raum gegeben ist, wobei Polarkoordinaten verwendet werden sollen. Hier wird die Sinus-Funktion mit Ausgabewerten zwischen -1 und +1 gebraucht, die zunächst nicht mit „Integers" berechenbar ist.

Trigonometrische Funktionen lassen sich mit Integer-Operationen programmieren, wenn man sich geeignete Datenstrukturen für die Argumente und die Zahlen aus dem Wertebereich der Funktion zurechtlegt. Jede der vorkommenden Größen ist in der Praxis beschränkt, es gibt einen Wert mit maximalem Betrag (wie dieser Wert für eine Variable aussieht, das ist je nach Einsatzfall verschieden, aber immer für eine solche Programmieraufgabe ermittelbar). Eine Lösung für die Darstellung solcher Größen verläuft so: Man nimmt die Integerzahl mit maximalem Betrag (sie sei mit Im bezeichnet), die man im gewählten Computersystem darstellen kann, ordnet ihr das Maximum m aus dem Wertebereich der Variablen zu und stellt alle kleineren Zahlen x aus dem Wertebereich im Computer dar, indem man diejenige Integerzahl Ix im Computer aufsucht, die sich aus der Formel

^Ix/_Im = ^x/_m ergibt.

Also: Ix = ^xIm/_m,

wobei hier einfach später die Nachkommastellen weggelassen werden. Das ist so etwas wie „Normalisierung". In Bezug auf Trigonometrie ergibt sich, wenn man mit 16-Bit-Worten umgeht (was beim Z80 naheliegend ist), nach diesem Schema eine natürliche Datenstruktur sowohl für Argumentwerte als auch für die Werte der trigonometrischen Funktionen:

1. Die Darstellung der Eingangsgröße (Argument = Winkel) als Wort liefert, wenn man die zugelassenen Werte auf 0...360 Grad beschränkt, eine Genauigkeit von 360/65536 Grad, was 19,8 Sekunden entspricht. Hier ist die Periodizität der trigonometrischen Funktionen sogar noch bei der Addition und Subtraktion mitberücksichtigt. Beispiel:
   10 Grad wird dargestellt als die „Maschinenzahl", die dem Wert
   (10/360)*65536 = 1820,4393
   entspricht, also als 1820 dez = 071CH im Computer.

2. Die Werte von sin(x) und cos(x) müssen ebenfalls normalisiert dargestellt werden. Wenn hier 1 das Maximum m ist und zum Beispiel x = sin(10 Grad), dann ist
   (sin(10 Grad)/1) = 0,001011000111010000.... bin,
   was als
   0001011000111010 bin = 163AH
   gespeichert wird. m = 2^15, weil cos(0°) = 1 dargestellt werden muß. Der relative Fehler beträgt dann -1,8 E-5.

Zur praktischen Berechnung genügt es, wenn jeder vorgegebene Winkel auf den ersten Quadranten reduziert wird. Bild 1 zeigt, daß dort alle vorkommenden Werte positiv sind (0 =< x < PI/2 und 0 =< sin(x), cos(x) =< 1).

Bild 1. Alle Winkel lassen sich auf einen positiven Hauptbereich von 0 bis PI/2 reduzieren

Die Algorithmen zur Berechnung des Betrages der beiden Funktionen gestalten sich hier besonders einfach. Die notwendigen Vorzeichen werden gesondert hinzugefügt. Im vorliegenden Konzept werden nun sin und cos simultan berechnet. Das ist sinnvoll, weil in den meisten Fällen ohnehin beide Werte berechnet werden müssen und sowohl Code-Länge als auch Zeitverhalten in diesem Fall günstig sind.
Zum Beispiel dient die Formel

sin(x) = - cos(x-270Grad)

als Ausgangspunkt für Berechnungen von sin, wenn x im 4. Quadranten liegt (anhand einer Tabelle werden alle vorkommenden Fälle später bearbeitet), woraus sich Reduktion und anschließend Vorzeichen und Vertauschung der Ergebnisse ableiten lassen. Tan(x) kann von hier aus mit einer einzigen Division errechnet werden. Die zur Berechnung noch benötigten weiteren Variablen werden in ähnlicher Weise normiert, wie schon geschehen.

Sin und cos kann man nach Lehrbuch in eine Reihe entwickeln:

sin(x) = x - (x^3)/3! + (x^5)/5! - (x^7)/7! + ...

Das kann man umformen in:

sin(x) = x*(1 - (x^2)/3! + (x^4)/5! - (x^6)/7! + ... )

während für cos gilt:

cos(x) = 1 - (x^2)/2! + (x^4)/4! - (x^6)/6! + ...

wobei x immer im Bogenmaß anzugeben ist. In den letzten beiden Formeln taucht x immer mit geraden Exponenten auf, was die gemeinsame Berechnung von sin und cos sehr vereinfacht.
Die Genauigkeit bei den einzelnen Gliedern hängt in erster Linie von der verwendeten Stellenzahl bei der Quadrierung ab. Normalisiert man hier wie vorhin geschildert, ergeben sich die in Bild 2 gezeigten Faktoren.

Wertebereich Maximum binär Exponent zur Normalisierung 0 <= x <= 1,57... x^2 <= 2,46... x^4 <= 6,08... x^6 <= 15,0... x^8 <= 37,0... 1, ... 10, ... 110, ... 1111, ... 100101, ... 2^15 2^14 2^13 2^12 2^10 Bild 2. Die Faktoren der Normalisierung ergeben sich aus der Wertebeschränkung des Winkels x

Diese Variablen seien im folgenden mit xibog bezeichnet, wobei i die Potenz bei x angebe. Die Fehler bei den xibog sind dann recht gering, wenn sie nach Bild 2 dargestellt werden. Der Gesamtfehler einer Berechnung von sin und cos mit Reihenentwicklung hängt dann im wesentlichen vom Abbruchkriterium während der iterativen Berechnung der Reihe ab.
Da das Glied mit x^10 in beiden Reihen oben wegen des starken Anwachsens der Fakultät bei kleinen x nur noch Beiträge von der Größenordnung 2,5 E-5 (beim sin sogar 2,3 E-6) liefert, werden alle Summanden ab einschließlich x10bog nicht mehr berücksichtigt.

Mit den Variablen

x1bog = x * 2 ^ 15
x2bog = x ^ 2 * 2 ^ 14
x4bog = x ^ 4 * 2 ^ 13
x6bog = x ^ 6 * 2 ^ 12
x8bog = x ^ 8 * 2 ^ 10

ist

sx:=(1/(2^15))*((x1bog)*((2^15)-(x2bog)/3+(x4bog)/(2*3*5)-(x6bog)/(3*5*6*7)+(x8bog)/(2*3*5*6*7*9))) ~ sin(x)*2^15

und

cx:=(2^15)-(x2bog)+(x4bog)/(2*3)-(x6bog)/(3*5*6)+(x8bog)/(2*3*5*6*7) ~ cos(x)*2^15

Diese Formeln sind in den Listings realisiert. Dabei sind die Klammern so gesetzt, daß bei der Umsetzung in das Programm die Genauigkeit erhalten bleibt.

Entsprechend den drei Teilproblemen gibt es drei Module, XICALC, SINUS und COSINUS, die linear aufgebaut sind. Dazu kommen noch einige Hilfs-Routinen, die etwas Rechenarbeit leisten. Bild 3 zeigt einen Ablaufplan.

| (a) Winkel in de bereitstellen | (b) 'XICALC' : xibog := fi(<de>) | (c) 'SINUS' : de := sx(xibog) | (c) push de | (d) 'COSINUS' : hl := cx(xibog) | (d) pop de      {'sx'} | (e) Auswertung : sx / cx -> sin(x) / cos(x) | Bild 3. Der Ablaufplan

Im Z80-Registerpaar DE wird der Winkel erwartet, wobei die Normalisierung wie besprochen durchgeführt sein soll. XICALC berechnet daraus die fünf Variablen xibog. Die Reduzierung der Argumente auf den ersten Quadranten besteht hier aufgrund der Darstellung der Argumentwerte in der Ausmaskierung der beiden höchstwertigen Bits des Winkels, also der Bits 7 und 6 in Register D des Z80. Der Winkel x wird mit

x [Grad] = <DE> * 360 / 2^16
in
x1bog = x * 2^15 = <DE> * PI = (<DE> * (PI * 2^14))/2^14

umgewandelt, wobei auch PI normalisiert angegeben sein muß.

Die Divisionen im Programm sind soweit als möglich auf Rotationen der Register zurückgeführt. Die Unterprogramme MULTI und DIVI halten Parameterübergabebedingungen ein, die alles sehr schnell machen. Das Ergebnis der Multiplikation beim Umrechnen nach den Formeln oben steht in DE-HL (32 Bit). Division durch 2^14 bedeutet Stellenverschiebung um 14 Positionen nach rechts, was auch so aufgefaßt werden kann, daß die Bits 14 bis 29 aus DE-HL in ein 16 Bit breites Register geladen werden (die höchstwertigen beiden Bits werden wegmaskiert). Mit Linksverschiebung um zwei Stellen steht das Ergebnis in DE. Bild 4 zeigt das Listing. Bild 5 zeigt die Routine SINUS, die die Werte xibog aus dem RAM nimmt und dann sx in DE ausgibt. Das Ergebnis sx muß vor dem Aufruf von COSINUS zwischengespeichert werden, was einfach wegen der einheitlichen Struktur der Module nicht automatisch von SINUS gemacht wird.
Bild 6 zeigt die Routine COSINUS, die ebenfalls auf die xibog zurückgreift. Das Ergebnis cx wird in HL geliefert. Nachdem sx wieder in DE installiert ist, stehen beide Werte in den Z80-Registern zur Verfügung.

Bei der Reduktion auf den ersten Quadranten müssen vier Fälle beachtet werden, die den vier Quadranten entsprechen und die aus den bei x wegmaskierten beiden höchstwertigen Bits abgeleitet werden können. Bild 7 zeigt tabellarisch, wie sx und cx bezüglich Vorzeichen und Zuordnung zu sin und cos interpretiert werden müssen.

Quadrant Bit 7 / 6 2^15*sin(x) 2^15*cos(x) 1 2 3 4 0 0 0 1 1 0 1 1 <de> + <hl> - <de> - <hl> <hl> - <de> - <hl> + <de> Bild 7. Die Zuordnung der Ergebniswerte zum Argumentwert

Das Beispiel am Anfang mit sin(x) = - cos (x - 270 Grad) könnte verbal so interpretiert werden: Der Sinus eines Winkels im 4. Quadranten ist mit dem negativen Cosinus des auf den 1. Quadranten reduzierten Winkels identisch. Sin(x) ist also in solch einem Fall der Wert cx in HL, der noch dazu mit dem negativen Vorzeichen versehen werden muß.

Tests mit den hier abgebildeten Programmen haben gezeigt, daß Sinus und Cosinus bei 4 MHz Taktfrequenz in 4 ms berechnet werden können. Eine Gleitkomma-Rechnung auf dem 68000 nach dem Gleitkommapaket aus mc, Ausgabe 2/1985, dauerte 19 mal so lang (allerdings bei doppelter Zahl gültiger Stellen). Das Beispiel mit sin(10 Grad) liefert als Ergebnis

<DE> = 1638H = 5688 dez.

Umgerechnet ist das 0,17358398 ..., was einen relativen Fehler von -3.7 E-4 bedeutet. Der mitberechnete Cosinus lag mit -2,5 E-5 Abweichung noch genauer. Die Länge der Quelle ist mit 200 Byte noch recht handlich. Die verwendeten externen Module zur Multiplikation und Division von 16-Bit-Zahlen sind ohnehin notwendig, müssen aber bestimmten Ein- und Ausgabebedingungen genügen. Bild 8 zeigt mögliche Varianten mit den vorgeschriebenen Merkmalen.