Robert Sedgewick: Algorithmen

[ Inhaltsverzeichnis ] [ vorhergehende Seite ] [ nächste Seite ] [ Stichwort ]


35. Zufallszahlen



Methode der linearen Kongruenz

Das bekannteste Verfahren für die Erzeugung von Zufallszahlen, welches seit seiner Vorstellung durch D. Lehmer im Jahre 1951 fast ausschließlich verwendet wird, ist die sogenannte Methode der linearen Kongruenz (Restmethode von Lehmer). Falls seed eine beliebige Zahl enthält, so erstellt die folgende Anweisung unter Benutzung dieser Methode ein Feld mit N Zufallszahlen:

    a[0]:=seed;
    for i:=1 to N do
      a[i]:=(a[i-1]*b+1) mod m;

Das heißt, um eine neue Zufallszahl zu erhalten, nehme man die vorangegangene, multipliziere sie mit einer Konstanten b, addiere 1 und berechne den Rest, der sich bei der Division durch eine zweite Konstante m ergibt. Das Ergebnis ist stets eine ganze Zahl zwischen 0 und m - 1. Dies ist für die Anwendung in Computern günstig, da sich die Funktion mod gewöhnlich leicht implementieren läßt: Wenn wir den Überlauf bei den Rechenoperationen unterdrücken, so entfernt die Hardware der meisten Computer die Überlauf-Bits, und führt somit effizient eine mod-Operation aus, wobei m um eins größer ist als die größte ganze Zahl, die in einem Computerwort dargestellt werden kann. Auch diesmal sind die Zahlen nicht wirklich zufällig; das Programm erzeugt lediglich Zahlen, von denen wir hoffen, daß sie einem bestimmten anderen Prozeß als zufällig erscheinen.

So einfach die Methode zu sein scheint, war der Zufallszahlengenerator der linearen Kongruenz doch Gegenstand einer Bücher füllenden, ausführlichen und schwierigen mathematischen Analyse. Aus diesen Ergebnissen können wir gewisse Anhaltspunkte für die Wahl der Konstanten seed, b und m erhalten. Einige dem »gesunden Menschenverstand« entsprechende Prinzipien kommen hierbei zur Anwendung, doch in diesem Falle ist der gesunde Menschenverstand nicht ausreichend, um gute Zufallszahlen zu gewährleisten. Erstens sollte m groß sein; es kann die Größe des Computerwortes sein, muß jedoch nicht ganz so groß sein, wenn dies unzweckmäßig ist (siehe untenstehende Implementation). Normalerweise wird es zweckmäßig sein, m als eine Zweier- oder Zehnerpotenz zu wählen. Zweitens sollte b weder zu groß noch zu klein sein; eine günstige Wahl ist die Verwendung einer Zahl, die eine Ziffer weniger hat als m. Drittens sollte b eine beliebige Konstante ohne besonderes Muster in ihren Ziffern sein, abgesehen davon, daß sie auf ...x 21 enden sollte, wobei x gerade sein soll: diese letzte Forderung ist sicherlich eigenartig, doch sie verhindert das Auftreten bestimmter, möglicherweise problematischer Fälle, die durch die mathematische Analyse gefunden wurden.

Die oben beschriebenen Regeln wurden von D. E. Knuth entwickelt, in dessen Lehrbuch diese Fragen recht gründlich behandelt werden. Knuth zeigt, daß eine solche Auswahl bewirkt, daß die Methode der linearen Kongruenz gute Zufallszahlen erzeugt, die verschiedenen komplizierten statistischen Tests genügen. Das ernsthafteste potentielle Problem, das schnell sichtbar werden kann, besteht darin, daß der Generator in einen Zyklus geraten kann und viel früher, als es der Fall sein sollte, Zahlen erzeugt, die er bereits erzeugt hat. Zum Beispiel führt die Wahl von b = 19, m = 381 und seed = 0 zur Erzeugung der Folge 0, 1, 20, 0, 1, 20, ..., was eine nicht sehr zufällige Folge ganzer Zahlen zwischen 0 und 380 ist. Leider sind nicht alle derartigen Schwierigkeiten so leicht zu erkennen, so daß man gut beraten ist, wenn man den von Knuth angegebenen Richtlinien folgt, wodurch man viele der von ihm entdeckten versteckten Fallen umgeht.

Um den Zufallszahlengenerator zu starten, kann man einen beliebigen Anfangswert benutzen, ohne daß dies besondere Auswirkungen hätte (abgesehen natürlich davon, daß verschiedene Anfangswerte zur Erzeugung unterschiedlicher zufälliger Folgen führen). Oft ist es unnötig, wie im obigen Programm die gesamte Folge zu speichern. Stattdessen verwenden wir einfach eine globale Variable a, die mit einem bestimmten Wert initialisiert und dann durch die Berechnung a:=(a * b + 1) mod m aktualisiert wird.

