Implementere Mersenne Twister

Introduksjon

Python sin innebygde random-modul bruker generatoren Mersenne Twister (MT19937). Dersom dere er ferdig med den andre oppgaven og ønsker å lage en robust generator kan dere prøve dere på denne. Implementer MersenneTwister-klassen i mersenne.py, og kontroller at den gir de samme tallene som Python. Dere kan bruke plott.py for å sjekke om denne utgaven har de samme egenskapene som LCG-generatoren.

Kort om Mersenne Twister

En LCG har bare ett tall som tilstand. Det gjør den rask, men også lett å se mønstre i. Mersenne Twister har i stedet en tilstand på 624 tall på 32 bit hver. Der perioden til LCG er maksimalt størrelsen på modulusen, har denne en periode på 2^19937 - 1 steg (et Mersenne-primtall, derav navnet). Generatoren fungerer slikt:

1. Initialisering: Frøet fyller tilstanden mt[0..623], hvor hvert tall regnes ut fra det forrige.
2. Vri (twist): Når alle de 624 tallene er laget opp, regnes hele tilstanden ut på nytt ved å blande hvert tall med to andre i listen.
3. Temper: Hvert tall som hentes ut, blandes litt med bitskift og xor for å forbedre fordelingen av bitene.

Dere trenger bare disse operasjonene på heltall:

Konstanter

N = 624
M = 397
A = 0x9908B0DF
F = 1812433253
OVRE_MASKE  = 0x80000000   # bit 31
NEDRE_MASKE = 0x7FFFFFFF   # bit 0-30

Algoritmen

Initialisering med frø froe:

mt[0] = froe
mt[i] = (F * (mt[i-1] xor (mt[i-1] >> 30)) + i)  % 2^32      for i = 1 ... 623

Vri (kjøres når alle 624 tall er brukt):

for i = 0 ... 623:
    y     = (mt[i] & OVRE_MASKE) | (mt[(i+1) mod 624] & NEDRE_MASKE)
    mt[i] = mt[(i+M) mod 624] xor (y >> 1)
    hvis y er oddetall: mt[i] = mt[i] xor A

Temper av neste tall y = mt[indeks]:

y = y xor (y >> 11)
y = y xor ((y << 7)  & 0x9D2C5680)
y = y xor ((y << 15) & 0xEFC60000)
y = y xor (y >> 18)

Klassegrensesnitt

from mersenne import MersenneTwister

generator = MersenneTwister(5489)
print(generator.rand())     # 3499211612
print(generator.rand())     # 581869302
print(generator.random())
print(generator.randint(1, 6))

Del 1: Tilstand og initialisering

1. Implementer __init__ med initialiseringen over. Husk å kappe til 32 bit etter hver utregning.
2. Sett self.indeks = N. Det betyr at tilstanden er brukt opp, og rand vil vri den før første tall hentes.
3. Kontroller at mt[1] blir 1 og mt[2] blir 1812433255 når frøet er 0.

Del 2: Vri og temper

1. Implementer _vri. Merk at tilstanden endres på stedet, så tall tidlig i listen er allerede oppdatert når i kommer nær slutten.
2. Implementer rand: vri om nødvendig, hent mt[indeks], temper, øk indeks og returner.
3. Kjør med frø 5489. De tre første tallene skal være 3499211612, 581869302 og 3890346734. Tall nummer 10000 skal være 4123659995.
4. Med frø 1 skal de to første tallene være 1791095845 og 4282876139.

Del 3: Flyttall og heltall

1. Implementer random ved å dele på self.modulus.
2. Implementer randint som i LCG. Hvilket problem får dere med de laveste bitene her sammenlignet med LCG?
3. Kjør plott_histogram_og_par(generator) og plott_tripler(generator) med MersenneTwister i stedet for LCG. Ser dere noe struktur?

Del 4: Sammenlign med Python

Python lar dere hente tilstanden til Random-klassen med getstate(). Tilstanden er en tuppel med 625 tall: de 624 første er selve tilstandslisten mt, og det siste er indeksen. Dere kan derfor kopiere tilstanden inn i deres egen generator og sjekke at den er riktig, uten å kjenne Pythons frøalgoritme.

import random
from mersenne import MersenneTwister, N      # N = 624, lengden på tilstandslisten

ekte = random.Random(42)
versjon, tilstand, _ = ekte.getstate()      # tilstand er 624 tall + indeksen
mt = MersenneTwister(0)
mt.mt = list(tilstand[:N])                  # kopier tilstandslisten
mt.indeks = tilstand[N]                     # kopier indeksen

for _ in range(1000):
    assert mt.rand() == ekte.getrandbits(32)

1. Kjør koden over. Hvis det ikke fungerer, sammenlign mt.mt og tilstand etter første vri.
2. Python bruker en annen algoritme for å sette frøverdien, så vi kan dessverre ikke bruke MersenneTwister(42) for å få samme tallfølge som random.Random(42). Hvordan Python gjør det, ser dere på i den første utvidelsen.

Tips

• Bruk & 0xFFFFFFFF etter multiplikasjon og skift til venstre. Dette kapper tallene til 32-bit. Skift til høyre trenger det ikke.
• Bruk % N for å holde indeksene innenfor listen: mt[(i + 1) % N]. N er størrelsen til mt-listen.
• Skriv ut tilstanden med print(mt.mt[:5]) hvis noe ikke stemmer. Sjekk det første tallet som avviker.
• Vri ikke hele tilstanden før den er brukt opp. Det er indeks-sjekken i rand som styrer dette.

Utvidelser

• Python bruker ikke frøalgoritmen over, men init_by_array, som tar en liste med frø. Implementer den (se Matsumoto og Nishimuras referansekode, mt19937ar.c) og start med init_genrand(19650218). Da skal MersenneTwister(42) gi de samme tallene som random.Random(42).getrandbits(32).
• Python lager flyttall med 53 bit på denne måten: a = rand() >> 5, b = rand() >> 6, (a * 67108864 + b) / 9007199254740992. Implementer random slik og sammenlign med random.Random(42).random().
• Tilstanden kan gjenskapes fra 624 påfølgende tallene. Hent 624 tall fra en generator, avtemper hvert tall (alle tempersteg kan snus) og lag en ny generator som forutsier alle senere tall. Hva sier dette om bruk av Mersenne Twister til passord og kryptografi?
• Mål hvor lang tid 1 million kall tar for LCG, MersenneTwister og Pythons random. Hvorfor er Pythons egen så mye raskere?
• Sammenlign med PCG eller xorshift, som er raskere og har mindre tilstand. Hva taper dere?

Tips: For å måle hvor lang tid noe tar i Python bruker man typisk time.perf_counter().

import time # En modul for å håndtere tidsobjekter. Inneholder mange andre nyttige funskjoner.
def mal_tid(funksjon, n=1_000_000):
    """Returnerer hvor mange sekunder n kall til funksjonen tar. Brukes til å sammenligne generatorene."""
    start = time.perf_counter() # Lagrer starttidspunktet
    for _ in range(n):
        funksjon() # Gjør n-kall på en funksjon
    return time.perf_counter() - start # Returnerer tidsdeltaet fra start til slutt