Algorytm Mersenne Twister - generator liczb pseudolosowych

Z pewnością zalicza się do jednego z ciekawszych i lepszych generatorów liczb pseudolosowych. Został stworzony i opracowany przez Makoto Matsumoto i Takuji Nishimura z Uniwersytetu Keio. Okres algorytmu wynosi \( 2^{19937}-1 \), natomiast stopień równomiernego rozmieszczenia jest 623-wymiarowy. Algorytm stosuje dwa techniczne zabiegi polepszające jego wydajność i zdolność do generowania pseudolosowości. Pierwszą z nich są tzw. niekompletne tablice służące do znajdowania okresu liczby pierwszej Mersenne'a. Drugą z nich jest szybki algorytm testujący zdolność do rozkładu wielomianu poprzez rekurencję liniową. Do życia algorytm potrzebuje tylko zdefiniowania rekurencji i szybkiego algorytmu obliczania wektora stanu z 1-bitowego strumienia wyjścia. Priorytetem Mersenne Twister'a jest złożoność obliczeniowa jego algorytmu odwrotnego rozkładu wielomianu pomnożona przez stopień danego wielomianu. W celu osiągnięcia bardziej równomiernego rozmieszczenia stosuje się m.in. algorytm Lenstry.

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

Algorytm Euklidesa - obliczanie NWD

Jest to jeden z najstarszych algorytmów służących do obliczania największego wspólnego dzielnika (NWD) i zarazem jeden z najprostszych algorytmów, jakie tylko sobie można wymarzyć. Został wymyślony przez Euklidesa, słynnego greckiego matematyka.
Danymi wejściowymi są dwie liczby naturalne, dla których szukamy NWD. Najpierw obliczamy resztę z dzielenia dwóch wybranych liczb. Następnie wartość drugiej z tych liczb przypisujemy, jako wartość pierwszej, a obliczoną wcześniej resztę z dzielenia zapisujemy jako drugą z tych liczb. Otrzymujemy dwie nowe liczby, dla których powtarzamy powyższy cykl operacji, dopóki druga z liczb nie będzie równa 0. Na koniec działania algorytmu największym wspólnym dzielnikiem jest pierwsza z tych liczb.
Aby można było łatwiej to sobie wyobrazić zapiszemy dwie liczby wejściowe jako a oraz b, gdzie a,b należą do zbioru liczb naturalnych.

• Oblicz resztę z dzielenia liczb a i b.
• Przypisz liczbie a liczbę b.
• Przypisz liczbie b resztę z dzielenia liczb a i b obliczoną w pierwszym kroku.
• Jeśli b jest różne od zera, przejdź do kroku 1.

Poniżej prezentuję program w języku C++ implementujący algorytm Euklidesa wykorzystujący operator modulo (reszty z dzielenia):

#include <iostream>

using namespace std;

int main()
{
    int a, b, c;
    cout << "Wprowadz dwie liczby: ";
    cin >> a;
    cin >> b;
    while (b != 0) {
        c = a%b;
        a = b;
        b = c;
    }

    cout << "NWD tych liczb to: " << a << endl;

    return 0;
}

Jest to pierwsza możliwa implementacja algorytmu Euklidesa stosującego resztę z dzielenia. Istnieje jeszcze jedna możliwa stosująca odejmowanie, ale o niej za chwilę. Zajmiemy się omówieniem programu. Najpierw tworzymy 3 zmienne a, b, c reprezentujące liczby całkowite. Następnie wartości liczb a oraz b pobieramy ze standardowego wejścia za pomocą obiektu cin. Wykonujemy pętlę while dopóki liczba b nie będzie równa 0. Wewnątrz pętli do zmiennej c przypisujemy resztę z dzielenia liczb a i b. Potem do liczby a przypisujemy liczbę b, a do liczby b przypisujemy wartość zmiennej c, czyli naszą resztę z dzielenia. Jeśli b równe jest zeru wychodzimy z pętli, w przeciwnym wypadku wracamy na jej początek. Po wyjściu z pętli pozostaje nam tylko wypisanie liczby a będącej największym wspólnym dzielnikiem na standardowe wyjście przy pomocy obiektu cout.

Druga implementacja tego algorytmu korzysta nie z operatora modulo i reszty z dzielenia, lecz z odejmowania liczb od siebie. Jest nieco wolniejsza, gdyż trzeba wykonać sporą liczbę odejmowań dla liczb, które się sporo od siebie różnią, jednak warto o niej wspomnieć. Danymi wejściowymi tutaj również są dwie liczby naturalne a, b. Algorytm polega na tym, że dopóki któraś z liczb nie osiągnie zera odejmujemy większą od mniejszej. Po czym większą z liczb zastępujemy różnicą tych liczb i tak do zera.

• Oblicz różnicę liczb a i b.
• Przypisz większej liczbie różnicę obliczoną w kroku pierwszym.
• Jeśli a i b są różne od zera, przejdź do kroku 1.

Program w języku C++ implementujący algorytm Euklidesa stosujący odejmowanie:

#include <iostream>

using namespace std;

int main()
{
    int a, b, c;
    cout << "Wprowadz dwie liczby: " << endl;
    cin >> a;
    cin >> b;

    while (a != 0 && b != 0) {
        if (a > b) {
            c = a - b;
            a = c;
        } else {
            c = b - a;
            b = c;
        }
    }

    cout << "NWD: " << a + b << endl;

    return 0;
}

Tak jak w poprzedniej implementacji i tutaj będą nam potrzebne trzy zmienne a, b i c będące typem całkowitym. Wprowadzamy je ze standardowego wejścia. Następnie deklarujemy pętlę while, która wykonuje się dopóki a i b są różne od zera. Wewnątrz pętli stosujemy instrukcję warunkową. W przypadku, gdy a jest większe od b, przypisujemy do zmiennej c różnicę liczb a i b. Następnie do zmiennej a przypisujemy zmienną c, czyli otrzymaną różnicę. W innych przypadkach tj. gdy a jest mniejsze lub równe b, przypisujemy do zmiennej c różnicę liczb b i a, po czym inicjujemy zmienną b za pomocą wartości tej różnicy. Gdy pętla zakończy swe działanie, któraś ze zmiennych (a lub b) przechowuje największy wspólny dzielnik wobec czego najprostszą metodą, aby go odczytać nie stosując instrukcji warunkowych jest dodanie obu wartości, gdyż jedna z nich na pewno jest równa zeru. Sumę obu wartości wypisujemy na standardowe wyjście, co jest równoznaczne z wypisaniem NWD.