Porówanie wydajności Mersenne Twister z innymi generatorami liczb pseudolosowych:
Jak widać w powyższej tabeli algorytm MT19937 ma znaczną przewagę nad resztą, jeśli chodzi o stosunek okresu losowości do czasu wygenerowania \( 10^7 \) liczb pseudolosowych w danej przestrzeni roboczej.
Warto zaznaczyć, że algorytm Mersenne Twister nie nadaje się do użycia w kryptografii. Liczby przez niego generowane da się łatwo zgadnąć transformując liniowo macierz "hartującą" używaną przez algorytm, na jej odwrotność. Wtedy to używając liniowej rekurencji możemy zgadnąć dalsze ciągi generowanych liczb poprzez analizę wygenerowanych liczb.
Algorytm Mersenne Twister jest mniej więcej tak samo szybki jak algorytm rand znany z ANSI-C, co więcej zawiera o wiele większy okres i stopień równomiernego rozmieszczenia, co daje nam pseudolosowość o lepszej jakości. Poniżej zaprezentuje kod Mersenne Twister'a w C++, który pozwala zarówno na generowanie pseudolosowych liczb całkowitych, jak i też zmiennoprzecinkowych.
#include <iostream> #include <ctime> using namespace std; // Parametry okresu const unsigned long N = 624; const unsigned long M = 397; const unsigned long MATRIX_A = 0x9908b0df; const unsigned long UPPER_MASK = 0x80000000; const unsigned long LOWER_MASK = 0x7fffffff; // Parametry hartowania macierzy const unsigned long TEMPERING_MASK_B = 0x9d2c5680; const unsigned long TEMPERING_MASK_C = 0xefc60000; // Makra przesuwajace bity #define TEMPERING_SHIFT_U(y) (y >> 11) #define TEMPERING_SHIFT_S(y) (y << 7) #define TEMPERING_SHIFT_T(y) (y << 15) #define TEMPERING_SHIFT_L(y) (y >> 18) static unsigned long mt[N]; // tablica wektora stanu static unsigned int mti = N+1; // Funkcja inicjalizujaca tablice sola niezerowa void sgenrand(unsigned long seed) { mt[0] = seed & 0xffffffff; for (mti=1; mti<N; mti++) { mt[mti] = (69069 * mt[mti-1]) & 0xffffffff; } } // Glowna funkcja losujaca zwracajaca typ calkowity unsigned long genrand() { unsigned long y; static unsigned long mag01[2] = {0x0, MATRIX_A}; if (mti >= N) { unsigned int kk; // Jesli sgenrand nie zostalo wywolane inicjalizuje tablice domyslnie if (mti == N+1) sgenrand(4357); for (kk=0; kk<N-M; kk++) { y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK); mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1]; } for (; kk<N-1; kk++) { y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK); mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1]; } y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK); mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1]; mti = 0; } y = mt[mti++]; y ^= TEMPERING_SHIFT_U(y); y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B; y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C; y ^= TEMPERING_SHIFT_L(y); return y; } // Funkcja losujaca zwracajca typ zmiennoprzecinkowy double genrand_d() { return (static_cast<double>(genrand()) / static_cast<unsigned long>(0xffffffff)); } // Program wypisujacy pierwsze 100 wylosowanych liczb z przedzialu 0-100 int main() { int min = 0; int max = 100; // Inicjalizacja generatora aktualnym czasem sgenrand(time(NULL)); for (int i=0; i<100; i++) { cout << (genrand()%(max-min))+min << "\t"; if (i%10==9) cout << "\n"; } return 0; }
Działanie algorytmu polega na początkowym utworzeniu 624 wektorów słów o długości 32 (oznaczmy ją w) lub 64 bitów i wartościach od 0 do \(2^{32}-1\) lub \(2^{64}-1\). To na tych słowach operuje Mersenne Twister w swych późniejszych fazach wykonania. Algorytm bazuje na rekurencji liniowej określonej wzorem:
$$ x_{k+n}=x_{k+m} ⨁ {(x_{k}^{u}│x_{k+1}^{l})A} $$
k to liczby całkowite postaci 0, 1, 2, 3, ... x oznacza wektor składający się z wyrazów 0 lub 1 i reprezentujący słowo. Liczba całkowita n to stopień rekurencji, liczba r (\(0≤r≤w-1\)), ukryta w implementacji \(x_{k}^{u}\) dzieli wektor na pół. Liczba m taka że \(1≤m≤n\). Określona jest także stała macierz wielkości \(w×w\) zapisana jako A składająca się ze zbioru zer i jedynek. Każdy z wektorów od zerowego do 623-go inicjalizujemy wartością początkową. Następnie generowane są wartości \(x_{k+n}\) poprzez podstawianie k w postaci 0, 1, 2, ...; Zapis \((x_{k}^{u}│x_{k+1}^{l})\) interpretujemy jako połączenie górnych bitów o długości w - r należących do wektora \(x_{k}\) i dolnych bitów o długości r należących do wektora \(x_{k+1}\). Tak połączony wektor mnożymy przez macierz A, zaś iloczyn dodajemy modulo dwa (XOR) do wektora \(x_{k+m}\). Tak otrzymany wynik generuje wektor \(x_{k+n}\). Macierz A ma postać: $$ A=\begin{pmatrix} &1&&&\\&&1&&\\&&&⋱&\\&&&&1\\a_{w-1}&a_{w-2}&…&…&a_{0} \end{pmatrix} $$ Pozwala to na proste mnożenie macierzy przebiegające tak: $$ xA=\begin{cases} \text{przesun_bity_w_prawo(x)} & x_0=0 \\ \text{przesun_bity_w_prawo(x)} ⨁ a & x_0=1 \end{cases} $$ , gdzie \( a=(a_{w-1}, a_{w-2}, ..., a_0) \), \( x=(x_{w-1}, x_{w-2}, ..., x_0) \) Rekurencję obliczamy z użyciem operatorów bitowych XOR, OR i AND.
Aby ulepszyć równomierne rozmieszczenie algorytmu stosuje się macierz odwrotną T o rozmiarze \(w×w\). Aby użyć jej w algorytmie należy zastosować poniższe operacje: $$ y=x⨁(x>>u) \\ y=y⨁(y>>s) AND b \\ y=y⨁(y>>t) AND c \\ y=y⨁(y>>l) $$, gdzie u, s, t, l to liczby całkowite, natomiast b i c to maski bitowe. Aby rekurencja działała należy stworzyć tablicę liczb całkowitych o długości n. Zrozumienie algorytmu wymaga najpierw zrozumienia jego strony matematycznej, szerzej opisanej w dokumencie przygotowanym przez samych autorów Mersenne Twister'a. Postarałem się ją nieco przybliżyć, lecz osoby, które chcą wejść głębiej w strukturę algorytmu muszą zajrzeć do źródła.
Źródło: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
Brak komentarzy:
Prześlij komentarz