Arithmetik mit 64-Bit-Zahlen
Der Datentyp long in also 32-Bit-Variablen, umfasst einen Bereich von ca. 4,3 Milliarden Zahlen. Das sieht auf den ersten Blick unermesslich viel aus, bei bestimmten Aufgabenstellungen aber reicht selbst das nicht. Deshalb wird hier ein Routinenpaket zum Umgang mit 64-Bit-Zahlen vorgestellt. Dieser Typ dublong hat einen Wertebereich von ca. 18,5 Trillionen Zahlen!
Sie werden sich bestimmt fragen, wie man es eigentlich anstellt, an die Grenzen des long-Typs zu stoßen. Aber es geht schneller, als man denkt. Stellen Sie sich folgende Berechnung vor: Es soll ermittelt werden, wieviel Prozent x Mark von y Mark sind. Um Rundungsfehler zu vermeiden, verwendet man nicht float-Variablen, sondern longs. Da man so keine Nachkommastellen darstellen kann, werden die Beträge in Pfennig angegeben. Das Ergebnis soll auf drei Stellen nach dem Komma genau sein, d.h. ist mit 1000 multipliziert (12.345%->12345). Die Multiplikationen müssen vor der Division erfolgen, damit sich keine Abschneidefehler ergeben. Die Formel zur Berechnung lautet dann:
erg = x * 100 * 1000/y
Nehmen wir jetzt als Zahlenbeispiel x = 300 DM, in unserer Darstellung 30000. Der größte Zahlenwert, der in der Berechnung erreicht wird, ist dann:
30000 * 100 * 1000 = 3e9
Das paßt schon nicht mehr in eine long-Variable, denn die faßt nur Werte bis +2.147e9, es ergibt sich somit ein Überlaufsfehler. Schon bei diesen relativ kleinen Beispielzahlen wird also der Wertebereich von 32 Bit „gesprengt“.
Grundlegende Operationen mit dublongs
Um einen neuen Datentyp dublong für 64-Bit-Variablen verwenden zu können, brauchen wir als erstes eine Struktur dafür. Dazu nimmt man einfach 2 longs, wobei aber die Low-Hälfte ein unsigned long sein muß, damit Vergleiche richtig funktionieren.
Beachtenswert ist, daß dub-long-Strukturen direkt als Parameter übergeben werden können, nicht nur Zeiger auf sie. Die meisten C-Compiler behandeln das so, als würden da getrennt 2 longs stehen. Ihre Reihenfolge ist genauso, wie sie in der Struktur stehen, d.h. zuerst High, dann Low. Im Paket wird das auch ab und zu umgekehrt und 2 longs statt einem dublong als Parameter verwendet. Dadurch kann man sich oft dl_assign() sparen.
Nicht ganz so einfach ist es beim Rückgabewert. Ist das Ergebnis einer Funktion ein dublong, kann es nicht direkt mit return() zurückgegeben werden. Folglich muß der Funktion zusätzlich ein Zeiger mitgegeben werden, der ihr sagt, wohin das Resultat geschrieben werden soll. Zur Vereinfachung von Rechenausdrücken wird dieser Zeiger (freilich unverändert) als Funktionswert zurückgeliefert, so daß im selben Term mit dem Ergebnis weitergerechnet werden kann. Ein Beispiel wird klarer machen, was ich meine:
erg = dl_div(*dl_add(dub1, dub2, &dub3 ), dub4, &rest);
entspricht (in „normalen“ Ty-pen):
erg = (n1+n2)/n4;
duh3 wird als Platz für das Ergebnis der Addition benötigt. Es muß jedoch von dl_div() nicht direkt angesprochen werden, so daß man die Operationen ineinanderschachteln kann. Das ist so nur möglich, wenn das Resultat von dladd() mit ‘*dl_add' angesprochen werden kann.
Als nächsten Schritt bei der Einführung eines neuen Typs muß man Funktionen bereitstellen, die die Umwandlungen zwischen dem neuen und bekannten Typen erledigen. Im dublong-Paket nehmen diese Aufgaben dl_todl() und dl_tolong() wahr. dl_todl() wandelt long in dublong, indem es den long-Wert als Low-Hälfte der dublong betrachtet. Die High-Hälfte ist 0, wenn die long-Zahl größer oder gleich Null ist, sonst -1. Dieser Vorgang, den man Vorzeichenerweiterung nennt, geschieht analog zu den EXT-Befehlen der 68000-Maschinensprache, die Byte auf Wort oder Wort auf Long erweitern. dl_tolong() gibt den Wert einer dublong als long zurück. Damit das funktionieren kann, muß der Wert natürlich in eine long-Variable passen, sonst ist das Ergebnis ohne Sinn.
Eine weitere einfache Operation ist das Negieren. Man könnte dafür zwar auch die Subtraktion benutzen ( -x = 0-x ), aber schneller geht es, wenn man beide Hälften mit dem ‘~*-Operator NOT-iert, d.h. alle Bits umdreht, und dann 1 addiert. Das ist genau die Konvention, die intern zur Darstellung von negativen Zahlen benutzt wird, das sog. Zweierkomplement.
Addition und Subtraktion
Addition und Subtraktion von dublongs gestalten sich recht einfach, da die benötigten Operationen bereits in der 68000-Maschinensprache enthalten sind. Es gibt neben den einfachen Addier- bzw. Subtrahierbefehlen ADD und SUB Varianten, die den Übertrag einer vorangegangenen Operation berücksichtigen: ADDX und SUBX. Diese dürfen aber nur auf zwei Datenregister angewandt werden.
Zum Addieren zweier dublongs addiert man also zuerst die Low-Hälften mit ADD.L, wobei ein Übertrag entstehen kann. Dann addiert man die High-Hälften mit ADDX.L, wobei der vorherige Übertrag berücksichtigt wird. Das Subtrahieren erfolgt entsprechend. Mit weiteren ADDX/SUBX könnte man auch beliebig lange Zahlen addieren bzw. subtrahieren.
Multiplikation
Etwas kniffliger wird es bei der Multiplikation. Die Funktion dl_mul() soll zwei longs zu einem dublong verknüpfen, entspricht also dem Maschinenbefehl MUL in doppelter Breite.
Um unsere Multiplikation durchzuführen, gehen wir vor wie jemand, der nur das kleine Einmaleins kennt und zweistellige Zahlen multiplizieren will (und keinen Taschenrechner besitzt): Man multipliziert die Stellen einzeln miteinander und addiert die Ergebnisse entsprechend ihrer Wertigkeit. Übertragen heißt das: das kleine Einmaleins sind die 16-Bit-Zahlen, deren Produkte mit MUL berechenbar sind. Die Analogie sehen Sie in Bild 1.
Um uns das Leben etwas einfacher zu machen, gehen wir von positiven Faktoren aus. Das heißt, wir ermitteln, welches Vorzeichen das Ergebnis haben wird, und machen die Faktoren einfach positiv. Am Ende der Multiplikation muß dann eventuell das Ergebnis noch negativ gemacht werden.
Um das Addieren der Teilergebnisse etwas cleverer zu gestalten, gehe ich dabei folgendermaßen vor: Zuerst wird High1 * High2 und Low 1 * Low2 berechnet und daraus ein dublong gebildet. Dann müssen noch High1 * Low2 und High2 * Low1 zu der Mitte des Ergebnisses addiert werden. Dabei kann aber jeweils ein Übertrag in das höchstwertigste Wort entstehen. Um diesen zu addieren, holen wir das oberste Wort in ein Register (D1) und addieren mit ADDX den Wert 0. Das bewirkt, daß Dl um 1 erhöht wird, wenn aus der vorangegangenen Addition ein Übertrag entstanden ist, sonst bleibt es unverändert. Da bei ADDX nur Register erlaubt sind, benutze ich D2 als Nullregister. Zum Schluß muß D1 natürlich noch zurückgeschrieben werden.
Da dl_mul() nur long mal long rechnen kann, gibt es noch eine zweite Multiplikationsroutine im Paket: dl_mul2(). Damit können zwei dublongs verknüpft werden, jedoch hat das Ergebnis nicht immer in einem dublong Platz. Auch z.B. das Resultat von int mal int paßt ja nicht immer in eine int-Variable. Man muß also hier aufpassen, daß man keine zu großen Zahlen nimmt. Das liegt -ganz C-Philosophie - in der Verantwortung des Programmierers.
Doch zur Arbeitsweise von dl_mul2(): Hier kann man das ganze leider nicht so einfach auf ein paar MULs zurückführen. Es müßten nämlich sage und schreibe 16 MUL-Befehle sein, dazu kommen noch die Additionen und Übertragsrechnungen. Schneller geht es, indem man die Multiplikation auf Addition und Schieben zurückführt. Das Prinzip ist dasselbe wie oben, nur daß jetzt das Einmaleins recht klein ist: entweder addieren wir den ersten Operator zum Ergebnis oder nicht. Ob das passiert, hängt vom höchstwertigsten Bit des 2. Operators ab. Nach jeder Abfrage schieben wir diesen nach links, so daß das nächste Bit ganz links steht. Zusätzlich nehmen wir das Ergebnis mal 2, sprich schieben es auch nach links, und wiederholen das ganze für alle 64 Bits. Damit entsteht so etwas wie in Bild 2 zu sehen ist. Dort wird es zwar nur für 8 Bit dargestellt, aber so spart man sich lange Zahlenkolonnen, und der Vorgang wird trotzdem klar.
Die Division
Als die schwierigste Aufgabe erweist sich die Division. Sie kann nämlich nicht wie die Multiplikation auf die vorhandenen Maschinenbefehle zurückgeführt werden. Es muß also eine Methode eingesetzt werden, die die auf Subtraktionen zurückführt. Die Idee, wie vorher ein ähnliches Prinzip wie die schriftliche Division anzuwenden, erweist sich leider als recht ineffizient.
Die schnellste Lösung ergab ein Algorithmus aus [1], den ich auf 64-Bit-Zahlen und 68000er-Assembler anpaßte. Er arbeitet nach folgendem Schema: Von der High-Hälfte des Dividenden wird der Divisor abgezogen. Nach jeder Subtraktion muß das Übertrags-Flag X negiert werden. Dieses Bit wird dann von rechts in das Ergebnis hineingeschoben, und auch der Dividend wird nach links ge-shiftet. Abhängig vom X-Flag wird der Divisor entweder addiert oder subtrahiert. Dieser Vorgang wird 32mal wiederholt und zum Schluß - wenn nötig - der Rest positiv gemacht. Das Prinzip dieses Algorithmus’ ist in Bild 3 dargestellt.
Das ist - offen gesagt - etwas kompliziert, auch die mathematischen Grundlagen sind nicht ohne weiteres offensichtlich, aber es funktioniert.
Umwandlung ASCII -> dublong
Als kleines Extra ist im dublong-Paket noch eine Funktion atodl() enthalten, die ähnlich atoi() oder atol() einen String in eine dublong-Variable liest. Sie funktioniert nach dem gleichen Algorithmus wie alle diese Routinen: Eine Ziffer wird gelesen und zum bisherigen Wert addiert. Bevor die nächste Zahl betrachtet wird, muß der aktuelle Wert noch mit 10 multipliziert werden. Abgebrochen wird der Vorgang, wenn ein Zeichen kommt, das keine Zahl ist. Falls am Anfang des Strings ein stand, muß das Gesamtergebnis noch negiert werden. Mit Hilfe dieser Prozedur kann man somit Konstanten verwenden, die nicht in longs passen, oder dublongs aus Benutzereingaben lesen.
Literatur:
[1J Zaks, Rodney: Programming the 6502, Sybex Inc., Berkeley 1980
/* Headerfile zu 'dublong.c' */
typedef struct {
long high;
unsigned long low;
} duhlong;
/* Parameter : */
extern dublong *dl_todl(); /* ( 1, &dl ) */
extern long dl_tolong(); /* ( dl ) */
extern dublong *dl neg(); /* ( &dl ) */
extern dublong *dl_add(); /* ( dl,dl,&dl ) */
extern dublong *dl_sub(); /* ( dl,dl,&dl ) */
extern dublong *dl_mul(); /* ( 1, 1, &dl ) */
extern dublong *dl_mul2(); /* ( dl,dl,&dl ) */
extern long dl_div(); /* ( dl, 1, &1 ) */
extern dublong *atodl(); /* ( str, &dl ) */
/* Implementierung der essentiellen Funktionen*/
/* für 'double longs’, d.h. vorzeichenbehaf- */
/* tete 64-Bit-Zahlen */
/* (c) 1991 MAXON Computer */
/* by Roman Hodek */
/* Funktionsübersicht: */
/* dl_todl() : erweitert long zu dublong */
/* dl_tolong() : verwandelt dublong in long; */
/* wenn der Zahlenbereich von */
/* long überschritten ist, ist */
/* das Ergebnis ohne Aussage */
/* dl_neg() : negiert eine dublong */
/* dl_add() : addiert : dl = dl+dl */
/* dl_sub() : subtrahiert : dl = dl-dl */
/* dl_mul() : multipliziert : dl = 1*1 */
/* dl_mul2() : multipliziert : dl = dl * dl */
/* dl_div() : dividiert l = dl/1 und */
/* l = dl%l */
/* atodl() : liest ein dublong aus einem */
/* String */
/* Alle Funktionen, die ihr dublong-Ergebnis */
/* über einen Zeiger ablegen, liefern diesen */
/* Zeiger auch als Rückgabewert. */
/* die dublong-Struktur : */
typedef struct {
long high,low;
} dublong;
/* Grenzwerte des Typs dublong : */
#define DLONG_BIT 64
#define DLONG_MIN -9.223372037e18
#define DLONG_MAX 9.223372037e18
#ifndef abs
#define abs(x) ( ((x)<0) ? -(x) : (x) )
#endif
/* erweitert 1 auf 64 Bit und weist es dl zu */
dublong *dl_todl( l, dl )
dublong *dl;
long l;
{ dl->low = 1;
dl->high = { 1 >=0 ) ? 0 : -1;
return( dl );
}
/* wandelt dublong ’dl' in long */
long dl_tolong( dl )
dublong dl;
{ return( dl.low ); }
/* negativiert ’op' ( = monadisches Minus ) */
dublong *dl_neg( op )
dublong *op;
{ dublong *dl_add();
/* beide Hälften negieren */
op->low = ~op->low; op->high = ~op->high;
/* + 1 wegen Zweierkon^lement */
dl_add( 0L, 1L, op->high, op->low, op );
return( op );
}
/* addiert die dublongs 1op1’ und ’op2’, Er- */
/* gebnis nach 'erg' */
/* *erg = op1 + op2 */
dublong *dl_add{ op1, op2, erg )
dublong op1, op2, *erg;
{ asm { move.l erg(A6),A0
move.l op1+4(A6),D0 ; Low's addieren
add.l op2+4<A6),D0
move.l D0,4(A0)
move.l op1(A6),D0 ; High's mit Über-
move.l op2(A6),D1 ; trag addieren
addx.1 D1,D0
move.l D0,(A0)
}
return( erg );
}
/* Subtrahiert dublong 'op2' von dublong 'op2'*/
/* Ergebnis nach 'erg' */
/* *erg = op1 - op2 */
dublong *dl_sub( op1, op2, erg )
dublong op1, op2, *erg;
{ asm { move.l erg(A6),A0
move.l op1+4(A6),D0 ; Low's subtra-
sub.l op2+4(A6),D0 ; hieren
move.l D0,4(A0)
move.l op1(A6),D0 ; High's mit Über-
move.l op2(A6),D1 ; trag subtrahie-
subx.l D1,D0 ; ren
move.l D0,(A0)
}
return( erg );
}
/* Multipliziert zwei longs zu einem dublong */
/* *erg = op1 * op2 */
dublong *dl_mul( op1, op2, erg )
long op1, op2;
dublong *erg;
{ register int negflag;
/* negflag zeigt an, ob das Ergebnis */
/* negativ ist */
negflag = ( (op1>0) ^ (op2>0) );
/* beide Faktoren positiv machen */
op1 = abs(op1); op2 = abs(op2);
/* Registerbelegung : */
/* d1 : höchstwertigstes Wort des Ergebnis- */
/* ses */
/* d2 : immer 0, zum addieren des Übertrags */
asm { move.l erg(A6),A0
clr.w D2
move.w op1+2(A6),D0 ; Low 1 *
mulu op2+2(A6),D0 ; Low 2
move.l D0,4(A0)
move.w op1(A6),D0 ; High 1 *
mulu op2(A6),D0 ; High 2
move.l D0,(A0)
move.w (A0),D1 ; höchstes Wort
move.w op1+2(A6),D0 ; Low 1 *
mulu op2(A6),D0 ; High 2
add.l D0,2(A0) ; addieren
addx.w D2,D1 ; Übertag add.
move.w op1(A6),D0 ; High 1 *
mulu op2+2(A6),D0 ; Low 2
add.l D0,2(A0) ; addieren
addx.w D2,D1 ; Übertrag add
move.w D1,(A0) ; höchstes Wort
} /* zurückschreiben */
/* Vorzeichen des Ergebnisses berichtigen */
if (negflag) dl_neg( erg );
return( erg );
}
/* Multipliziert 2 dublong's zu einem dublong */
/* Überschreitet das Ergebnis den Zahlenbe- */
/* reich, ist das Ergebnis ohne Aussage ! */
/* *erg = op1 * op2 */
dublong *dl_mul2( op1, op2, erg )
dublong op1, op2, *erg;
{ register long work1, work2, counter;
/* 3 Registervariablen reservieren (d5-d7) */
register int negflag = 0;
/* negflag zeigt an, ob das Ergebnis nega- */
/* ti ist; beide Operatoren werden posi- */
/* tiv gemacht */
if ( op1.high<0 ) { dl_neg( &op1 );
negflag ^= 1; }
if ( op2.high<0 ) { dl_neg( &op2 );
negflag ^= 1; }
/* Registerbelegung ; */
/* d0/dl : Multiplikand */
/* d1/d2 : Multiplikator (wird geschoben) */
/* d5 : Bitzähler */
/* d6/d7 : Ergebnis */
asm{ move.l op1(A6),D0 ; Register laden
move.l op1+4(A6),D1
move.l op2(A6),D2
move.l op2+4(A6),D3
clr.l D6 ; Ergebnis auf 0 initialis.
clr.l D7
moveq #63,D5 ; 64 Bits mal
loop: lsl.l #1,D3 ; Multiplikator schieben
roxl.1 #1,D2
bcc no_add ; wenn herausgeschobenes
add.l D1,D7 ; Bit = 1, Multiplikand
addx.1 D0,D6 ; zum Eregbnis addieren
no_add: lsl.l #1,D7 ; Ergebnis schieben
roxl.l #1,D6 ; erg *= 2
dbf D5,loop
movea.l erg(A6),A0 ; Ergebnis ablegen
move.l D6,(A0)
move.l D7,4(A0)
}
/* Vorzeichen des Ergebnisses berichtigen */
if (negflag) dl_neg( erg );
return( erg );
}
/* dividiert dublong 'op1' durch long 'op2' */
/* Quotient in long 'erg', Rest in 'rest' */
long dl_div( op1, op2, rest )
long op2, *rest;
dublong op1;
{ register long divisor, save_ccr;
/* 2 Registervar. reservieren (d6 und d7) */
register int negflag;
/* Divison durch 0 -> Exception auslosen */
if (op2==0) asm{ divu #0,D0 }
/* negflag zeigt an, ob das Eregbnis nega- */
/* tiv wird */
negflag = ((op1.high>=0) ^ (op2>0));
/* beide Operanden positiv machen */
op2 = abs(op2);
if ( op1.high<0 ) dl_neg( &op1 );
/* Registerbelegung : */
/* d0/d1 : Dividend */
/* d2 : Bitzähler */
/* d3 : Ergebnis */
/* d6 : Zwischenspeicher für CCR */
/* d7 : Divisor */
asm{ move.l op1(A6),D0 ; Register laden
move.l op1+4(A6),D1
move.l op2(A6),divisor
moveq #31,D2
sub.l divisor,D0 ; Subtrk. vorab
move SR,save_ccr
bchg #4,save_ccr
loop:
move save_ccr,CCR ; X-Flag ins Erg.
roxl.l #1,D3 ; schieben
lsl.l #1,D1 ; Dividend rot-
roxl.l #1,D0 ; tieren
btst #4,save_ccr
beq _add
sub.l divisor,D0 ; wenn X = 1
move SR,save_ccr ; Dvsr von Dvnd
bchg #4,save_ccr ; subtrahieren —
bra _loop
_add: add.l divisor,D0 ; wenn X = 0
move SR,save_ccr ; addieren
_loop: dbf D2,loop
btst #4,save_ccr ; wenn X = 0
bne no_add ; den Rest posi-
add.l divisor,D0 ; tiv machen
bclr #4, save_ccr
no_add:
move save_ccr,CCR ; letztes Bit ins
roxi.1 #1,D3 ; Erg. schieben
move.l D3,divisor ; 'divisor' ist jetzt das Ergeb.
move.l rest(A6),A0
move.l D0,(A0)
}
/* Ergebnis mit richtigem Vorzeichen zu- */
/* rückgeben */
return( negflag ? -divisor : divisor );
}
dublong *atodl{ str, erg )
register char *str;
register dublong *erg;
{ register int negflag;
/* Vorzeichen da ? */
if (*str==’-’) { ++str; negflag = 1; }
/* Ergebnis auf 0 initialisieren */
dl_assign( 0L, erg );
while( *str>='0' && *str<='9' ) {
dl_add( *dl_mul2( erg->high, erg->low, 0L, 10L, erg ),
0L, (long)(*str-'0'), erg );
++str;
};
if (negflag) dl_neg( erg );
return( erg );
}
Roman Hodek