Mersenne Twister

Der Mersenne-Twister ist ein Pseudozufallszahlengenerator, der 1997 von Makoto Matsumoto und Takuji Nishimura entwickelt wurde. Er ermöglicht die schnelle Erzeugung hochwertiger Sequenzen von Pseudozufallszahlen und wurde extra darauf zugeschnitten, die Probleme älterer Algorithmen zu überwinden.

• Er hat die extrem lange Periode von 2¹⁹⁹³⁷-1 (≈4,3·10⁶⁰⁰¹). Diese Periodenlänge erklärt auch den Namen des Algorithmus: Sie ist eine Mersenne-Primzahl, und einige Eigenschaften des Algorithmus resultieren daraus.
• Er liefert hochgradig gleichverteilte Sequenzen (bewiesene Gleichverteilung bis zur Dimension 623, siehe unten). Daraus folgt eine extrem kleine Korrelation zwischen aufeinanderfolgenden Wertefolgen der Ausgabesequenz.
• Er ist schneller als jeder andere bekannte (hinreichend gute) Algorithmus.
• Alle Bits der Ausgabesequenz sind für sich gleichverteilt.

Andererseits hat er den Nachteil, auf einer großen Datenmenge von 2,5 kB (624 Wörter mit je 32 Bits) zu arbeiten. Mit den heutigen Rechnerarchitekturen mit schnellem, aber relativ kleinem Cache und langsamerem Arbeitsspeicher kann sich daraus ein Geschwindigkeitsnachteil ergeben.

Der Mersenne-Twister ist ein verdrillter Generalized Feedback Shift Register Algorithmus, kurz (GFSR), von rationaler Normalenform (TGFSR(R)), mit Statusbit Reflection und Tempering, mit dem die Multiplikation einer Matrix T gemeint ist. Er wird charakterisiert durch die folgenden Einheiten:

• w: Wortgröße (in Anzahl an Bits)
• n: Grad der Rekursionen
• m: mittleres Wort, oder die Anzahl an parallelen Sequenzen, 1 ≤ m ≤ n
• r: Trennpunkt eines Wortes, oder die Anzahl an Bits der niederen Bitmaske, 0 ≤ r ≤ w - 1
• a: Koeffizienten der verdrillten Matrix in rationaler Normalenform
• b, c: TGFSR(R) Tempering Bitmasken
• s, t: TGFSR(R) Tempering Bit Shifts
• u, l: zusätzliche Mersenne Twister Tempering Bit Shifts

mit der Einschränkung das 2^nw − r − 1 eine Mersenne Primzahl ist. Diese Wahl vereinfacht den Primitivitätstest und k-Verteilungstest, welche für die Parametersuche benötigt werden.

Für ein Wort x mit der Bitlänge w, wird es als rekursives Verhältnis, wie folgt ausgedrückt:

mit | als das bitweise ODER und ⊕ als das bitweise Exklusiv ODER (XOR), während x^u, x^l eine Anwendung der höheren und niederen Bitmasken auf das x darstellen. Die Twist Transformation A ist in der rationalen Normalenform definiert.

mit I_n − 1 als die (n − 1) × (n − 1) Identitätsmatrix (und im Gegensatz zu der normalen Matrixmultiplikation, ersetzt das bitweise XOR die Addition).

Ähnlich wie der TGFSR(R), ist der Mersenne Twister in einer Kaskade geschalten mit einer Tempering Transformation um die reduzierte Dimensionsgröße der Gleichverteilung (Weylsche “Gleichverteilung von Zahlen mod 1”) zu kompensieren. Diese ist äquivalent zu der Transformation A = R → A = T⁻¹RT, T umkehrbar. Das Tempering ist beim Mersenne Twister als

y := x ⊕ (x >> u)

y := :y ⊕ ((y << s) & b)

y := :y ⊕ ((y << t) & c)

z := y ⊕ (y >> l)

definiert.

Zusammenfassend kann man den Mersenne Twister in 6 einfache Einzelschritte unterteilen. Zunächst wird eine Maske für die höheren und niederen Bits generiert. Anschließend wird ein Array x mit einem Seed initialisiert. Danach werden die höheren Bits y von x[i] mit den niederen von x[i+1] verknüpft. Nun wird y mit der Matrix A multipliziert. A wird dabei vorsichtig gewählt, damit es leicht mit Rechtsshifts und EXOR berechnet werden kann. Nun kommt das Tempering ins Spiel, bei dem mit T multipliziert wird um eine bessere Weylsche Gleichverteilung und Genauigkeit zu erreichen. Im letzten Schritt wird i um 1 inkrementiert und der Vorgang wiederholt.

Der Algorithmus ist im Detail auf der offiziellen Seite beschrieben. Dort findet sich auch der originale, in C implementierte Quellcode.

Die Kombination der gewonnenen Erkenntnisse

