In der „PC AMSTRAD INTERNATIONAL 3/90" wurde der folgende Artikel abgedruckt. (Der Artikel ist ursprünglich für den CPC geschrieben)
Aufgaben mit reellen Zahlen ausgeführt mit ganzen Zahlen.
(Als Zusatz ein Ausschnitt aus einem Artikel der CPC 7/87, der die Berechnung von Quadratwurzeln vorstellt. Auch im Internet findet sich eine Lösung.)
Grundlegende Artikel zur Berechnung von Sinus- und Cosinus-Funktionen mittels Integerarithmetik finden sich hier und hier. Ein anderer Artikel berechnet die Wurzel-Funktion.

Die Assemblerecke

Schnelle Tabellen

Der Z80-Prozessor ist leider kein Rechengenie. Eine Addition oder Subtraktion wird zwar mit einem einzigen Maschinenbefehl erledigt. Bereits eine schlichte Ganzzahl-Multiplikation erfordert jedoch eine komplizierte Routine, die bis zum Endergebnis mehrere Schleifendurchläufe braucht und den Betrieb dementsprechend lange aufhält. Noch schlimmer sieht es bei höheren mathematischen Funktionen wie Sinus und Cosinus aus, die nicht ohne zeitaufwendige Fließkomma-Arithmetik auskommen. Es gibt jedoch ein Verfahren, mit dem sich in vielen Fällen eine dramatische Temposteigerung erreichen läßt: tabellenorientierte Arithmetik.

Die Grundidee ist sehr einfach. Dem Programm werden die Rechenergebnisse von vornherein in Tabellenform zur Verfügung gestellt. Es braucht dann nur noch nachzuschlagen, was sich problemlos mit wenigen Maschinenbefehlen erledigen läßt. Die Tabelle belegt natürlich einigen Speicherplatz, so daß hier letztendlich ein Tauschhandel 'Mikrosekunden gegen Kilobytes' stattfindet. Angenommen, in einer zeitkrischen Schleife muß wiederholt das Quadrat einer 8-Bit-Zahl berechnet werden. Anstatt jedes Mal eine Multiplikationsroutine aufzurufen, kann gleich zu Beginn einmalig eine Tabelle aller möglichen Ergebnisse von 0 * 0 bis 255 * 255 erstellt werden. Da die Resulate zwei Bytes in Anspruch nehmen, kostet dieses Vorhaben insgesamt 512 Byte Arbeitsspeicher.
Listing 1 zeigt die praktische Vorgehensweise und stellt zwei Routinen zur Verfügung, die Sie in eigenen Programmen einsetzen können. INITSQ wird am Anfang des Programms aufgerufen und dient dazu, die Tabellenwerte zu berechnen. Die Routine kommt ohne Multiplikationen aus, da sie eine spezielle Eigenart der Quadratzahlen ausnutzt: Der Wert n2 entspricht exakt der Summe der ersten n ungeraden Zahlen!
Hier ein paar Beispiele:
1 = 1
4 = 1 + 3
9 = 1 + 3 + 5
16 = 1 + 3 + 5 + 7
...und so weiter. Die Assemblerschleife muß also nur die ungerade Zahl laufend um 2 erhöhen,zur bisherigen Summe addieren und das Resultat in die Tabelle eintragen.
Die Routine GETSQ ist für die eigentliche Quadratermittlung zuständig. Sie braucht nur noch aus dem Faktor im A-Register die Tabellenadresse zu berechnen und den dort befindlichen Wert ins DE-Registerpaar einzulösen. Einschließlich RET-Befehl begnügt sie sich mit nur 84 Taktzyklen (21 MikroSekunden) und ist damit wesentlich schneller als jede herkömmliche Multiplikationsroutine.
Die Winkelfunktionen kommen häufig bei der Berechnung von Grafikkoordinaten zum Einsatz, die abschließend auf Ganzzahlen gerundet werden. Die hohe Genauigkeit der Fließkomma-Arithmetik ist in diesem Fall jedoch überhaupt nicht erforderlich — wie läßt sich dieser unnötige Aufwand vermeiden? Ideal wäre eine Routine, die von vornherein nur mit Integer-Ganzzahlen im Bereich von -32768 bis +32767 rechnet. Die damit verbundene vierstellige Genauigkeit ist zumindest für Grafikkordinaten absolut ausreichend. Leider wollen sich die Sinus- und Cosinusfunktion mit ihrem Weitebereich von -1...+ 1 und dem langen "Nachkommaschwanz" nicht in dieses System einfügen.
Das Problem läßt sich jedoch umgehen, wenn man die Sinuswerte mit einem konstanten Faktor multipliziert und auf Ganzzahlen rundet. Als Faktor bietet sich die größte positive Integerzahl (32767) an, um eine maximale Genauigkeit zu erzielen. Beschränkt man sich zusätzlich auf ganzzahlige Funktionsargumente, lassen sich die dazugehörigen Sinuswerte leicht in einer Tabelle erfassen. Ein Blick auf Listing 2 zeigt die praktische Realisierung: Ab Zeile 1090 finden Sie 91 Ganzzahl-Sinuswerte für die Winkel von 0 bis 90 Grad. Die Tabelle wird mit Hilfe der Assemblerdirektive DW (Define Word) in das Programm integriert. Für die Berechnung des Integer-Sinus ist die Routine INTSIN ab Zeile 1360 zuständig. Sie nimmt im HL-Registerpaar den Winkel in Empfang, ermittelt den passenden Tabellenwert und gibt das Ergebnis in HL zurück. Die Zeilen 1360 bis 1450 beschäftigen sich damit, den Winkel in den Bereich von 0...360 Grad zu bringen. Eine Addition oder Subtraktion von 360 hat keinen Einfluß auf das Resultat, da sich die Sinusfunktion nach jeweils 360 Grad periodisch wiederholt. Ab Zeile 1470 sorgt eine Fallunterscheidung für das richtige Ergebnis bei Winkeln, die größer als 90 Grad sind, indem sie die Symmetrie der Sinusfunktion ausnutzt. Um das Assembler-Chaos etwas zu entwirren, hier der analoge Vorgang in BASIC (w ist der Winkel, t die Tabelle als Array DIM(90), e das Ergebnis):
IF w < 90 THEN e = t(w)
ELSE IF w < 180 THEN e = t(180-w)
ELSE IF w < 270 THEN e = -t(w-180)
ELSE e = -t(360-w)
Als Resultat können also auch Werte < 0 auftreten. Die Routine GETSIN gibt jedoch grundsätzlich positive Werte zurück und signalisiert das Vorzeichen separat im A-Register (0 = positiv, &FF = negativ). Auf die übliche Zweierkomplement-Codierung für vorzeichenbehaftete Ganzzahlen wird hier bewußt verzichtet, um Probleme mit der Multiplikationsroutine zu vermeiden, die gleich noch ins Spiel kommt.

