Serwis Edukacyjny
w I-LO w Tarnowie
obrazek

Materiały dla uczniów liceum

  Wyjście       Spis treści       Wstecz       Dalej  

Autor artykułu: mgr Jerzy Wałaszek

©2024 mgr Jerzy Wałaszek
I LO w Tarnowie

Generator pseudolosowy Mersenne Twister

SPIS TREŚCI

Mersenne Twister

Opisane w poprzednim rozdziale liniowe generatory kongruencyjne liczb pseudolosowych posiadają wiele różnych wad eliminujących je z niektórych zastosowań. W celu usunięcia tych niedogodności, w roku 1997 dwaj japończycy, Makoto Matsumoto Takuji Nishimura, opracowali nowy rodzaj generatora liczb pseudolosowych o nazwie Mersenne Twister – w skrócie MT. Generator oparty został na liniowej rekurencji macierzowej (ang. matrix linear recursion)skończonej dziedzinie binarnej (ang. finite binary field). Umożliwia on generację wysokiej jakości liczb pseudolosowych o bardzo dużym okresie powtarzania, będącym liczbą pierwszą Mersenne'a.

Liczby Mersenne'a są potęgami liczby 2 pomniejszonymi o 1:

Mn = 2n-1

Oto kilka z nich:

1, 3, 7, 15, 31, 63, 127, 255, 511, 1023, 2047, 4095, 8191, 16383, 32767, 65535, …

Liczbą pierwszą Mersenne'a (ang. Mersenne Prime) nazywamy liczbę Mersenne'a, która jest pierwsza. Nie wszystkie liczby Mersenne'a są liczbami pierwszymi. Jedno z założeń mówi, iż  liczba Mersenne'a może być liczbą pierwszą tylko wtedy, gdy wykładnik n sam jest liczbą pierwszą. Jeśli wykładnik jest liczbą złożoną, to 2n-1 nie jest pierwsze. Pozwala to szybko eliminować z ciągu liczb Mersenne'a liczby złożone. Jednakże pozostałe liczby należy sprawdzać testem pierwszości, gdyż nie ma pewności, iż są one rzeczywiście pierwsze. Na przykład dla n = 11 otrzymujemy:

211-142047 = 23×89 – liczba złożona

Zadanie to jest trudne obliczeniowo, ponieważ liczby Mersenne'a bardzo szybko rosną. Na dzień dzisiejszy największą znaną liczbą Mersenne'a jest:

232582657-1

Jest to 44-ta liczba pierwsza Mersenne'a i posiada 9808358 cyfr dziesiętnych. Znaleziono ją 4 września 2006 roku dzięki projektowi GIMPS (ang. Great Internet Mersenne Prime Search), który polegał na  rozproszonym (tj. wykorzystującym moc obliczeniową komputerów podłączonych do sieci Internet) sprawdzaniu pierwszości liczb Mersenne'a.

Omawiany tutaj generator MT wykorzystuje 24-tą liczbę Mersenne'a o wartości 19937 – stąd bierze się jego oznaczenie MT19937. Generator posiada olbrzymi okres powtarzania równy 219937-1. Ponieważ teoria będąca podstawą działania generatora MT jest bardzo skomplikowana matematycznie (wykracza poza poziom liceum), ograniczymy się jedynie do podania samego algorytmu, bez wnikania w podstawy matematyczne jego działania. Wszystkie opisane tutaj fakty pochodzą z oryginalnej publikacji:

Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator

Makoto Matsumoto and Takuji Nishimura

Keio University/Max Planck Institut für Mathematik
Keio University

Cały dokument możesz pobrać stąd.

Generator MT jest jednym, dużym rejestrem przesuwnym ze sprzężeniami zwrotnymi powodującymi rotację oraz mieszanie się bitów. Rejestr ma długość 19937 bitów (liczba Mersenne), które w pamięci zajmują 624 słowa 32-bitowe. Generowanie liczb pseudolosowych składa się z trzech etapów.

  1. Pierwszy etap polega na zainicjowaniu 624 słów 32-bitowych rejestru wartościami pseudolosowymi. Do tego celu można wykorzystać zwykły liniowy generator LCG (my wykorzystujemy generator Super-Duper LCG(232, 69069, 0, X0)).
  2. W drugim etapie generujemy liczbę pseudolosową na podstawie zawartości rejestru. Rejestr jest odpowiednio modyfikowany przy tworzeniu każdej nowej liczby pseudolosowej.
  3. W trzecim etapie wygenerowaną liczbę odpowiednio dopasowujemy, aby otrzymać równomierny rozkład bitów.

