»Lange Zahlen« in Pascal

Wer ärgert sich nicht darüber: Da hat man nun einen Rechner für ein paar tausend Mark - einen wahren Datenfresser vor sich, lädt Pascal hinein und was passiert beim jonglieren mit großen Zahlen?
»Overflow!«
Der Rechner streikt!

Wer will es genau wissen?

444! =

215849822966416938124093484678763977522
688933423911139272822344760723304263104
294496589577688721302429607721233905499
865884989120401810449286699449594741700
138829359308216743908805909331418696422
726219252910161110657569776306015863272
974081788704178176125928967124961489694
143231665834254132403506160786008478880
703187870540711394181810496944119553978
895775427586762396443695195620412156896
954424431104417491076774630170333975915
684392008987506995167749373363243965364
378536133789289474389082982363959239343
809192095847737811603628491131611260552
556744347257254652901445750019966944662
656339166877115737365363234874392035875
855364646100671808307470789777194448586
202523372013451343806787571920976593756
808399474310410313111150943567072536841
149261470999161741165789270579980199429
220358271673522887660057254252337382994
128874850457319557218710438670021146981
188082426807373004500000000000000000000
000000000000000000000000000000000000000
00000000000000000000000000000000000000
0000000000

Bis 100 Stellen rechnet LONGINT.PAS ja
noch recht flott, doch mit mehr...
Aber man sage nichts gegen die
Genauigkeit!

Spätestens jetzt wird Snobisten klar: Ein Taschenrechner für 50 Mark kann besser rechnen! Während ein Taschenrechner noch mühelos 69! = 1,711224 * 10⁹⁸ berechnet, geben die meisten Programmiersprachen der Computer nur noch die Fehlermeldung »Overflow« aus.
Ganz schlimm wird es, will man das exakte Ergebnis bekommen. Schafft ein Taschenrechner hier nur 11!, so muß der PC schon bei 8! die Flagge streichen.
Abhilfe schafft hier »LONGINT.PAS« - eine Unterprogrammbibliothek, die die Operationen +,-,*,=,<, »shl« und »shr« für einen von Ihnen zu bestimmenden, neuen Datentyp »LONGINTEGER« beinhaltet.
Maschinenspracheprogrammierer wird die Vorgehensweise nicht unbekannt sein: Die Zahlen werden in einem »Array of Byte« dargestellt, auf welches dann die Dualarithmetik angewandt wird. Für Pascalisten, deren Compiler den Datentyp »Byte« nicht kennt, sei vermerkt, daß es ein »Array of Char« in Verbindung mit den Funktionen »CHR« und »ORD« auch tut.
Allerdings leidet darunter etwas die Geschwindigkeit!
Die Aufrufe der,für den Typ LONGINTEGER zuständigen Operationen sind etwas gewöhnungsbedürftig. Da ein Pascal-Compiler vom Benutzer nicht ohne weiteres erweitert werden kann, wurden die LONGINTEGER-Operationen als Prozeduren realisiert.
Die Tabelle zeigt eine Gegenüberstellung der »normalen« Integer-Operationen und der für den LONGINTEGER-Typ.
Wie groß und genau Sie mit Ihrem neuen Typ rechnen können, bestimmt die Konstante »LONGINT_Len«, mit der die Länge des die langen Zahlen aufnehmenden Byte Arrays bestimmt wird. Diese ergibt sich aus der nächstgrößeren Zahl größer oder gleich »LONTGINT_Digits/logm(256)«. Als Beispiel für die Leistungsfähigkeit der Bibliothek wurde ein kurzes Progrämmchen angefügt, das Fakultäten berechnet.

(Dietmar Bückart)

Tabelle:

Turbo-Anweisung	Longinteger	Anmerkungen
z:=x + y	LAdd(z, x, y)	K. Überlaufprüfung
z:=x - y	LSub(z, x, y)	K. Überlaufprüfung
z:=x * y	LMul(z, x, y)	K. Überlaufprüfung
z:=x div y	LDiv(z, r, x, y)
r:=x mod y
x < y		LLess(x, y)	bool'sche Funktion
x:=x shl 1	LShl(x)		schneller als x * 2
x:=x shr 1	LShr(x)		schneller als x Div 2
readln(x)	LRead(x)	Eingabe
write(x)	LWrite(x)	Ausgabe, langsam
---		IntToLong(x, i)	Wandelt Integer i in Langinteger x um

Listing von LONGINT.PAS