Schnelle Kurven

Wir verfügen also jetzt über eine Routine, die aus einem beliebigen ganzahligen Winkel einen mit 32767 multiplizierten Integer-Sinus ermittelt. Was läßt sich damit anfangen? Angenommen, ein Programm soll die Sinuskurve auf dem Bildschirm zeichnen. Der Winkel wird entlang der x-Achse aufgetragen. Da sich die y-Achse des CPC-Grafikbildschirms von 0 bis 399 erstreckt, wird man den Sinuswert mit einem Vergrößerungsfaktor multiplizieren. Dazu kommt noch ein Summand für die Zentrierung der Kurve:
y = 150 * SIN(x) + 200
Wenn wir die Ganzzahl-Funktion INTSIN verwenden, ist eine Division durch 32767 erforderlich, um ein korrektes Resultat zu erhalten:
y = (150 * INTSIN(x))/32767 + 200
Dabei tauchen einige Probleme auf:
- Das Produkt 150 * INTSIN(x) setzt sich aus zwei 16-Bit-Werten zusammen und führt zu einem 32-Bit-Resultat (Long Integer). Wir brauchen also eine entsprechende Multiplikationsroutine.
- Die Division eines 32-Bit-Wertes durch 32767 ist programmtechnisch sehr aufwendig. Wesentlich einfacher ist eine Division durch die Zweierpotenz 32768 (215), die sich mit einfachen Schiebeoperationen erledigen läßt.
- Um das Ergebnis so genau wie möglich zu halten, sollte korrekt gerundet werden.
Die 32-Bit-Multiplikation finden Sie ab Zeile 1850. Die Kommentare geben Hinweise auf die Funktionsweise; eine ausführliche Erläuterung würde hier jedoch zu weit führen. Betrachten Sie die Routine MUL32 einfach als Black Box, die in HL und DE zwei positive 16-Bit-Faktoren in Empfang nimmt. Die ersten 16 Bit (Highword) werden in HL und die zweiten 16 Bit (Lowword) in DE zurückgegeben.
Die kleine Ungenauigkeit bei der Division (32768 statt 32767) nehmen wir in Kauf, da sie praktisch kaum Auswirkungen hat. Bleibt also noch die Rundung des Ergebnisses: Eine normale Ganzzahldivision (im CPC-BASIC der umgekehrte Schrägstrich) ist unbefriedigend, da sie den Divisionsrest einfach unterschlägt. 39 \ 10 ergibt zum Beispiel 3 statt 4. Es gibt jedoch eine simple Technik, die hier sofort für Abhilfe sorgt. Man muß nur vor der Division die Hälfte des Divisors addieren, also zum Beispiel (39 + 5) \ 10 = 4.
Die Grundlage für die Berechnung von Ausdrücken der Form f * SIN(x) ist also die folgende Formel:
y = (f * INTSIN(x) + 16384) \ 32768
Die Division durch 32768 (215) entspricht einem Rechtsschieben um 15 Bit. Noch einfacher ist es, den Dividenden erst mit 2 zu multiplizieren (das heißt, ein Bit nach links zu schieben) und dann durch 216 zu teilen. Das Rechtsschieben um 16 Bit können wir uns nämlich komplett sparen, indem wir einfach die oberen 16 Bit des 32-Bit-Wertes als Endergebnis übernehmen!