Algorytm inicjowania rejestru MT

Wejście:

MT : tablica 624 liczb 32-bitowych przechowująca 19937 bitów rejestru przesuwnego.
X0 : ziarno dla generatora LCG(4294967296, 69069, 0), który posłuży do generacji zawartości rejestru przesuwnego. Wartość tę można pozyskiwać np. z zegara czasu rzeczywistego.

Wyjście:

Tablica MT wypełniona liczbami pseudolosowymi.

Zmienna pomocnicza

i : indeksuje kolejne liczby tablicy MT, i ∈ N.

Lista kroków:

K01: MT[0] ← X0 and $FFFFFFFF ; ziarno ustawiamy w pierwszym elemencie
K02: Dla i = 1, 2, …, 623
     wykonuj krok K03
K03:   MT[i] ← (69069×MT[i-1]) and $FFFFFFFF ; kolejne elementy wypełniamy
                                             ; liczbami pseudolosowymi
K04: Zakończ

Algorytm generatora pseudolosowego Mersenne Twister

Wejście:

MT  : tablica 624 liczb 32-bitowych przechowująca 19937 bitów rejestru przesuwnego.
mti : indeks, wg którego algorytm wybiera elementy MT do modyfikacji oraz do generacji wartości wyjściowej. Jest to zmienna statyczna, która zachowuje swoją wartość pomiędzy kolejnymi wywołaniami funkcji generującej liczby pseudolosowe. Powinna być zainicjowana wartością 0.

Wyjście:

Liczba pseudolosowa z zakresu od 0 do 232-1.

Elementy pomocnicze:

MA : dwuelementowa tablica, która służy do mnożenia elementów MT przez wybraną macierz.
    A[0] = 0
    A[1] = $9908B0DF
mti : indeks, wg którego algorytm wybiera elementy MT  do modyfikacji oraz do generacji wartości wyjściowej. Jest to zmienna statyczna, która zachowuje swoją wartość pomiędzy kolejnymi wywołaniami funkcji generującej liczby pseudolosowe. Powinna być zainicjowana wartością 0.
y : generowana liczba pseudolosowa.

Lista kroków:

K01: y ← (MT[mti] and $80000000) or (MT[(mti+1) mod 624] and $7FFFFFFF) ; łączymy dwa wyrazy z MT
K02: MT[mti] ← MT[(mti+397) mod 624] xor (y shr 1) xor MA[y and 1] ; tworzymy sprzężenie zwrotne
K03: y ← MT[mti] ; wyznaczamy liczbę pseudolosową
K04: y ← y xor (y shr 11)
K05: y ← y xor ((y shl 7) and $9D2C5680)
K07: y ← y xor ((y shl 15) and $EFC60000)
K08: y ← y xor (y shr 18)
K09: mti ← (mti+1) mod 624 ; obliczamy następny indeks
K10: Zakończ z wynikiem y

Przykładowe programy

Uwaga:

Zanim uruchomisz program, przeczytaj wstęp do tego artykułu, w którym wyjaśniamy funkcje tych programów oraz sposób korzystania z nich.

Program odczytuje z pierwszego wiersza trzy liczby a, b oraz n. Liczby ab określają przedział całkowity, w którym mają być wygenerowane liczby pseudolosowe. Liczba n  określa ile liczb pseudolosowych należy wygenerować. Liczby pseudolosowe są generowane w kolejnych wierszach.
Pascal
// Generator Mersenne Twister
// Data   : 16.04.2008
// (C)2020 mgr Jerzy Wałaszek
//---------------------------

program prg;

uses SysUtils, DateUtils;

var

  MT  : array [0..623] of longword;
  mti : integer;

// Inicjuje MT []
//--------------
procedure InicjujMT (x0 : longword);

var
  i : integer;
  x : qword;
begin
  MT[0] := x0;
  for i := 1 to 623 do
  begin
    x := MT[i-1];
    x := (23023*x) and $00000000ffffffff;
    x := (    3*x) and $00000000ffffffff;
    MT[i] := x;
  end;
  mti := 0;