In Pascal (und vielen anderen Programmiersprachen) sind wir von einer anwendungsfähigen Implementation immer noch einen Schritt entfernt, da wir nicht die Möglichkeit haben, den Überlauf zu unterdrücken; er ist als Fehlerbedingung definiert, die zu nicht vorhersagbaren Ergebnissen führen kann. Nehmen wir an, daß unser Computer ein 32-Bit-Wort hat, und wählen wir m = 100000000, b = 31415821 und als Anfangswert a = 1234567. Alle diese Werte sind deutlich kleiner als die größte darstellbare ganze Zahl, doch bereits die erste Operation a * b + 1 führt zu einem Überlauf. Der Teil des Produkts, der den Überlauf verursacht, ist für unsere Berechnung jedoch nicht von Bedeutung; uns interessieren nur die letzten acht Ziffern. Der Trick besteht darin, den Überlauf zu vermeiden, indem die Multiplikation in Teile zerlegt wird. Um p mit q zu multiplizieren, schreiben wir p = 104p1 + p0 und q = 104q1 + q0, so daß das Produkt

pq = (104p1 + p0)(104q1 + q0)
= 108p1q1 + 104(p1q0 + p0q1) + p0q0.
lautet.

Da wir nun für das Ergebnis nur acht Ziffern benötigen, können wir den ersten Term und alle außer den letzten vier Ziffern des zweiten Terms (vor der Multiplikation) ignorieren. Dies führt zu dem folgenden Programm:

    program random(input,output);
    const m=100000000; m1=10000; b=31415821;
    var i,a,N: integer;
    function mult(p,q: integer): integer;
      var p1,p0,q1,q0: integer;
      begin
      p1:=p div m1; p0:=p mod m1;
      q1:=q div m1; q0:=q mod m1;
      mult:=(((p0*q1+p1*q0) mod m1)*m1+p0*q0) mod m;
      end;
    function random: integer;
      begin
      a:=(mult(a,b)+1) mod m;
      random:=a;
      end;
    begin
    read(N,a);
    for i:=1 to N do writeln(random)
    end.

Die Funktion mult in diesem Programm berechnet p * q mod m, wobei kein Überlauf auftritt, solange m kleiner ist als die Hälfte der größten darstellbaren ganzen Zahl. Das Verfahren kann offensichtlich mit m = m1 * m1 auch für andere Werte von m angewandt werden.

Wenn dieses Programm mit den Eingabedaten N = 10 und a = 1234567 abgearbeitet wird, gibt es die folgenden zehn Zahlen aus: 35884508, 80001069, 63512650, 43635651, 1034472, 87181513, 6917174, 209855, 67115956, 59939877. Diesen Zahlen wohnt offensichtlich ein nicht zufälliges Element inne: Zum Beispiel durchlaufen die Endziffern zyklisch die Ziffernfolge 0-9. Anhand der Formel läßt sich leicht beweisen, daß dies der Fall sein muß. Allgemein gesagt, sind die rechts stehenden Ziffern nicht besonders zufällig, eine Tatsache, die die Ursache für einen verbreiteten und ernsthaften Fehler bei der Benutzung von derartigen Zufallszahlengeneratoren ist. Das folgende ist ein schlechtes Programm für die Erzeugung von Zufallszahlen im Bereich [0,r-1]:

    function randombad(r: integer): integer;
      begin
      a:=(mult(b,a)+1) mod m;
      randombad:=a mod r;
      end;

Die nicht zufälligen Ziffern auf der rechten Seite sind die einzigen verwendeten Ziffern, so daß die sich ergebende Folge wenige der gewünschten Eigenschaften besitzt. Dieses Problem läßt sich leicht beheben, indem man stattdessen die links stehenden Ziffern benutzt. Wir können eine Zahl zwischen 0 und r-1 berechnen, indem wir a * r div m bilden, doch auch diesmal muß der Überlauf umgangen werden:

    function randomint(r: integer): integer;
      begin
      a:=(mult(a,b)+1) mod m;
      randomint:=((a div m1)*r) div m1
      end;

Ein anderes gebräuchliches Verfahren besteht darin, zufällige reelle Zahlen zwischen 0 und 1 zu erzeugen, indem man die obigen Zahlen als Nachkommastellen eines Dezimalbruchs auffaßt. Dies läßt sich implementieren, indem einfach der reelle Wert a/m statt der ganzen Zahl a zurückgegeben wird. Ein Anwender kann dann eine ganze Zahl im Bereich [0,r] erhalten, indem er diesen Wert einfach mit r multipliziert und das Ergebnis auf die nächstkleinere ganze Zahl abrundet. Es ist auch möglich, daß eine reelle Zahl zwischen 0 und 1 genau das ist, was benötigt wird.


[ Inhaltsverzeichnis ] [ vorhergehende Seite ] [ nächste Seite ] [ Stichwort ]