(* -------------------------------------------------- *)
(*                   LONGINT.PAS                      *)
(* impl. beliebig langer positiver Integers in Pascal *)
(* -------------------------------------------------- *)

CONST LONGINT_Digits= 150; (* Anz. Stellen f. LONGINT *)
      (* Anzahl Bytes zur internen Darstellung der
         gewuenschten Stellenzahl
         = LONGINT_Digits/log(256):                   *)
      LONGINT_Len   = 63;

TYPE  LONGINT = ARRAY[1..LONGINT_Len] OF BYTE;

(* -------------------------------------------------- *)
(*  norm. INTEGER in eine LONGINTEGER konvertieren:   *)

PROCEDURE IntToLong (VAR a: LONGINT; b: INTEGER);

VAR i: INTEGER;

BEGIN
  FOR i := 3 TO LONGINT_Len DO a[I]:=0;
  a[2] := HI(b);  a[1] := LO(b);
END;

(* -------------------------------------------------- *)
(*                    a < b ?                         *)

FUNCTION LLess (a, b: LONGINT): BOOLEAN;

VAR i: INTEGER;

BEGIN
  i := LONGINT_Len; LLess := FALSE;
  WHILE (a[i] = b[i]) AND (i > 1) DO i := Pred(i);
  IF a[i] < b[i] THEN LLess := TRUE;
END;

(* -------------------------------------------------- *)
(*                  sum := a + b                      *)

PROCEDURE LAdd (VAR sum: LONGINT; a, b: LONGINT);

VAR i, j: INTEGER;
    mov : BYTE;

BEGIN
  mov := 0;
  FOR i := 1 TO LONGINT_Len DO
  BEGIN
    j := a[i] + b[i] + mov;  sum[i] := LO(j);
    mov := HI(j);
  END;
END;

(* -------------------------------------------------- *)
(*                 diff := a - b                      *)

PROCEDURE LSub (VAR diff: LONGINT; a, b: LONGINT);

VAR i: INTEGER;
    one: LONGINT;

BEGIN
  (* Zweierkomplement des Subtrahenden b bilden: *)
  FOR i := 1 TO LONGINT_Len DO b[i] := NOT b[i];
  IntToLong(one, 1);  LAdd(b, b, one);
  LAdd(diff, a, b);        (* und zu a addieren! *)
END;

(* -------------------------------------------------- *)
(*    Links-Shift, entspricht Multiplikation mit 2    *)

PROCEDURE LShl (VAR a:LONGINT);

VAR i: INTEGER;
    mov, schieb: BYTE;

BEGIN
  mov := 0;
  FOR i := 1 TO LONGINT_Len DO
  BEGIN
    schieb := (a[i] SHL 1) + mov;
    mov := a[i] SHR 7; a[i] := schieb;
  END;
END;

(* -------------------------------------------------- *)
(*     Rechts-Shift, entspricht Division durch 2      *)

PROCEDURE LShr (VAR a: LONGINT);

VAR i: INTEGER;
    mov, schieb: BYTE;
BEGIN
  mov := 0;
  FOR i := LONGINT_Len DOWNTO 1 DO
  BEGIN
    schieb := (a[i] SHR 1) + mov;
    mov := a[i] SHL 7; a[i] := schieb;
  END;
END;

(* -------------------------------------------------- *)
(*                 prod := a * b                      *)

PROCEDURE LMul (VAR prod: LONGINT; a,b: LONGINT);

VAR i,j,x, alen, blen: INTEGER;

BEGIN
  IntToLong(prod,0);
  alen := LONGINT_Len;  blen := LONGINT_Len;
  WHILE (a[alen] = 0) AND (alen > 1) DO
    alen := Pred(alen);
  WHILE (b[blen] = 0) AND (blen > 1) DO
    blen := Pred(blen);
  FOR i := 1 TO alen DO
  BEGIN
    x := 0;
    FOR j := 1 TO blen DO
    BEGIN
      x := a[i]*b[j] + HI(x) + prod[pred(i+j)];
      prod[pred(i+j)] := LO(x);
    END;
    prod[i+blen] := HI(x);
  END;
END;

(* -------------------------------------------------- *)
(*         quot := a DIV b, rest := a MOD b           *)

PROCEDURE LDiv (VAR quot, rest :LONGINT; a, b: LONGINT);

VAR p, one: LONGINT;

