Generatoren
Die Rekursions- Relation
Alle üblichen Methoden zur Erzeugung von Pseudozufallszahlen basieren auf der folgenden Rekursions- Relation:
\begin{displaymath}
x_{n+1} = mod(a_{0} x_{n} + a_{1} x_{n-1} + ..... + a_{j} x_{n-j} + b, M)
\end{displaymath} (1)

Hierbei sind die $x_{n}, x_{n-1}, ..., x_{n-j}$ die in den $j$ vorherigen Aufrufen berechneten Zufallszahlen. Alle auf der rechten Seite der Relation auftretenden Parameter sowie die Startwerte $x_{0}, x_{1}, ... , x_{j}$ dürfen keine gemeinsamen Teiler außer der 1 haben. Diese Forderung folgt aus der Zahlentheorie. Beschränkt man sich auf Primzahlen, so ist diese Bedingung immer erfüllt. $M$ sollte möglichst groß gewählt werden, am besten gleich der größten darstellbaren Zahl auf dem jeweiligen Rechner. Durch diese Wahl wird nicht nur die Dimension des Ereignisraumes möglichst groß, sondern, wie wir im einzelnen noch sehen werden, auch die Qualität des Generators wird verbessert. Auf alle diese genannten Bedingungen werden wir später noch detaillierter eingehen. Die aus der Rekursions- Relation folgende Folge ganzer Zahlen im Intervall $\{ 0, M \}$ kann auf eine Folge im Intervall $\{ 0, 1 \}$ abgebildet werden gemäß $x=x_{n}/M$.

Wie oben schon gesagt, sollte $M$ möglichst groß gewählt werden. Die größte darstellbare Zahl des Rechners ist $2^{n-1}-1$, wobei $n$ also die Anzahl der bits in der arithmetischen Einheit ist, einschließlich des Vorzeichenbits. Wählt man $R=2^{n-1}$, so nutzt man den vollen Zahlenbereich der positiven ganzen Zahlen aus. Das mathematische Resultat $C$ der Summe

\begin{displaymath}
C = \sum_{i=0}^{j} a_{i} x_{n-i} +b
\end{displaymath}

ist im allgemeinen größer als $M$, also nicht mehr als Zahl auf dem Rechner darstellbar. Hier kann man jetzt die Überlauf- Arithmetik ausnutzen, die wir im vorigen Kapitel ausführlich diskutiert haben. In der 2-er Komplement Darstellung hatten wir gezeigt, daß

\begin{displaymath}
mod(C,2^{n-1}) = C_{R} \; \; \; f''ur \; \; C_{R} \geq 0
\end{displaymath}


\begin{displaymath}
mod(C,2^{n-1}) = C_{R} + (2^{n-1} -1) +1 \; \; \; f''ur \; \; C_{R} < 0 ,
\end{displaymath}

wobei $C_{R}$ das vom Rechner gelieferte Resultat ist. Auf der linken Seite dieser Gleichung steht aber jetzt genau unsere Rekursions- Relation (1). Wählt man also $M=2^{n-1}$, so ergibt sich folgende Programmier- Vorschrift: 1. wenn $C_{R} \geq 0$, dann ist $x_{n+1} = C_{R}$; 2. wenn $C_{R} < 0$, dann bilde zunächst $C_{R} + (2^{n-1}-1)$ und addiere danach noch einmal eine 1 hinzu. Das Ergebnis ist $mod(C,2^{n-1})$. Der Ausschnitt aus einem Simulationsprogramm zur Erzeugung von Pseudozufallszahlen auf einem 32-bit Rechner würde wie in Tabelle 1 aussehen.

Tabelle 1: Programm einer Rekusionsrelation
......
int a = 69060;
int ran = 1;
......
ran *= a;
if(ran $<$ 0) ran += 2147483647 + 1;

Auf normalen Taschenrechnern gibt es allerdings die Festkomma- Darstellung ganzer Zahlen nur für die Adressierung (z.B. in Schleifen- Anweisungen). Taschenrechner verwenden ausschließlich Gleitkomma- Darstellung. Hier können wir die Überlauf- Arithmetik nicht anwenden. Stattdessen arbeiten Taschenrechner aber mit doppelter Genauigkeit bei der Darstellung der Mantisse, nämlich mit 64 bit oder 12 signifikanten Dezimalstellen. Solange das mathematische Resultat $C$ noch ohne größere Rundungsfehler von unserem Taschenrechner berechnet wird, können wir einfach die Modulo- Funktion auf das Resultat $C$ anwenden. Der obige Programm- Auschnitt wird dann ersetzt durch die Zeile
$C = C - 32768 * INT(C/32768) $
oder, falls die Modulo- Funktion als Systemroutine auf dem Rechner vorhanden ist, einfach durch
$C = C \% 32768$.
Wir werden später in diesem Abschnitt eine Reihe von weiteren Zufallszahlengeneratoren diskutieren. Zunächst jedoch müssen wir zwei statistische Besonderheiten kennenlernen, nämlich Histogramme und den $\chi^{2}$- Test.

Histogramme
Wir hatten gelernt, daß die Dimension des Ereignisraumes bei der Simulation von Pseudozufallszahlen durch den Parameter $M$ in der Rekursions- Relation (1) gegeben war. $M$ sollte möglichst groß gewählt werden, sodaß die Wahrscheinlichkeit für die Simulation einer bestimmten Zahl $x_{0}$ außerordentlich klein ist, nämlich $P(x_{0})=1/M$. Wir können aber das Zahlenintervall $[0,M]$ bzw. $[0,1]$ in $m$ kleinere Intervalle aufteilen, sodaß alle x mit

\begin{displaymath}
x_{\mu -1} < x \leq x_{\mu}, \; \; \mu = 1,2,3,...,m
\end{displaymath}

zum $\mu$-ten Intervall gehören. Die Wahrscheinlichkeit $P_{\mu}$ dafür, daß ein bestimmtes $x$ im $\mu$-ten Intervall liegt, ist dann gegeben durch
\begin{displaymath}
P_{\mu} = \sum_{x_{\mu -1} < x_{j} \leq x_{\mu}} P(x_{j}) \approx
\int_{x_{\mu -1}}^{x_{\mu}} dx p(x)
\end{displaymath} (2)