Um einen zuverlässigen PRNG zu schreiben sind die bisher beschriebenen Punkte von großer Bedeutung. Aus der Kombination von Wissen über die Wahrscheinlichkeitsverteilungen und moderne Algorithmen zur Generierung von Zufallszahlen, lässt sich das nun bewerkstelligen. Die gezeigten Methoden zur Berechnung statistisch verteilter Zufallszahlen werden bereits in vielen Bibliotheken vorimplementiert. So verfügt C++ Boost über eine eigene ”Boost Random Number Library“, die diverse Zufallszahlengeneratoren und Verteilungen zur Verfügung stellt.

Unter Nutzung dieser Bibliothek lässt sich beispielsweise ein Mersenne Twister Generator schreiben, der mithilfe der Gaußschen Normalverteilung Zufallszahlen erzeugt. Die folgende Funktion mit dem Namen SampleNormal gibt Proben von Zufallszahlen aus einer Normalenverteilung zurück. Der Mittel- oder Erwartungswert und die Standardabweichung der Verteilung werden als Parameter an die Funktion übergeben.

#include <boost/random.hpp>
#include <ctime>
 
using namespace boost;
 
double SampleNormal(double mean, double sigma)
{
    // Create a Mersenne twister random number generator
    // that is seeded once with #seconds since 1970
    static mt19937 rng(static_cast<unsigned>(std::time(0)));
 
    // select Gaussian probability distribution
    normal_distribution<double> norm_dist(mean, sigma);
 
    // bind random number generator to distribution, forming a function
    variate_generator<mt19937&, normal_distribution<double>> normal_sampler(rng, norm_dist);
 
    // sample from the distribution
    return normal_sampler();
}

Aus den Proben lässt sich ein sogenanntes Histogramm generieren, mit dem die Häufigkeitsverteilung der Zufallszahlen anschaulich graphisch dargestellt werden kann. Das folgende Histogramm zeigt sehr unregelmäßige Proben die eine Normalenverteilung aufweisen.

Bei einem echten Zufallszahlengenerator wären die hellblau gefärbten Proben über das gesamte Spektrum absolut gleichverteilt und gleich hoch. In unserem Fall wird die Verteilung der Werte, aufgrund der gewählten Verteilungsfunktion, sinngemäß ebenfalls der Gaußschen Glockenkurve nahe kommen.
Jeder PRNG erfüllt in keinem Fall das Kriterium eines einheitlichen Histogramms. Er kann diesem Kriterium zwar nahe kommen, allerdings wird er es nie zu 100 % erfüllen können. Nichts desto trotz liefert der Mersenne Twister sehr gleichverteilte Werte über das gesamte Spektrum.

Weitere Generatoren

Neben dem eingehend beschriebenen Mersenne Twister Zufallszahlengenerator existieren auch noch weitere bekannte Generatoren, die in der Praxis oft zum Einsatz kommen. Dazu zählen der Lagged Fibonacci Generator (LFG), der auf Fibonacci Sequenzen basiert, der Linear Congruential Generator, der zu den ältesten und bekanntesten Zufallszahlengeneratoren gehört bis hin zum alten Blum Blum Shub (BBS) Algorithmus, der sich heute aufgrund seiner Trägheit kaum noch für Simulationen, sondern hauptsächlich für Verschlüsselungen eignet.

In einem letzten Codebeispiel generieren wir einen Lagged Fibonacci Generator und wenden eine Normalenverteilung auf die Zufallszahlen an.

#include <boost/random/normal_distribution.hpp>
#include <boost/random/lagged_fibonacci.hpp>
#include <iostream>
#include <vector>
 
using namespace boost;
 
int main()
{
    const double mean = 10.0;
    const double sigma = 1.0;
    normal_distribution<double> norm_dist(mean, sigma);
    lagged_fibonacci44497 engine;
 
    // Make a histogram
    const int size = static_cast<int>(mean) * 2;
    std::vector<int> histo(size, 0);
    for(int i = 0; i != 1000000; ++i) {
        const double value = norm_dist.operator()<lagged_fibonacci44497>((engine));
        int index = static_cast<int>(value);
        index = index < 0 ? 0 : index;
        index = index > size - 1 ? size - 1 : index;
        ++histo[index];
    }
 
    // Output histogram
    for(int i = 0; i != size; ++i) {
        std::cout << histo[i] << std::endl;
    }
 
    return 0;
}

Schluss

Damit sind wir am Ende dieses Tutorials angekommen. Im Anhang finden Sie einen umfassenden Quellcode zu dem Thema. Es handelt sich dabei um eine kleine C++ Bibliothek, die unter anderem die Gammaverteilung und den Mersenne Twister implementiert. Sie können die Funktionen der Bibliothek frei nutzen und sich auf diese Weise beliebige Zufallszahlen für ihr Projekt generieren lassen.

Weitergehende Informationen

• Analysis of the Linux Random Number Generator: Eine ausführliche Analyse des Linux Zufallszahlengenerators von Zvi Gutterman (Safend and The Hebrew University of Jerusalem), Benny Pinkas (University of Haifa) und Tzachy Reinman (The Hebrew University of Jerusalem)