BEGIN
  p := b;  IntToLong(quot, 0);  IntToLong(one,1);
  REPEAT LShl(b) UNTIL NOT LLess(b, a);
  REPEAT
    LShl(quot);  LShl(b);
    IF NOT LLess(a ,b) THEN
    BEGIN
      LSub(a, a, b);  LAdd(quot, quot, one);
    END;
  UNTIL NOT LLess(p, b);
  rest := a;
END;

(* -------------------------------------------------- *)

PROCEDURE LRead (VAR a: LONGINT);

VAR i, j, m: INTEGER;
    lin: STRING[LONGINT_Digits];
    mov: BYTE;

BEGIN
  IntToLong(a, 0);  ReadLn(lin);  i := 1;
  WHILE i <= Length(lin) DO
  BEGIN
    IF lin[i] IN ['0'..'9'] THEN
      BEGIN
        mov := Ord(lin[i]) - Ord('0');
        FOR j := 1 to LONGINT_Len DO
        BEGIN
          m := a[j] * 10 + mov;
          a[j] := LO(m);  mov := HI(m);
        END;
        i := Succ(i);
      END
    ELSE
      i := Succ(Length(lin));
  END;
END;

(* -------------------------------------------------- *)

PROCEDURE LWrite(a: LONGINT);

VAR i, j, k, m, mov: INTEGER;
    puffer: ARRAY[1..LONGINT_Digits] OF CHAR;

BEGIN
  j := LONGINT_Len;  k := 0;
  WHILE (a[j]=0) DO j := Pred(j);
  REPEAT
    k := Succ(k);   mov := 0;
    FOR i := j DOWNTO 1 DO
    BEGIN
      m := a[i] + mov * 256;
      a[i] := m div 10;  mov := m MOD 10;
    END;
    puffer[k] := CHAR(Ord('0') +mov);
    IF a[j] = 0 THEN j := Pred(j);
  UNTIL j = 0;
  FOR i := k DOWNTO 1 DO Write(puffer[i]);
END;

(* -------------------------------------------------- *)
(*              Ende von LONGINT.PAS                  *)

Beispielprogramm FAK.PAS

Im Magazin „PASCAL", Ausgabe 8'87, wurden zum obigen Artikel folgende Korrekturen abgedruckt:

»Bugs« bei langen Zahlen

In dem Artikel »Lange Zahlen« (PASCAL 6/87, S.81f) ist die Prozedur »LDiv« nicht korrekt. Bei der Division von Zweier-Potenzen treten Fehler auf. Demnach wäre für a = 8 und b = 2 der Quotient 3 mit Rest 2.
Eine korrigierte Fassung von »LDiv« finden Sie in Listing 3. In Listing 4 ist die Funktion »LGleich« abgedruckt, mit der lange Integers auf Gleichheit überprüft werden können. Eine Anwendung der Funktionen sehen Sie in dem Programm Primzahl (Listing 5), wo man mit Hilfe dieser Bibliothek sehr große Primzahlen berechnen kann.

Martin Ost
Karlsruhe 1

Listing 3

Listing 4

PROGRAM Primzahl(Input,Output);
  (* Demonstriert  die  LONG-INTEGER-LIBRARY  *)
  (* mit einem Primzahlen-Sieb.	              *)

(*$I LONGINT.LIB *)

VAR n, Null, Zwei, Rest, Quot, z_a, z_b, Quad : LONGINT;
    Teiler_gefunden : BOOLEAN;

BEGIN
  ClrScr;
  WriteLn(' n >= 10 ');
  Write(' n = '); LRead(n), WriteLn;
  Write(' Ausgabe aller Primzahlen <= ');
  LWrite(n); WriteLn; WriteLn;
  Write('2  3  5  7  ');
  IntToLong(Null,0);
  IntToLong(Zwei,2);
  IntToLong(z_a,7);
  IntToLong(z_b,1);
  REPEAT
    LAdd(z_a,z_a,Zwei);
    REPEAT
      LAdd(z_b,z_b,Zwei);
      LMul(Quad,z_b,z_b);
      LDiv(Quot,Rest,z_a,z_b);
      Teiler_gefunden := LGleich(Rest,Null)
    UNTIL Teiler_gefunden OR (NOT LLess(Quad,z_a));
    IF NOT Teiler_gefunden THEN BEGIN
      Write('  ');
      LWrite(z_a)
    ENDs
    IntToLong(z_b,1)
  UNTIL NOT LLess(z_a,n)
END.

Listing 5