end;

// Inicjuje MT wartościami losowymi
//---------------------------------
procedure UprzypadkowijMT;

var
  t  : TDateTime;
  x0 : longword;
begin
  t  := Now;
  x0 := HourOf(t);
  x0 := x0*  60+MinuteOf(t);
  x0 := x0*  60+SecondOf(t);
  x0 := x0*1000+MillisecondOf(t);
  InicjujMT(x0);
end;

// Generator Mersenne Twister
//--------------------------
function MersenneTwister : longword;
const
  MA : array [0..1] of longword = (0, $9908b0df);
var
  y : longword;
  i1, i397 : integer;
begin
  i1   := mti+  1; if   i1 > 623 then i1 := 0;
  i397 := mti+397; if i397 > 623 then dec(i397, 624);
  y    := (MT[mti] and $80000000) or (MT[i1] and $7fffffff);
  MT[mti] := MT[i397] xor (y shr 1) xor MA[y and 1];
  y    := MT[mti];
  y    := y xor (y shr 11);
  y    := y xor ((y shl  7) and $9d2c5680);
  y    := y xor ((y shl 15) and $efc60000);
  y    := y xor (y shr 18);
  mti  := i1;
  MersenneTwister := y;
end;

var
  a, b, i, n : longint;
begin
  UprzypadkowijMT;
  readln(a, b, n);
  for i := 1 to n do
    write(a+(MersenneTwister mod (b-a+1)), ' ');
  writeln;
end.
C++
// Generator Mersenne Twister
// Data   : 16.04.2008
// (C)2020 mgr Jerzy Wałaszek
//---------------------------

#include <iostream>

using namespace std;

typedef unsigned long long ulong;

unsigned int MT[624];
int mti = 0;

// Inicjuje MT[]
//-------------
void InicjujMT(unsigned int x0)
{
  ulong x;

  MT[0] = x0;
  for(int i = 1; i < 623; i++)
  {
    x = MT[i-1];
    x = (23023*x) & 0xffffffffull;
    x = (    3*x) & 0xffffffffull;
    MT[i] = x;
  }
}

// Inicjuje MT wartościami losowymi
//---------------------------------
void UprzypadkowijMT()
{
  InicjujMT((unsigned int)time(NULL));
}

// Generator Mersenne Twister
//--------------------------
unsigned int MersenneTwister()
{
  const unsigned int MA[] = {0, 0x9908b0df};
  long int y;
  int i1, i397;

  i1      = mti+  1; if(  i1 > 623) i1    = 0;
  i397    = mti+397; if(i397 > 623) i397 -= 624;
  y       = (MT[mti] & 0x80000000) | (MT[i1] & 0x7fffffff);
  MT[mti] = MT[i397] ^ (y >> 1) ^ MA[y & 1];
  y       = MT[mti];
  y      ^=  y >> 11;
  y      ^= (y <<  7) & 0x9d2c5680;
  y      ^= (y << 15) & 0xefc60000;
  y      ^=  y >> 18;
  mti     = i1;
  return y;
}

int main()
{
  int a, b, i, n;

  UprzypadkowijMT();
  cin >> a >> b >> n;
  for(i = 1; i <= n; i++)
    cout << (a+(MersenneTwister() % (b-a+1))) << " ";
  cout << endl << endl;

  return 0;
}
Basic
' Generator Mersenne Twister
' Data   : 16.04.2008
' (C)2020 mgr Jerzy Wałaszek
'---------------------------

Dim Shared MT(623) As Uinteger
Dim Shared mti As Integer = 0

' Inicjuje MT[]
'-------------
Sub InicjujMT(Byval x0 As Uinteger)

  Dim i As Integer, x As Ulongint

  MT (0) = x0
  For i = 1 To 623
    x = MT(i-1)
    x = (23023*x) And &Hffffffff
    x = (    3*x) And &Hffffffff
    MT(i) = x
  Next 
End Sub

' Inicjuje MT wartościami losowymi
'---------------------------------
Sub UprzypadkowijMT()
  InicjujMT(Timer*1000)
End Sub