Sinus mit Faktor

Mit diesen Grundlagen sollte die Funktionsweise der Routine SINMUL ab Zeile 2400 einigermaßen verständlich sein. Sie berechnet einen Ausdruck der Form f * SIN(x) mit korrekter Ganzzahl-Rundung. Der Faktor f (der immer ≥ 0 sein muß) wird in DE und der Winkel x in HL übergeben; das Ergebnis ist im HL-Registerpaar zu finden, und zwar als vorzeichenbehaftete Integerzahl. Damit der Rundungs-Trick korrekt funktioniert, wird die Umwandlung in die Zweierkomplementdarstellung bereits vor der Division durchgeführt, wenn das A-Register ein negatives Vorzeichen signalisiert (siehe Zeile 2440, 2450). Diese Aufgabe erledigt das Unterprogramm VZW32 (32-Bit-Vorzeichenwechsel) ab Zeile 2110, indem es den Wert einfach byteweise von 0 abzieht.
Zum Abschluß bleibt noch die Frage, was zu tun ist, wenn nicht ein Sinus-, sondern ein Cosinuswert berechnet werden soll. Nichts einfacher als das! Nach den Rechengesetzen für die Winkelfunktionen gilt nämlich
COS(x) = SIN(x + 90°)
Entsprechend einfach fällt die Routine COSMUL ab Zeile 2310 aus, die einen Ausdruck der Form f * COS(x) berechnet. Sie erhöht nur den Winkel in HL um 90 und überläßt alles weitere der Routine SINMUL.

Hier die BASIC-Version:
10 MODE1:DEG
20 FOR x=0 TO 638 STEP 2
30 PLOT x,150*SIN(x)+200
40 PLOT x,150*COS(x)+200
50 NEXT x
Während die BASIC-Kurven mehr oder weniger über den Bildschirm schleichen, sorgt die reine Integer-Arithmetik in Maschinensprache für einen rasanten Bildaufbau. Assemblieren Sie das Listing, und überzeugen Sie sich selbst!
(Matthias Uphoff/cd)
Hier als Z80-Quelle.
Hier als Z80-Quelle.

Für den PCW funktioniert die Testroutine zum Zeichnen einer Sinus- und Cosinuskurve auch, wenn entsprechende Änderungen vorgenommen werden:
  1. Bildschirmtreiber: Der Plot-Befehl aus dem CPC-Beispiel existiert beim PCW nicht. Ich habe hierzu die entsprechende Routine aus dem Artikel Turbo-Grafik für den JOYCE verwendet.
  2. Bildschirmauflösung: Beim CPC beträgt die Anzahl der Vertikalen Zeilen 400, beim PCW 256. Die Konstante YMAX im Programm passt den Unterschied an.
  3. Laufzeitanpassung: Im Programm für den PCW wird
    • der Code zum plotten in das COMMON Memory kopiert,
    • der Bildschirm gelöscht,
    • Sinus und Cosinus dargestellt,
    • nach Ferigstellung auf eine Tastatureingabe gewartet und
    • der Bildschirm erneut gelöscht.
Hier das angepasste Programm und das Ergebnis, das sich durchaus sehen lassen kann!

Als Zusatz zu den vorgestellten Integer-Routinen hier noch zwei Unterprogramme zur Berechnung der Quadratwurzel, abgedruckt in der PC AMSTRAD INTERNATIONAL 7/87:
Zitat: Wurzeln im Assembler zu berechnen, gehört sicherlich mit zu den schwierigsten Aufgaben in Assembler, noch dazu, wenn das Programm extrem kurz sein soll.

Während die erste Routine das Ergebnis im Akkumulator berechnet, wird in der zweiten Routine das Ergebnis im Register DE berechnet. In der ersten Routine können nur Quadratwurzeln im Bereich von 0 bis 65025 berechnet werden, da das Ergebnis in einem Byte (0-255) gespeichert wird: 2552 = 65025!
Hier als Z80-Quelle. Hier als Z80-Quelle.

Unter z80-heaven.wikidot.com finden sich u.a. folgende Unterprogramme zur Berechnung von Wurzeln mit Integer-Zahlen.
Ergibt die Wurzel von E, zur nächsten Ganzzahl gerundet: Ergibt die Wurzel von E (nach unten gerundet).
Hier als Z80-Quelle. Hier als Z80-Quelle.
Ergibt die Wurzel von HL (nach unten gerundet): Ergibt die Wurzel von L, nach unten gerundet:
Hier als Z80-Quelle. Hier als Z80-Quelle

Eingescanned von Werner Cirsovius
April 2006, Januar 2014
© CPC