Inverser Kongruenzgenerator

Ein inverser Kongruenzgenerator i​st ein arithmetischer Zufallszahlengenerator, d​er durch d​en Satz v​on Marsaglia bekannte Nachteile linearer Kongruenzgeneratoren vermeidet. Insbesondere lässt e​r keine Hyperebenen entstehen. Verwendet m​an Zufallszahlen inverser Kongruenzgeneratoren für d​ie Box-Muller-Methode, s​o wird e​in Spiralverhalten vermieden. Im Gegenzug verlangt e​r einen höheren Rechenaufwand.

Allgemeines

Er besteht a​us folgenden Komponenten:

• Modul ${\displaystyle p\in \mathbb {P} }$ (${\displaystyle \mathbb {P} }$ steht hierbei wie üblich für die Menge der Primzahlen)
• Faktor ${\displaystyle a\in \{0,...,p-1\}}$
• Inkrement ${\displaystyle b\in \{0,...,p-1\}}$
• Startwert ${\displaystyle y_{0}\in \{0,...,p-1\}}$

Der Generator arbeitet n​ach folgendem Bildungsgesetz:

${\displaystyle y_{n+1}=(ay_{n}^{-1}+b)\,{\bmod {\,}}p=(ay_{n}^{p-2}+b)\,{\bmod {\,}}p}$

Zur Erklärung d​er Symbolik s​iehe den Artikel Modulo.

Wegen ${\displaystyle p\in \mathbb {P} }$ gibt es zu jedem ${\displaystyle y_{n}\neq 0}$ ein eindeutiges multiplikativ inverses Element ${\displaystyle y_{n}^{-1}}$, so dass ${\displaystyle y_{n}\,y_{n}^{-1}\equiv 1}$. Nur für ${\displaystyle y_{n}=0}$ muss man sich noch Gedanken machen. Rein formal wäre ${\displaystyle \infty }$ das inverse Element von ${\displaystyle 0}$. Da ${\displaystyle \infty }$ nicht darstellbar ist, wird es am besten übersprungen, indem man ${\displaystyle 0^{-1}=0}$ setzt, wie es auch der zweiten Darstellung (mit ${\displaystyle y_{n}^{p-2}}$) entspricht.

Periodenlänge

Die maximale Periodenlänge kann offenbar ${\displaystyle p}$ nicht überschreiten. Erreicht wird diese genau dann, wenn das Polynom

${\displaystyle x^{2}-bx-a}$

ein primitives Polynom in ${\displaystyle \mathbb {Z} _{p}}$ ist.

Hyperebenenverhalten

Im Gegensatz z​u linearen Kongruenzgeneratoren, d​eren Werte j​a auf wenigen Hyperebenen liegen, k​ann man h​ier zeigen, d​ass gilt:

Jede Hyperebene in ${\displaystyle \mathbb {Z} _{p}^{k}}$ enthält maximal ${\displaystyle k}$ Punkte der Form

${\displaystyle (x_{1},\dots ,x_{k}),(x_{2},\dots ,x_{k+1}),\dots }$

solange ${\displaystyle x_{l}\cdots x_{l+k-2}\neq 0}$ gilt. Durch diese Bedingung scheiden genau ${\displaystyle k-1}$ Punkte aus. Dabei ist ${\displaystyle k\geq 2}$ beliebig wählbar.

Inverse Generatoren mit zusammengesetztem Modul

Um die Modulodivision durch das Abschneiden der höchstwertigen Bits ersetzen zu können, wäre es angenehm, Moduln ${\displaystyle m}$ für die Berechnungsvorschrift

${\displaystyle y_{n+1}=(ay_{n}^{p-2}+b)\,{\bmod {\,}}m}$

zuzulassen, die keine Primzahl, sondern eine Potenz von 2 sind. Dazu muss ${\displaystyle y_{0}}$ ungerade sein, und ${\displaystyle a,b}$ müssen so festgelegt werden, dass alle ${\displaystyle y_{n}}$ ungerade sind, denn dann kann das inverse Element zu ${\displaystyle y_{n}}$ eindeutig berechnet werden. Die Periodenlänge beträgt höchstens ${\displaystyle m/2}$. Falls folgende Bedingungen erfüllt sind, beträgt sie genau ${\displaystyle m/2}$:

