Inverser Kongruenzgenerator

Ein inverser Kongruenzgenerator ist ein arithmetischer Zufallszahlengenerator, der durch den Satz von Marsaglia bekannte Nachteile linearer Kongruenzgeneratoren vermeidet. Insbesondere lässt er keine Hyperebenen entstehen. Verwendet man Zufallszahlen inverser Kongruenzgeneratoren für die Box-Muller-Methode, so wird ein Spiralverhalten vermieden. Im Gegenzug verlangt er einen höheren Rechenaufwand.

Er besteht aus 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 nach folgendem Bildungsgesetz:

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

Zur Erklärung der Symbolik siehe 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.

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.

Im Gegensatz zu linearen Kongruenzgeneratoren, deren Werte ja auf wenigen Hyperebenen liegen, kann man hier zeigen, dass 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.

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}})}$

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]; // Initialisiert 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
    }
}

Manchmal liest man auch 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 keine Verallgemeinerung dar; man erhält durch Ausmultiplizieren sofort die obige Gestalt.

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

• Kongruenzgenerator
• Erweiterter euklidischer Algorithmus

• 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).