' Generator Mersenne Twister
'--------------------------
Function MersenneTwister() As Uinteger

  Dim MA(1) As Uinteger => {0, &H9908b0df}
  Dim y As Uinteger
  Dim As Integer i1, i397

  i1      = mti+  1: If   i1 > 623 Then i1 = 0
  i397    = mti+397: If i397 > 623 Then i397 -= 624
  y       = (MT(mti) And &H80000000) Or (MT(i1) And &H7fffffff)
  MT(mti) = MT(i397) Xor (y Shr 1) Xor MA (y And 1)
  y       = MT (mti)
  y       = y Xor (y  Shr 11)
  y       = y Xor ((y Shl  7) And &H9d2c5680)
  y       = y Xor ((y Shl 15) And &Hefc60000)
  y       = y Xor (y  Shr 18)
  mti     = i1
  MersenneTwister = y
End Function

Dim As Integer a, b, i, n

UprzypadkowijMT()
Open Cons For Input As #1
Input #1, a, b, n
Close #1
For i = 1 To n
  Print a+(MersenneTwister() Mod (b-a+1));
Next
Print
End
Python (dodatek)
# Generator Mersenne Twister
# Data   : 16.02.2024
# (C)2024 mgr Jerzy Wałaszek
# ---------------------------

import time

# Zmienne globalne
mt = []
mti = 0


# Inicjuje MT[]
# -------------
def inicjuj_mt(x0):
    global mt
    mt.append(x0)
    for i in range(1, 624):
        x = mt[i-1]
        x = (23023*x) & 0xffffffff
        x = (    3*x) & 0xffffffff
        mt.append(x)
    return


# Inicjuje MT wartościami losowymi
# ---------------------------------
def uprzypadkowij_mt():
    inicjuj_mt(int(time.time()*1000))
    return


# Generator Mersenne Twister
# --------------------------
def mersenne_twister():
    global mti, mt
    ma = [0, 0x9908b0df]
    i1 = mti+1
    if i1 > 623: i1 = 0
    i397 = mti+397
    if i397 > 623: i397 -= 624
    y = (mt[mti] & 0x80000000) | (mt[i1] & 0x7fffffff)
    mt[mti] = mt[i397] ^ (y >> 1) ^ ma[y & 1]
    y = mt[mti]
    y ^= y >> 11
    y ^= (y << 7) & 0x9d2c5680
    y ^= (y << 15) & 0xefc60000
    y ^= y >> 18
    mti = i1
    return y

uprzypadkowij_mt()
arr = input().split()
a = int(arr[0])
b = int(arr[1])
n = int(arr[2])
for i in range(n):
    print(a+(mersenne_twister() % (b-a+1)), end=" ")
print()
Wynik:
1 6 200
6 2 3 6 2 5 1 5 4 3 2 3 2 4 1 1 3 6 4 2 4 5 5 4 4 5 3 4 4 3 3 6 1 3 5 1 1 4 4 3
4 3 2 5 1 2 2 1 6 3 1 1 1 3 4 6 6 4 1 5 1 3 4 2 5 1 6 3 2 5 5 5 5 6 5 6 3 2 5 3
5 5 4 6 5 1 5 4 4 6 5 6 1 5 3 1 6 3 4 5 3 6 5 3 5 6 5 5 6 2 2 1 1 4 6 1 5 3 4 5
2 1 1 3 6 1 1 2 4 5 4 5 4 4 1 1 1 1 4 2 3 3 1 1 5 1 2 3 2 1 4 4 6 5 3 6 6 4 2 3
6 5 1 6 4 1 5 2 2 6 2 1 6 2 6 2 4 3 4 1 6 3 2 6 2 3 3 4 4 1 6 3 2 5 5 3 3 4 2 6

Na początek:  podrozdziału   strony 

Zespół Przedmiotowy
Chemii-Fizyki-Informatyki

w I Liceum Ogólnokształcącym
im. Kazimierza Brodzińskiego
w Tarnowie
ul. Piłsudskiego 4
©2024 mgr Jerzy Wałaszek

Materiały tylko do użytku dydaktycznego. Ich kopiowanie i powielanie jest dozwolone
pod warunkiem podania źródła oraz niepobierania za to pieniędzy.

Pytania proszę przesyłać na adres email: i-lo@eduinf.waw.pl

Serwis wykorzystuje pliki cookies. Jeśli nie chcesz ich otrzymywać, zablokuj je w swojej przeglądarce.

Informacje dodatkowe.