mit

\begin{displaymath}
\sum_{\mu =1}^{m} P_{\mu} =1.
\end{displaymath}

Führen wir den Versuch V $n$-mal aus (bzw., rufen wir den Zufallszahlengenerator $n$-mal auf), so erhalten wir die $m$ Häufigkeiten

\begin{displaymath}
k_{1}, k_{2}, ..... , k_{m}
\end{displaymath}

mit

\begin{displaymath}
\sum_{\mu =1}^{m} k_{\mu} = n.
\end{displaymath}

Die mathematischen Erwartungswerte für die $k_{\mu}$ sind durch $n P_{\mu}$ gegeben. Die gemessenen Häufigkeiten $k_{\mu}$ und die mathematischen Erwartungswerte vergleichen wir in einem Histogramm .

Häufig werden in grafischen Darstellungen die theoretisch erwarteten Werte durch eine Kurve verbunden. Diese Kurven sind aber vom statistischen Standpunkt aus sehr problematisch, insbesondere dann, wenn die Intervall- Aufteilung grob ist. Durch die Summenbildung bzw. Integralbildung verlieren wir nämlich Informationen über die mathematische Wahrscheinlichkeitsverteilung, sodaß die theoretisch erwartete Verteilungsfunktion und die Kurve, die man durch Verbindung der Integralwerte erhält, verschieden sein können.

Eine Aufteilung in gleich große Intervalle ist eigentlich nur bei der Gleichverteilung vernünftig. Im allgemeinen sollte die Intervall- Aufteilung so durchgeführt werden, daß in jedem Intervall ungefähr gleich viele Versuche enthalten sind, oder, anders ausgedrückt, daß die Wahrscheinlichkeiten $P_{\mu}$ gleich sind für alle $\mu$. Dieses erreicht man mit der Intervalleinteilung

\begin{displaymath}
\int_{-\infty}^{x_{\mu}} dx p(x) = \frac{\mu}{m}, \; \; \mu = 1,2,...,m.
\end{displaymath} (3)

Mit Hilfe von (2) erhält man

\begin{displaymath}
P_{\mu} = \frac{1}{m}, \; \; \mu = 1,2,...,m.
\end{displaymath}

Für die gemessenen Häufigkeiten erwarten wir

\begin{displaymath}
k_{\mu} \approx \frac{n}{m}, \; \; \mu = 1,2,...,m.
\end{displaymath}

Der $\chi^{2}$- Test
Die gemessene Wahrscheinlichkeitsverteilung stimmt umso besser mit der mathematisch erwarteten Verteilung überein, je kleiner die Abweichungen zwischen mathematischer und gemessener Verteilung im Histogramm sind. Diese triviale qualitative Aussage müssen wir quantitativ formulieren. Hierzu dient der $\chi^{2}$- Test. Wir schreiben die Abweichungen

\begin{displaymath}
\delta k_{\mu} = k_{\mu} - P_{\mu} n, \; \; \mu = 1,2,...,m ,
\end{displaymath}

berechnen weiterhin die Varianzen $\sigma_{\mu}$ für die gemessenen Werte $k_{\mu}$,

\begin{displaymath}
\sigma_{\mu} = \sqrt{n P_{\mu} (1 - P_{\mu})} \approx \sqrt{n P_{\mu}},
\end{displaymath}

