Xorshift128+ ist ein moderner Algorithmus zur schnellen Erzeugung
von Zufallszahlen relativ hoher Güte, der die wichtigsten statistischen
Tests besteht. Der Generator soll eine Periodenlänge von
2128−1 besitzen und »BigCrush« aus der Test-Suite
»TestU01« bestehen. Zu beachten ist, dass es sich keinesfalls um einen
kryptografischen Generator handelt. Die folgende Implementierung
ist wie in Vignas Artikel [1] beschrieben, mit der dort als
günstig beschriebenen
Belegung a=23, b=17, c=26 und
wird so auch im Firefox-Quellcode [2] und in Wikipedia [3][4]
aufgeführt. Der interne Zustand state darf beliebig sein,
jedoch nicht (0, 0).
Da die niederwertigen 32 Bit eines Zufallswertes von
Xorshift128+ einige statistische Tests nicht bestehen sollen,
extrahiert man für rand_u32 besser die höherwertigen
32 Bit.
pub struct RNG {
state: (u64, u64)
}
impl RNG {
pub fn from_seed(seed: u64) -> Self {
Self {state: (
seed ^ 0xf4dbdf2183dcefb7, // [crc32(b"0"), crc32(b"1")]
seed ^ 0x1ad5be0d6dd28e9b // [crc32(b"2"), crc32(b"3")]
)}
}
pub fn rand_u64(&mut self) -> u64 {
let (mut x, y) = self.state;
self.state.0 = y;
x ^= x << 23;
self.state.1 = x ^ y ^ (x >> 17) ^ (y >> 26);
self.state.1.wrapping_add(y)
}
pub fn rand_u32(&mut self) -> u32 {
(self.rand_u64() >> 32) as u32
}
}
fn main() {
let mut rng = RNG::from_seed(0);
for _ in 0..10 {
println!("{:016x}", rng.rand_u64());
}
}
Auch das Zeroland von Xorshift128+ ist verschwindend klein. Als Zeroland bezeichnet man Bereiche der Zustandstrajektorie, die schlechte Zufallszahlen liefern. Bei vielen Generatoren ist das damit verbunden, dass viele Bits null sind, daher die Bezeichnung Zeroland. Der Mersenne-Twister soll z. B. im schlechtesten Fall, bei dem nur ein einziges Bit gesetzt ist, eine Aufwärmphase von über 700 000 Iterationen benötigen. Das kann man als problematisch sehen, denn der Generator könnte das Zeroland, wenn auch mit einer geringen Wahrscheinlichkeit, je nach Saat zufällig durchqueren.
[1] Sebastiano Vigna: »Further scramblings of Marsaglia's xorshift generators«. Journal of Computational and Applied Mathematics (Mai 2017). arXiv:1404.0390. doi:10.1016/j.cam.2016.11.006.
[2] »mfbt/XorShift128PlusRNG.h« in »gecko-dev«, abgerufen am 8. Feb. 2020. Gecko ist die aktuelle Engine von Firefox.
[3] »Xorshift«. Englische Wikipedia, abgerufen am 8. Feb. 2020.
[4] »Xorshift«. Deutsche Wikipedia, abgerufen am 8. Feb. 2020.
Zur Erzeugung von gleichverteilten Zufallszahlen
x ∈ [i..j) genügt es, Zahlen
u ∈ [0..m) zu erzeugen. Man kann nämlich
m := j-i
setzen und dann i zur erzeugten Zahl u
addieren. Ausgehend von
0 ≤ u < j-i,
bringt Addition von i auf beiden Seiten der jeweiligen
Ungleichung die äquivalenten Ungleichungen
i ≤ u+i < j-i+i.
Mit x := u+i ergibt sich
i ≤ x < j.
Naive Methode für Zufallszahlen aus [0..m):
pub fn rand_bounded(&mut self, m: u32) -> u32 {
self.rand_u32() % m
}
Die mit dieser Methode erzeugten Zahlen sind jedoch nicht
gleichverteilt. Dazu müssten alle Restklassen modulo m
gleich groß sein, dies ist jedoch nicht der Fall. Die mit dieser
Methode erzeugte Zahl u ist kleinster Repräsentant seiner
Restklasse. Die Restklasse von u ist das Urbild von
u bezüglich der Operation f(r) := r%m.
Man nennt das Urbild auch Faser, weil der Definitionsbereich
[0..232) in disjunkte Urbilder zerlegt wird.
Wie bekommen wir nun auf [0..m) gleichverteilte
Zufallszahlen? Eine wichtiges Grundverfahren ist die
Verwerfungsmethode. Das geht so: Erzeugen wir zunächst gleichverteilte
Zufallszahlen aus einem größeren Bereich und verwerfen davon alle
außerhalb [0..m), verbleiben die restlichen Zufallszahlen
gleichverteilt.
Da die Zufallzahlen jedoch aus [0..232)
entstammen, enstehen bei der Methode für ein kleines m
viel zu viele Verwerfungen, was zu einem enormen Rechenaufwand
führt. Die Idee: Schneiden wir Bits vom Ergebnis von
rand_u32() ab, verbleiben die übrigen Bits
gleichverteilt. Somit sind wir im Besitz von gleichverteilten
Zufallszahlen aus [0..2k)
zu jeder Zweierpotenz 2k. Zu m
suchen wir die nächste Zweiterpotenz
2k ≥ m.
Infolge sinkt der Rechenaufwand bei der Verwerfungsmethode
drastisch ab. Die Standardbibliothek besitzt sogar die Funktion
next_power_of_two.
pub fn rand_bounded(&mut self, m: u32) -> u32 {
let mask = m.next_power_of_two() - 1;
loop {
let x = mask & self.rand_u32();
if x < m {return x;}
}
}
Für m > 231 kommt es
bei der Berechnung der nächsten Zweiterpotenz zu einem Überlauf.
Korrektes Verhalten für diesen Fall wird so berücksichtigt:
let mask = m.wrapping_next_power_of_two().wrapping_sub(1);
Außerdem ist die wiederholte Berechnung der Bitmaske
bei festem m eigentlich nicht notwendig.
Partielle Applikation gestattet die Abtrennung dieser Berechnung.
pub fn rand_bounded<'a>(&'a mut self, m: u32)
-> impl 'a + FnMut()->u32
{
let mask = m.next_power_of_two() - 1;
move || {
loop {
let x = mask & self.rand_u32();
if x < m {return x;}
}
}
}
Dividiert man gleichverteilte Zufallszahlen aus
[0..232) durch 232,
resultiert dies näherungsweise in gleichverteilten
Zufallszahlen aus [0,1).
pub fn rand_float(&mut self) -> f64 {
f64::from(self.rand_u32())*2.3283064365386963E-10
}
Dies liefert 32 Bit Genauigkeit, – umgerechnet sind das
ca. neun dezimale Nachkommastellen. Warum nicht eine Zahl
aus [0..264) durch 264
teilen? Dazu muss man bedenken, dass bei f64
die Mantisse 53 Bit lang ist. Die Mantisse kann die erzeugten
64 Bit folglich nicht aufnehmen. Dies könnte zu einer
Verzerrung der Verteilung führen.
Bei u32 haben wir den großen Vorteil, eine Berechnung
über den gesamten Definitionsbereich laufen lassen zu können.
Bei u64 geht dies aufgrund der schieren Größe nicht
mehr. Erhält die Zuordnung
fn float(x: u32) -> f64 {
f64::from(x)*2.3283064365386963E-10
}
wirklich die Gleichverteilung, müssen zwei
gleich lange Teilintervalle von [0,1) beim Durchlaufen
aller Werte aus [0..232) gleich häufig getroffen
werden. Oder anders ausgedrückt, die Urbilder von gleich langen
Teilintervallen müssen gleich groß sein. Diese Experimente können
wir machen, ohne die genaue Implementation von f64::from
zu kennen.
Gegeben sei eine Verteilungsfunktion F. Ist
u eine Zufallszahl bezüglich Gleichverteilung auf
[0,1], so ist F−1(u)
eine Zufallszahl bezüglich Verteilungsfunktion F.
Hierbei ist F−1 die Umkehrfunktion
von F. Dieses Verfahren zur Erzeugung von Zufallszahlen
einer bestimmten Verteilung nennt sich Inversionsmethode.
Nun ist die Umkehrfunktion nicht bei jeder Verteilungsfunktion einfach bestimmbar. Wir sind aber gierig und wollen doch irgendwie an die Zufallszahlen herankommen. Den Schwierigkeiten aus dem Weg geht die Berechnung mit dem Bisektionsverfahren.
// Bestimme y ∈ [a, b], so dass x = f(y).
fn bisection(f: impl Fn(f64)->f64,
x: f64, mut a: f64, mut b: f64
) -> f64
{
for _ in 0..80 {
let m = 0.5*(a + b);
if f(m) - x < 0.0 {a = m;} else {b = m;}
}
a
}
// Zufallszahlengenerator zu einer Verteilungsfunktion,
// engl. cumulative distribution function.
fn rng_from_cdf(mut rng: RNG,
a: f64, b: f64, cdf: impl Fn(f64)->f64
) -> impl FnMut()->f64
{
move || bisection(&cdf, rng.rand_float(), a, b)
}
Als Beispiel diene die Exponentialverteilung mit Erwartungswert
1/λ.
fn rng_from_exp(rng: RNG, lambda: f64) -> impl FnMut()->f64 {
rng_from_cdf(rng, 0.0, 100.0/lambda, move |x|
if x <= 0.0 {0.0} else {1.0 - f64::exp(-lambda*x)})
}
fn main() {
let rng = RNG::from_seed(0);
let mut rand = rng_from_exp(rng, 1.0);
for _ in 0..10 {
println!("{:.6}", rand());
}
}
Oft möchte man Zufallszahlen in Form eines Iterators haben, wofür
from_fn nützlich ist.
pub fn iter_u32<'a>(&'a mut self)
-> impl 'a + Iterator<Item = u32>
{
std::iter::from_fn(move || Some(self.rand_u32()))
}
Wir können Methoden auch als gewöhnliche Funktionen aufrufen. Benutzung eines Funktionenzeigers lässt eine Anpassung mit generischer Funktionalität zu.
pub fn iter<'a, T: 'static>(&'a mut self, rand: fn(&mut Self) -> T)
-> impl 'a + Iterator<Item = T>
{
std::iter::from_fn(move || Some(rand(self)))
}
Die Benutzung schaut so aus:
fn main() {
let mut rng = RNG::from_seed(0);
for x in rng.iter(RNG::rand_u32).take(10) {
println!("{:08x}", x);
}
}
Es tut sich noch die Frage auf, wie die Mischung eines Arrays am besten vonstatten geht. Gemeint ist, wie sich zu einer Länge eine zufällige Permutation der Elemente erzeugen lässt. Dies wird erreicht durch die durstenfeldsche Variante des Fisher-Yates-Algorithmus:
fn shuffle<T>(&mut self, a: &mut [T]) {
for i in (1..a.len()).rev() {
let j = self.rand_bounded_usize(i + 1);
a.swap(i, j);
}
}
Man kann beweisen, dass die mit diesem Algorithmus erzeugten Permutationen gleichverteilt zufällig sind, sofern dies für die benutzten Zufallszahlen gilt.