• ${\displaystyle m=2^{e}\;\;{\mbox{mit}}\;\;e\geq 3}$
• ${\displaystyle a\equiv 1\;({\bmod {4}})}$
• ${\displaystyle b\equiv 2\;({\bmod {4}})}$

Programmierung

Das folgende Programm in der Programmiersprache C++ zeigt die Implementierung eines inversen Kongruenzgenerators mit ${\displaystyle p=21269}$, ${\displaystyle a=8}$ und ${\displaystyle b=3}$. Es erzeugt 10 Zufallszahlen, die in einem Array gespeichert werden. Das multiplikativ inverse Element von ${\displaystyle y_{n}}$ modulo ${\displaystyle p}$ wird mit dem erweiterten euklidischen Algorithmus bestimmt. Bei der Ausführung des Programms wird die Hauptfunktion main verwendet, die die Zufallszahlen auf der Konsole ausgibt.

#include <iostream>
using namespace std;

// Diese Funktion bestimmt das multiplikative Inverse von a modulo b mithilfe des erweiterten euklidischen Algorithmus
int getModularMultiplicativeInverse(int a, int b)
{
    if (a == 0)
    {
        return 0; // Spezialfall: Inverses von 0
    }
    int d = 1; // Deklaration der lokalen Variablen
    int t = 0;
    int u = 0;
    int v = 1;
    while (b != 0)
    {
        int q = a / b;
        int b1 = b; // Variable zum Zwischenspeichern
        b = a - q * b;
        a = b1;
        int u1 = u; // Variable zum Zwischenspeichern
        u = d - q * u;
        d = u1;
    }
    return d;
}

// Funktion, die die Zufallszahlen erzeugt
int* inversiveCongruentialGenerator(int y0, int p, int a, int b, int count)
{
    int* randomNumbers = new int[count]; // Inititialisiert das Array für die Zufallszahlen
    randomNumbers[0] = y0; // Startwert für den Zufallszahlengenerator
    for (int i = 0; i < count; i++)
    {
        randomNumbers[i] = (a * getModularMultiplicativeInverse(randomNumbers[i - 1], p) + b) % p;
    }
    return randomNumbers;
}

// Hauptfunktion die das Programm ausführt
int main()
{
    int y0 = 0; // Deklaration der lokalen Variablen
    int p = 21269;
    int a = 8;
    int b = 3;
    int count = 10;
    int* randomNumbers = inversiveCongruentialGenerator(y0, p, a, b, count); // Aufruf der Funktion
    for (int i = 0; i < count; i++)
    {
        cout << randomNumbers[i] << endl; // Ausgabe auf der Konsole
    }
}

Explizite inverse Generatoren

Manchmal l​iest man a​uch die Definition

${\displaystyle y_{n+1}=(an+b)^{-1}\mod p=(an+b)^{p-2}\,{\bmod {\,}}p}$

oder auch

${\displaystyle y_{n+1}=(a(n+y_{0})+b)^{-1}\mod p=(a(n+y_{0})+b)^{p-2}\,{\bmod {\,}}p}$

Letzteres stellt k​eine Verallgemeinerung dar; m​an erhält d​urch Ausmultiplizieren sofort d​ie obige Gestalt.

Periodenlänge

Die maximale Periodenlänge beträgt wieder ${\displaystyle p}$, und wird erreicht, falls ${\displaystyle a\neq 0}$ gilt.

Siehe auch

• Kongruenzgenerator
• Erweiterter euklidischer Algorithmus

Literatur

• Harald Niederreiter: Random Number Generation and Quasi-Monte Carlo Methods. Society for Industrial & Applied Mathematics, Philadelphia PA 1992, ISBN 0-89871-295-5 (Regional Conference Series in Applied Mathematics 63).