und bilden die sogenannte mittlere quadratische Abweichung
\begin{displaymath}
\Theta = \frac{1}{m-1} \sum_{\mu =1}^{m}
\frac{(\delta k_{\...
...{m}
\frac{(k_{\mu} - P_{\mu} n)^{2}}{n P_{\mu} (1 - P_{\mu})}
\end{displaymath} (4)

Je größer also die theoretisch erwartete Streuung der Werte $k_{\mu}$ um den Mittelwert $P_{\mu} n$ ist, desto kleiner ist der Beitrag in der Summe der quadratischen Abweichungen. $\Theta$ ist eine statistische Größe, da sie aus statistischen Größen $k_{\mu}$ berechnet wird. $\Theta$ besitzt also eine Dichtefunktion. Wir werden später im Detail zeigen, daß der Mittelwert und die Varianz durch
$\displaystyle < \Theta >$ $\textstyle =$ $\displaystyle 1$ (5)
$\displaystyle < (\Theta - < \Theta > )^{2} >$ $\textstyle =$ $\displaystyle \frac{2}{m-1}$ (6)

gegeben sind. Es erweist sich daher als günstig, eine andere Größe zu definieren, nämlich
\begin{displaymath}
\Theta' = (\Theta -1) \sqrt{\frac{m-1}{2}}.
\end{displaymath} (7)

Diese Größe hat den Mittelwert und Varianz
$\displaystyle < \Theta' >$ $\textstyle =$ $\displaystyle 0$ (8)
$\displaystyle < ( \Theta' - < \Theta' > )^{2} >$ $\textstyle =$ $\displaystyle 1.$ (9)

Daraus folgt, daß ein $\Theta'$- Wert mit
\begin{displaymath}
\alpha_{1} < \Theta' < \alpha_{2}
\end{displaymath} (10)

und

\begin{displaymath}
- \alpha_{1} \approx \alpha_{2} \approx 2
\end{displaymath}

anzeigt, daß die hypothetisch postulierten Wahrscheinlichkeiten $p_{\mu}$ und die gemessenen relativen Häufigkeiten $k_{\mu}/n$ mit $\approx 95 \%$ Vertrauensbereich miteinander verträglich sind.

Diesen $\chi^{2}$- Test und weitere Tests zur Prüfung von Hypothesen werden in einem späteren Kapitel noch ausführlicher diskutiert. Wir werden dann zeigen daß $\alpha_{1} \approx -2$ und $\alpha_{2} \approx 4$ gewählt werden muss. Im folgenden werden wir eine Hypothese akzeptieren, wenn beim $\chi^{2}$- Test ein $\Theta'$- Wert mit

\begin{displaymath}
-2 < \Theta' < 4
\end{displaymath} (11)

erhalten wird.

Der Korrelationstest
Mit dem $\chi^{2}$- Test können wir prüfen, ob die zugrundeliegende Zahlenfolge des Zufallszahlen- Generators mit der Hypothese einer gleichverteilten Veränderlichen verträglich ist. Dieser Test sagt aber nichts über Korrelationen aus. So ist z.B. die Zahlenfolge (mit $a_{0} = b = 1, M = 100$ und $a_{i} = 0$ für $i > 1$, siehe (1)

\begin{displaymath}
x_{i+1} = mod(x_{i} + 1,100), \; \; \; x_{0} = 0,
\end{displaymath} (12)

perfekt gleichverteilt im Intervall $0 \leq x < 100$. Andererseits ist die daraus resultierende Zahlenfolge extrem korreliert. Man kann sich leicht davon überzeugen, daß diese Rekursion keinen guten Zufallszahlen- Generator liefert. Auch bei der allgemeinen Rekursions- Relation (1) wird jede folgende Zufallszahl aus den vorhergehenden Zahlen berechnet. Aufeinanderfolgende Zahlen können daher, zumindest im streng statistischen Sinne, nicht unabhängig voneinander sein. Was bedeutet also bei Zufallszahlen der Begriff ,,unabhängig'' ?

Im folgenden Reihen- Korrelationstest prüfen wir, ob die Zufallszahl $x_{i+1}$ mit der vorher generierten Zahl $x_{i}$ korreliert ist. Gegeben seien $n$ Zufallszahlen

\begin{displaymath}
x_{i} = \{ x_{1}, x_{2}, x_{3}, ....., x_{n-1}, x_{n} \} .
\end{displaymath}

Allgemein ist der Korrelationskoeffizient bei 2 diskreten Veränderlichen $x$ und $y$ durch

\begin{displaymath}
\rho = \frac{n \sum_{i=1}^{n} x_{i} y_{i} - \sum_{i=1}^{n} x...
... \sum_{i=1}^{n} y_{i}^{2} -
( \sum_{i=1}^{n} y_{i} )^{2} ] }}
\end{displaymath}

gegeben. Setzen wir hierin $x_{i} = x_{i}$ und $y_{i}=x_{i+1}$, so folgt mit der Festsetzung $x_{n+1} = x_{1}$
\begin{displaymath}
\rho = \frac{ n \sum_{i=1}^{n} x_{i} x_{i+1} - ( \sum_{i=1}^...
... n \sum_{i=1}^{n} x_{i}^{2} - ( \sum_{i=1}^{n} x_{i} )^{2} } .
\end{displaymath} (13)

Dieses ist die Definition des Reihen- Korrelationskoeffizienten. Da $\rho$ aus zufälligen Veränderlichen berechnet wurde, muß $\rho$ ebenfalls durch eine Verteilungsfunktion beschreibbar sein. Wir wissen, daß $\rho$ zwischen $-1$ und $+1$ liegt. Für unabhängige Veränderliche $x$ muß $\rho$ sogar in der Nähe von Null liegen. Wir erwarten also

\begin{displaymath}
< \rho > \approx 0 .
\end{displaymath}

Die exakte Verteilungsfunktion ist nicht bekannt. Die empirischen Erfahrungswerte für einen ''guten'' Zufallszahlen- Generator sind:
\begin{displaymath}
\mu_{n} - 2 \sigma_{n} < \rho < \mu_{n} + 2 \sigma_{n}
\end{displaymath} (14)

mit
$\displaystyle \mu_{n}$ $\textstyle =$ $\displaystyle \frac{1}{1-n}$  
$\displaystyle \sigma_{n}$ $\textstyle =$ $\displaystyle \frac{1}{n-1} \sqrt{ \frac{n(n-3)}{n+1} } .$  

Schauen wir uns darauf hin noch einmal die Rekursions- Relation (12) als Beispiel eines extrem schlechten Zufallszahlen- Generators an. Wir erhalten (mit n=99):

\begin{displaymath}
\rho = \frac{n \sum_{i=0}^{n} i ( i+1 ) - (\sum_{i=0}^{n} i ...
... \sum_{i=0}{n} i^{2} -
( \sum_{i=0}^{n} i )^{2} } \approx 1.
\end{displaymath}

Für den Generator (12) erhält man eine maximale Korrelation.

Wir fassen zusammen: Die Unabhängigkeit aufeinanderfolgender Zahlen in einer Folge von Zufallszahlen ist dadurch definiert, daß der Reihen- Korrelationskoeffizient (13) der Bedingung (14) genügt. Dieses heißt aber nicht, daß die Zahlenfolge unabhängige Veränderliche im exakt statistischen Sinne liefert. Das Verschwinden des Korrelationskoeffizienten ist notwendig, aber nicht hinreichend. Daher nennt man diese Zufallszahlen auch Pseudo- Zufallszahlen.

In einer viel beachteten Arbeit hat Marsaglia gezeigt, daß, obwohl manche Generatoren hervorragende Eigenschaften im Reihen- Korrelationstest zeigen, diese total bei der Simulation von zufälligen Punkten in mehrdimensionalen Räumen versagen. In einem $n$- dimensionalen Raum würfelt man nacheinander die Koordinaten auf den $n$ Achsen. Sei die Zahlenfolge

\begin{displaymath}
\{ x_{1}, x_{2}, ......, x_{n}, x_{n+1}, ..... \}
\end{displaymath}

vorgelegt, dann sind die zufälligen Vektoren im $n$- dimensionalen Raum durch
$\displaystyle \vec{X}_{1}$ $\textstyle =$ $\displaystyle ( x_{1}, x_{2}, ....., x_{n} )$  
$\displaystyle \vec{X}_{2}$ $\textstyle =$ $\displaystyle ( x_{n+1}, x_{n+2}, ....., x_{2n} )$  
$\displaystyle \vec{X}_{3}$ $\textstyle =$ $\displaystyle ( x_{2n+1}, x_{2n+2}, ....., x_{3n} )$  
  $\textstyle .$    
  $\textstyle .$    
$\displaystyle \vec{X}_{m}$ $\textstyle =$ $\displaystyle (x_{(m-1)n+1}, x_{(m-1)n+2}, ....., x_{m n} )$  

gegeben. Hierbei zeigte sich, daß fast alle Generatoren bevorzugt Punkte auf bestimmten Hyperebenen erzeugen. Dieses Schwäche von Zufallszahlen- Generatoren können wir mit dem Reihen- Korrelationstest nicht aufzeigen.

Der Phasentest
Ein weiterer Test für die serielle Unabhängigkeit von Zufallszahlen ist der Phasentest. Gegeben sei eine Zahlenfolge

\begin{displaymath}
x_{i} = \{ x_{1}, x_{2}, x_{3}, ....., x_{n-1}, x_{n} \} .
\end{displaymath}

In dieser Zahlenfolge prüfen wir die Längen $l$ der monoton steigenden Unterfolgen:

\begin{displaymath}
x_{i-1} \geq x_{i} < x_{i+1} < ..... < x_{i+j} \geq x_{i+j+1} .
\end{displaymath}

Zum Beispiel hat die einziffrige Zahlenfolge

\begin{displaymath}
x_{i} = \{ 5, 4, 2, 4, 5, 5, 6, 9, 8, 6, 1, 3 \}
\end{displaymath}

die monoton steigenden Unterfolgen (auch Phasen genannt)

\begin{displaymath}
\{ 5 \}, \{ 4 \}, \{2, 4, 5 \}, \{5, 6, 9\}, \{ 8 \}, \{ 6 \}, \{1, 3\} ,
\end{displaymath}

und damit die Längen

\begin{displaymath}
l = 1, 1, 3, 3, 1, 1, 2.
\end{displaymath}

Aus relativ schwierigen mathematischen Rechnungen folgt der Satz: Die diskreten Häufigkeiten für die Anzahlen der Phasen mit der Länge $j$ sind für eine Zahlenfolge von $n$ unabhängigen Zufallszahlen durch
\begin{displaymath}
k_{j} = \frac{(n+1) (j^{2}+j-1)-(j+2)(j^{2} - j - 1)}{(j+2)! }
\end{displaymath} (15)

gegeben.

Für eine vorgegebene Zahlenfolge kann man die Unabhängigkeit dann mit Hilfe des $\chi^{2}$- Tests prüfen.

Das Testprogramm
Ein Java Applet, das die soeben diskutierten Tests durchführt, können Sie hier aufrufen. Zuerst müssen Sie einen der bereits vorprogrammierten Generatoren aufrufen. Wählen Sie dazu im Menu Generatoren zunächst den Generator JavaRandom auf. Dieses ist der Standard- Generator der Java Programmiersprache. Den Source Code können Sie mit dem Menu- Button Source einsehen. Es wird ein internes Fenster mit dem Listing eingeblendet. Danach wählen Sie im Menu Plots die für Sie interessanten Histogramme auf. Das Histogramm GleichVert zeigt die Gesamtverteilung. Die erwarteten Mittelwerte und Varianz dieser Verteilung sind.

$\displaystyle < x >$ $\textstyle =$ $\displaystyle \frac{1}{2}$  
$\displaystyle \sigma_{x} = \sqrt{< (x-<x>)^{2} > ) }$ $\textstyle =$ $\displaystyle \frac{1}{\sqrt{12}} = 0.288675$  

Nach Betätigung des Buttons Start werden nacheinander Sequenzen mit je $n=100$ Zufallszahlen erzeugt. Die Verteilung der Sequence wird im Plot Sequence gezeigt. Für jede dieser Sequencen wird ein $\chi^{2}$- Test durchgeführt sowie der Reihenkorrelationskoeffizient berechnet und in die Plots Chi2-Test und Korrelation eingetragen. Zusätzlich wird die Verteilung der Phasen bestimmt und in den Plot Phasentest eingetragen. In diesem Plot zeigen wir auch die theoretisch erwartete Verteilung, die in schwarzer Farbe dargestellt ist und einen Offset von 10 auf der horizontalen Skala hat.

Die weiteren wählbaren Generatoren sind der RN32, der bereits in Tabelle 1 angegeben war, sowie die Generatoren Ranmar, Ranlux und Ranecu, die wir später noch diskutieren werden. Als Beispiel eines schlechten Generators (BadRandom) haben wir in einem an sich guten 16-bit Generator einen Schreibfehler eingebaut (siehe Source Code). Dieses Beispiel zeigt, dass es offensichtlich in der Rekursionsrelation (1) auf die richtige Wahl der Parameter ankommt. Dieses wird uns im Rest dieses Kapitels beschäftigen. Nach einigen Minuten Rechenzeit haben wir die Ergebnisse in Tabelle 2 erhalten.


Tabelle: Test- Ergebnisse für einige Generatoren bei 200 Sequenzen mit je 100 Zufallszahlen.
  Theorie JavaRandom RN32 Ranmar BadRandom    
               
Gleichverteilungs- Test:              
$<x>$ 0.500 0.0497 0.500 0.498 0.455    
$\sigma_{x}$ 0.289 0.290 0.288 0.289 0.271    
               
Chi2-Test:              
$< \Theta' >$ 0. 0.078 0.092 0.036 1.679    
$\sigma_{\Theta'}$ 1. 0.990 0.925 0.980 0.686    
               
Korrelations- Test:              
$< \rho >$ -0.010 -0.006 -0.005 -0.002 0.288    
$\sigma_{\rho}$ 0.099 0.102 0.109 0.103 0.060    
               
Phasen- Test:              
Übereinstimmung   gut sehr gut gut schlecht    

Natürliche Zahlen
Für die folgenden Erörterungen benötigen wir noch einige Begriffe aus der Theorie der natürlichen Zahlen. Der wohl wichtigste ist der Begriff der Primzahl. Eine Zahl $p$ aus der Menge der natürlichen Zahlen $n>0$ heißt Primzahl, wenn $p$ außer durch sich selbst und der Zahl 1 durch keine andere natürliche Zahl ohne Rest teilbar ist. Bereits Euklid hat 300 v.Chr. festgestellt, daß es unendlich viele Primzahlen gibt. Unter Primzahlzwillingen versteht man zwei Primzahlen, die durch genau eine Nicht- Primzahl voneinander getrennt sind, z.B. 5 und 7 oder 17 und 19. Ein wichtiger Satz in diesem Zusammenhang ist der folgende
Satz: Jede Zahl zwischen zwei Primzahlzwillingen ist durch 2 und 3 und damit auch durch 6 teilbar.

Will man daher alle Primzahlzwillinge aufsuchen, muß man nur alle Vielfachen von 6 bilden und dann die vorherige und nachfolgende Zahl auf die Eigenschaft ,,Primzahl'' prüfen. Ob es unendlich viele Primzahlzwillinge gibt, ist bis heute nicht gezeigt worden.

Unter einer kanonischen Zerlegung einer natürlichen Zahl $n$ versteht man eine Zerlegung in Primzahl- Potenzen

\begin{displaymath}
n = p_{1}^{a_{1}} p_{2}^{a_{2}} \cdot \cdot \cdot p_{r}^{a_{r}}
\end{displaymath} (16)

mit positiven ganzzahligen Exponenten $a_{1}, a_{2}, ...., a_{r}$ und den Primzahlen $p_{1}, p_{2}, ...., p_{r}$. Diese Formel gibt sofort die Anzahl der Teiler von $n$ an. Jede Zahl

\begin{displaymath}
m = p_{1}^{b_{1}} p_{2}^{b_{2}} \cdot \cdot \cdot p_{r}^{b_{r}}
\end{displaymath}

mit

\begin{displaymath}
0 \leq b_{i} < a_{i}, \; \; \; i=1,2,...,r
\end{displaymath}

ist Teiler von $n$. Denn es ist dann

\begin{displaymath}
\frac{n}{m} = p_{1}^{a_{1}-b_{1}} p_{2}^{a_{2}-b_{2}} \cdot \cdot \cdot
p_{r}^{a_{r}-b_{r}}
\end{displaymath}

ganzzahlig. Damit ist es einfach zu erkennen, daß

\begin{displaymath}
\tau (n) = (a_{1}+1)(a_{2}+1) \cdot \cdot \cdot (a_{r}+1)
\end{displaymath}

die Anzahl der echten positiven Teiler von $n$ ist. Als echten Teiler bezeichnet man hierbei alle Teiler mit Ausnahme von 1 und der Zahl $n$ selbst.

Dividiert man eine beliebige natürliche Zahl $b$ durch eine beliebige andere natürliche Zahl $a$ , so erhält man im allgemeinen einen Rest. Der kleinste positive Rest $r$ ist durch die Formel

\begin{displaymath}
b = q a + r
\end{displaymath}

oder

\begin{displaymath}
\frac{b}{a} = q + \frac{r}{a}
\end{displaymath}

gegeben. Dieses kann man auch mit Hilfe der bereits definierten Modulo- Funktion schreiben als

\begin{displaymath}
r = mod(b,a).
\end{displaymath}

Sind $m$ positive ganze Zahlen $a_{1}, a_{2}, ...., a_{m}$ gegeben, so gibt es $n$ positive ganze Zahlen $t_{j}, j=1,2,...,n$, die alle gemeinsame Teiler der $a_{i}$ sind., d.h.

\begin{displaymath}
a_{i} = q_{ij} t_{j}.
\end{displaymath}

Unter diesen gibt es einen größten gemeinsamen Teiler (GGT) $t$. Als Erläuterung hierzu möge das bereits von Euklid entwickelte Verfahren zur Bestimmung des größten gemeinsamen Teilers (GGT) dienen. Sind einige der $a_{i}$ identisch gleich Null, so können diese Zahlen gestrichen werden, weil jede natürliche Zahl Teiler der Null ist. Wir nehmen also an, daß alle $a_{i} > 0, i=1,2,...,m$. Wir sortieren die Zahlen nun so, daß $a_{1}$ die kleinste Zahl der $a_{i}$ ist und dividieren alle restlichen Zahlen durch $a_{1}$:

\begin{displaymath}
a_{j} = a_{1} q_{j} +r_{j}, \; \; \; j=2,3,...,m.
\end{displaymath}

Hierbei ist jetzt

\begin{displaymath}
0 \leq r_{j} < a_{1}, \; \; \; j=2,3,...,m.
\end{displaymath}

Dadurch erreicht man, daß

\begin{displaymath}
GGT \{ a_{1}, a_{2}, ..., a_{m} \} = GGT \{ a_{1}, r_{2}, ..., r_{m} \}.
\end{displaymath}

Man verfährt nun mit der reduzierten Zahlenfolge $ \{ r_{2}, r_{3}, ..., r_{m} \}$ weiter, sucht die kleinste der Zahlen $r_{i}, i=2,3,...,m$, und bildet wie oben die kleinsten positiven Reste, Das Verfahren bricht nach endlich vielen Schritten ab und liefert den $GGT \{ a_{1}, a_{2}, ..., a_{m} \}$.

Als praktisches Beispiel bestimmen wir den größten gemeinsamen Teiler der Zahlenfolge $ \{ 36, 0, 96, 120 \} $. Nach Streichen der Null und Division durch 36 ergibt der erste Schritt

\begin{displaymath}
GGT \{ 36,96,120 \} = GGT \{ 36, 24, 12 \}.
\end{displaymath}

Daher, nach Umstellen der 12,

\begin{displaymath}
GGT \{ 36, 96, 120 \} = GGT \{ 12, 24, 36 \} = 12.
\end{displaymath}

Der größte gemeinsame Teiler ist also 12, alle positiven Teiler der obigen Zahlenfolge sind durch $ \{ 1, 2, 3, 4, 6, 12 \}$ gegeben.

Multiplikative Kongruenzmethode
Wir kommen jetzt auf die Diskussion der allgemeinen Rekursionsrelation (1) zurück, die wir der Wichtigkeit halber noch einmal notieren:

\begin{displaymath}
x_{n+1} = mod( a_{0} x_{n} + a_{1} x_{n-1} + ..... + a_{j} x_{n-j} + b, M) ,
\end{displaymath}

mit nichtnegativen ganzen Zahlen $M, a_{0}, a_{1}, ..., a_{j}, b$ und $x_{0}, x_{1}, ... , x_{j}$. Wir setzen voraus, daß der größte gemeinsame Teiler aller Parameter gleich 1 ist,
\begin{displaymath}
GGT \{ x_{0}, x_{1}, ..., x_{j}, b, M, a_{0}, a_{1}, ..., a_{j} \} = 1.
\end{displaymath} (17)

Die Glieder der Folge können nicht größer als $M$ werden, notwendigerweise muß es also periodische Wiederholungen der Zahlenfolge geben. Jede Periode umfaßt $L(M)$ Zahlen, $L(M)$ selbst heißt Periodenlänge. Die Zahlenfolgen mit dem Modul $M$ bei sonst gleichen Parametern bezeichnen wir im folgenden mit Z(M). Es gilt der
Satz: Seien $M_{1}$ und $M_{2}$ zwei teilerfremde positive ganze Zahlen mit den Periodenlängen $L(M_{1})$ und $L(M_{2})$. Dann ist

\begin{displaymath}
L(M_{1} M_{2}) = KGV \{ L(M_{1}, L(M_{2}) \}.
\end{displaymath} (18)

Hierbei bedeutet $KGV$ das kleinste gemeinsame Vielfache der Zahlen in der geschweiften Klammer. Die Zahl $M$ läßt sich immer als kanonische Primzahlzerlegung schreiben:

\begin{displaymath}
M = p_{1}^{a_{1}} p_{2}^{a_{2}} \cdot \cdot \cdot p_{r}^{a_{r}} .
\end{displaymath}

Da Primzahlen nach Definition immer teilerfremd zueinander sind, gilt also auch
\begin{displaymath}
L(M) = KGV \{ L(p_{1}^{a_{1}}), L(p_{2}^{a_{2}}), ..., L(p_{r}^{a_{r}}) \}.
\end{displaymath} (19)

Zwei weitere Sätze sind für das folgende wichtig:
Satz: Wenn $Z(p^{m})$ die Periodenlänge $L(p^{m})$ besitzt, mit Primzahl $p$ und natürlicher Zahl $m$, dann ist entweder

\begin{displaymath}
L(p^{m+1}) = L(p) \end{displaymath} (20)

oder
\begin{displaymath}
L(p^{m+1}) = p L(p).
\end{displaymath} (21)

Satz: Ist $L(p^{m+1}) = p L(p)$, so ist auch
\begin{displaymath}
L(p^{m+l}) = p^{l} L(p^{m}),
\end{displaymath} (22)

mit natürlichen Zahlen $l$ und $m$.

Wir gehen jetzt auf den wichtigsten und gleichzeitig einfachsten Spezialfall ein, nämlich auf die spezielle multiplikative Kongruenzmethode mit Modul $2^{m}$:

\begin{displaymath}
x_{n+1} = mod( k x_{n}, 2^{m}).
\end{displaymath} (23)

Wir diskutieren im folgenden die einfachen Fälle $m=1,2,3,4$ mit dem Startwert $x_{0} = 1$. Für $m=1$, d.h. $M=2$, muß $k=1$ sein. Die Zahlenfolge besteht nur aus der Zahl 1, die Periodenlänge ist $L(2^{1}) = 1$. Für $m=2$, d.h. $M=4$, kann $k=1$ oder $k=3$ gesetzt werden. Im ersten Fall ergibt sich wiederum eine Zahlenfolge mit nur Einsen und der Periodenlänge 1, im zweiten Fall abwechselnd 3 und 1.

\begin{displaymath}
X(2^{2}) = \lbrace \begin{tabular}{ccc} \{ 1,1,1,..... \} & f''ur &
k=1 \\ \{ 3,1,3,1,..... \} & f''ur & k=3 \end{tabular}
\end{displaymath}


\begin{displaymath}
L(2^{2}) = \lbrace \begin{tabular}{ccc} 1 & f''ur & k=1 \\
2 & f''ur & k=3 \end{tabular}
\end{displaymath}

Für $m=3$, d.h. $M=8$, kann $k=1,3,5,7$ gewählt werden. Untersuchung der Rekursion mit $x_{0} = 1$ ergibt

\begin{displaymath}
X(2^{3}) = \lbrace \begin{tabular}{ccc} \{ 1,1,1,1,..... \} ...
....... \} & & k=5 \\
\{ 7,1,7,1,..... \} & & k=7 \end{tabular}
\end{displaymath}


\begin{displaymath}
L(2^{3}) = \lbrace \begin{tabular}{ccc} 1 & f''ur & k=1 \\
2 & & k=3 \\ 2 & & k=5 \\ 2 & & k=7 \end{tabular}
\end{displaymath}

Würden wir mit anderen Startwerten $x_{0}$ beginnen, die natürlich die Bedingung (17) erfüllem müssen, bekämen wir zwar andere Zahlenfolgen, an der Periodenlänge ändert sich aber nichts.

Interessanter wird es jetzt bei $m=4$, d.h. $M=16$. Erlaubte $k$- Werte sind für $x_{0} = 1$ die Zahlen 1, 3, 5, 7, 9, 11, 13 und 15. Wie der Leser leicht nachprüfen kann, ergeben sich für die Periodenlängen

\begin{displaymath}
L(2^{4}) = \lbrace \begin{tabular}{ccc} 1 & f''ur & k=1 \\
...
... k=9 \\ 4 & & k=11 \\
4 & & k=13 \\ 2 & & k=15 \end{tabular}
\end{displaymath}

Man sieht also, daß bestimmte $k$- Werte vor anderen durch eine maximale Periodenlänge ausgezeichnet sind. In Verallgemeinerung der soeben diskutierten Beispiele gilt der
Satz: Die maximale Periodenlänge der Rekursion (23) ist für $m \geq 3$ durch $2^{m-2}$ gegeben. Diese Länge wird für $k= 3 + 8n$ und $k= 5 + 8n$ mit $n \geq 0$ und für alle ungeraden $x_{0}$ erreicht.

Wir kommen jetzt noch einmal auf den bereits am Anfang dieses Kapitels diskutierten Generator

$\displaystyle x_{n+1}$ $\textstyle =$ $\displaystyle mod(k x_{n}, 2^{15}) \; f''ur \; 16-bit \; Rechner$  
$\displaystyle x_{n+1}$ $\textstyle =$ $\displaystyle mod(k x_{n}, 2^{31}) \; f''ur \; 32-bit \; Rechner$  

zurück. Für die maximale Periodenlänge erhält man $L(2^{15}) = 2^{13} = 8192$ b.z.w. $L(2^{31}) = 2^{29} = 536870912$. Diese maximalen Periodenlängen erhält man für $k= 3 + 8n$ und $k= 5 + 8n$ mit $n \geq 0$ und einem beliebigen, ungeraden Wert für $x_{0}$. Welcher dieser $k$- Werte den besten Generator liefert, ist bis heute ungelöst. Hier ist man auf umfangreiche Tests angewiesen. Getestet und für ,,gut'' befunden wurde für 32-bit Rechner der Wert $k= 5 + 8n$. Die Zahl $n$ sollte so gewählt werden. daß $k$ größenordnungsmäßig in der Nähe von $\sqrt{M}$ liegt. Ein bekannter Generator (RN32) für 32-bit Rechner benutzt die Zahl $k=69069$ (siehe das Applet in diesem Kapitel).

Eine Erweiterung dieses Generators ergibt sich durch Hinzunahme des Summanden $b$:

\begin{displaymath}
x_{n+1} = mod( k x_{n} + b, 2^{m}).
\end{displaymath} (24)

Über diesen Generator sind umfangreiche Arbeiten veröffentlicht worden. Ohne Beweis geben wir das wichtigste Resultat: Die maximale Periodenlänge des Generators (24) beträgt $L(2^{m}) = 2^{m}$. Diese Länge wird erreicht für $k = 1 + 4n$ und $b = 1 + 2l$, mit beliebigen $n,l \geq 0$. Eine weitere Einschränkung der Parameter ist auch hier nicht möglich, man ist auf umfangreiche Tests angewiesen.

Zusammenstellung der wichtigsten Kongruenzmethoden
Im folgenden geben wir eine Zusammenstellung der am häufigsten benutzten Parameter in professionellen Anwendungen. Alle Generatoren verwenden keine additive Konstante, d.h. $b=0$.
1. $k=23$, $M=10^{8}+1$: Dieses sind die Parameter, die Lehmer [ ] in seiner Originalarbeit über multiplikative Kongruenzmethoden vorgeschlagen hatte.
2. $k=65539$, $M=2^{29}$: Dieses ist ein Produkt der Computerfirma IBM für ihre /360- Serie. Die Parameter wurden optimiert, um einen möglichst kleinen Reihen- Korrelationskoeffizienten zu erhalten. Unglücklicherweise hat dieser Generator verheerende Korrelationen höherer Ordnung im $n$- dimensionalen Raum. Der Grund hierfür scheint zu sein, daß der $k$- Wert nicht die Relation

\begin{displaymath}
k= 5 + 8n, \; \; \; n \geq 0 ,
\end{displaymath}

sondern die Relation

\begin{displaymath}
k = 3 + 8n, \; \; \; n \geq 0
\end{displaymath}

erfüllt. Dieses war einer der Tests, die gezeigt hatten, daß diese zweite Bedingung für einen guten Generator vermieden werden sollte.
3. $k=69069$, $M = 2^{32}$: Dieses ist der am häufigsten benutzte Generator der Jahre 1980 bis 1990 und unter dem Namen RNDM bekannt. Gegenüber dem von uns ausgiebig diskutierten Generator RNDM unterscheidet er sich durch eine doppelt so große Periodenlänge von $L(2^{32}) = 2^{30} = 1073741824$. Da in diesem Generator alle 32 bits, einschließlich des Vorzeichenbits, benutzt werden, sind spezielle Programmtechniken notwendig (siehe Aufgaben).
4. $k=16807$, $M=2^{31}-1$: Ebenfalls ein Produkt der Firma IBM.
5. $k=1664525$, $M = 2^{32}$: Der hier verwendete Multiplier $k$ wurde von Knuth als der beste für den Modul $2^{32}$ vorgeschlagen.
6. $k=742938285$, $M=2^{31}-1$: Ein von L'Ecuyer vorgeschlagener Generator. In seiner Arbeit weist der Autor nach, daß der angegebene Multiplier der beste für den Modul $2^{31}-1$ ist.

Fibonacci- Generatoren
Wir hatten im Rahmen der Kombinatorik bereits eine Fibonacci- Sequenz kennengelernt:

\begin{displaymath}
F(k+1) = F(k) + F(k-1).
\end{displaymath}

In Verallgemeinerung dieser Sequenz nennt man alle Generatoren, die auf der Rekursion
\begin{displaymath}
x_{n} = mod( x_{n-p} \bigodot x_{n-q}, 2^{m})
\end{displaymath} (25)

basieren, Fibonacci- Generatoren. Die Verknüpfung $\bigodot$ bedeutet hierbei irgendeine binäre oder logische Operation zwischen den im Dualcode dargestellten Zahlen. Man kann aber auch für $\bigodot$ die normale Addition oder Subtraktion verwenden. Für bestimmte Wahl von $p$ und $q$ kann man zeigen, daß die maximale Periodenlänge durch

\begin{displaymath}
L = (2^{p}-1) 2^{m-1}
\end{displaymath}

gegeben ist.

Die Fibonacci- Generatoren sind zwar in unzähligen Arbeiten theoretisch untersucht worden, eine klare Antwort ist jedoch (zumindest für den Autor dieses Tutorials) nicht erkennbar. Es muß jedoch schon jetzt erwähnt werden, daß der beste Generator, der gegenwärtig auf dem Markt ist, ein Fibonacci- Generator ist (siehe später).

Vergrößerung der Periodenlänge
Alle von uns bisher diskutierten Generatoren hatten eine Periodenlänge in der Größenordnung ihres Moduls. Mit anderen Worten, nachdem ungefähr alle Zahlen aus dem Zahlenbereich $0 < x < M$ einmal generiert worden waren, wiederholte sich die Zahlenfolge in exakt gleicher Weise. Wenn also

\begin{displaymath}
Z_{1}(M) = \{ x_{1}, x_{2}, ....., x_{L(M)} \}
\end{displaymath}

die generierte Zahlenfolge beim ersten Durchlauf war, so ist

\begin{displaymath}
Z_{2}(M) = \{ x_{L(M)+1}, x_{L(M)+2}, ....., x_{2L(M)}
\end{displaymath}

mit
$\displaystyle x_{1}$ $\textstyle =$ $\displaystyle x_{L(M)+1}$  
$\displaystyle x_{2}$ $\textstyle =$ $\displaystyle x_{L(M)+2}$  
  $\textstyle .$    
  $\textstyle .$    
$\displaystyle x_{L(M)}$ $\textstyle =$ $\displaystyle x_{2L(M)}$  

die Zahlenfolge beim zweiten Durchlauf. Dieses ist die exakte Definition der Periodenlänge $L(M)$.

Man kann nun aber die Periodenlänge in einfachster Weise vergrößern, indem man in der zweiten Zahlenfolge eine Permutation der Zahlen erzeugt. Die verschiedenen Permutationen dürfen keine Korrelation zeigen, d.h. die Permutationen müssen mit einem zweiten, unabhängig vom ersten Generator laufenden, Zufallszahlen- Generator ausgewürfelt werden. Gemäß der 1. Grundaufgabe der Kombinatorik kann man dann, zumindest im Prinzip, eine Vergrößerung der Periodenlänge um den Faktor $L(M)! = 1 \cdot 2 \cdot \cdot \cdot L(M)$ erwarten. Wir wollen in diesem Zusammenhang zwei Generatoren diskutieren, die heutzutage (2004) als das beste angegeben werden, was es auf dem Markt gibt, nämlich der Generator RANECU von L'Ecuyer [ ] und der Generator RANMAR von Marsaglia [ ]. Beide Generatoren arbeiten nicht, wie bisher alle diskutierten Generatoren, mit der Überlauf- Arithmetik. Alle Rechenoperationen bleiben im Zahlenbereich des Rechners. Es gibt Versionen für 16-bit und 32-bit Rechner. Wir diskutieren hier die 32-bit Versionen.

Der RANECU Generator
L'Ecuyer beschreibt einen Generator, der aus 2 (b.z.w. 3) multiplikativen Kongruenzmethoden besteht, aber so, daß alle Rechenoperationen innerhalb des Zahlenbereichs des 32-bit (b.z.w. 16-bit) Rechners bleiben. Die Begründung dieses Generators ist nicht einfach zu verstehen.

Die Parameter in den beiden kombinierten Kongruenzen sind $k_{1}=40014$, $M_{1}=2147483563$ und $k_{2}=40692$, $M_{2}=2147483399$ und sollten nicht verändert werden. Man erkennt, daß keine dieser Zahlen etwas mit unseren bisherigen Kriterien für ,,gute'' Zufallszahlen- Generatoren zu tun hat. Trotzdem hat sich gezeigt, daß der Generator hervorragende Eigenschaften bezüglich der Gleichverteilung, der seriellen sowie der Unabhängigkeit höherer Ordnung besitzt. Die Periodenlänge ist $\approx 2^{60} \approx 10^{18}$.

Der RANMAR Generator
Dieser Generator von Marsaglia und Zamann basiert auf zwei Fibonacci- Sequenzen, eine mit einer Subtraktion und Modulo- Funktion, die andere auf einer einfachen Subtraktion. Auch in diesem Fall wollen wir nicht die Einzelheiten diskutieren, es geht über den Rahmen dieses Tutorials hinaus. Wichtig ist lediglich, daß die Startwerte IR und IS aus nicht mehr als 5 dezimalen Ziffern bestehen dürfen und die Relationen

$\displaystyle 0$ $\textstyle \leq$ $\displaystyle IR \leq 31328,$  
$\displaystyle 0$ $\textstyle \leq$ $\displaystyle IS \leq 30081$  

erfüllen müssen. Der Generator muß mit einem einmaligen Aufruf der Unterroutine RMARIN initialisiert werden. Der RANMAR Generator arbeitet ebenfalls im reinen 32-bit Zahlenbereich des Rechners und erzeugt eine Zahlenfolge mit der unglaublich großen Periodenlänge von $2^{144} \approx 10^{43}$.

Aufgaben

Aufgabe 1: Untersuchen Sie den System- Generator Ihres Rechners. Suchen Sie dazu im Source Code der Java Distribution den entscheidenden Code heraus, der beim einfache Aufruf Math.random() interpretiert wird. Welcher Algorithmus verbirgt sich hinter diesem Generator ?

Aufgabe 2: Diskutieren Sie die Periodenlänge des Generators

\begin{displaymath}
x_{n+1} = mod( k x_{n}, 2^{m}-1).
\end{displaymath}

Gehen Sie hierbei von $m$= 2, 3 und 4 aus und verwenden Sie die vollständige Induktion von $m$ auf $m+1$.
a) Welche $k$- Werte ergeben eine maximale Periodenlänge?
b) Schreiben Sie den Java Code für $m=31$ und $k=16807$.
c) Untersuchen Sie diesen Generator auf Gleichverteilung, serielle Korrelationen, $\chi^{2}$- Test und Phasentest.





Harm Fesefeldt
2005-03-29