Generator liczb pseudolosowych (cz.1)

No i stało się... teraz jakoś trzeba zacząć :) Przeglądając zakamarki starego dysku natknąłem się na pewen projekt. Kodzik powstał przy okazji zaliczania Metaheurystyk parę lat temu, a ponieważ wydał mi się ciekawy i dawno nie pisałem już w C++, postanowiłem wyciągnąć go na światło dzienne... Temat stary jak same komputery - generowanie liczb pseudolosowych. Nuda? :) Cóż, niekoniecznie... czasem przydaje się coś nieco bardziej wymyślnego (czyt. o lepszych właściwościach statystycznych) niż klasyczne:

#define frand() ((double) rand() / (RAND_MAX+1.0))

To co? Jedziemy ? ;-) Teorii sobie i Wam oszczędzę... dociekliwych i z zacięciem matematycznym odsyłam do tej książeczki [2]. W niej to opisane są poniższe algorytmy. Kod zamieszczam jedynie w celach poznaczych, a studentów WSISIZ przestrzegam przed metodą Kopiego-Paste'a, bo się skończy brakem zaliczenia ;-)

Chcemy uzyskać generator o rozkładzie normalnym oraz Cauchy'ego... Na początek jednak jako podstawa potrzebny nam będzie generator o rozkładzie równomiernym.

Zdefiniujmy sobie klasę RandomNumberGenerator:

class RandomNumberGenerator {
private:
   double a,b,c;
   static RandomNumberGenerator* instance;
protected:
   RandomNumberGenerator();
public:
   static RandomNumberGenerator* getInstance();
   void init();
   double getFromUniformDistribution();
   double getFromNormalDistribution();
   double getFromCauchyDistribution();
};

Klasa jak widać będzie miała trzy metody do pobierania losowych wartości z różnych rozkładów. Zaimplementujmy pierwszy rozkład, używając do tego celu tzw. ogólnego generatora liniowego. Generator ten realizuje schemat liniowy:

X_n+1 = (a₁X_n + a₂X_n-1 + ... + a_kX_n-k+1 + c) mod m

Okres tego generatora wynosi 2⁹⁶, generator spełnia też testy statystyczne DIEHARD, wydaje się więc całkiem OK. Inicjowany jest trzema nieujemnymi liczbami całkowitymi (zmienne prywatne a,b,c) - można do tego celu użyć funkcji time(), zapewniając różne wyniki za każdym uruchomieniem programu. Ustalenie tych wartości stałymi spowoduje oczywiście generowanie takiego samego ciągu liczb za każdym razem (algorytm jest deterministyczny). Ponieważ znane są metody określania parametrów początkowych na podstawie obserwacji generowanego ciągu, nie jest to algorytm bezpieczny do zastosowań kryptograficznych[1].

Przedstawiona w książce Wieczorkowskiego i Zielińskiego[2] implementacja generuje liczby z przedziału (-0.5, 0.5), dostosowałem ją więc tak aby przedział ten wynosił (0,1) - taki, jaki jest potrzebny do późniejszego użycia w generatorze o rozkładzie normalnym. Inicjacja generatora następuje w metodzie init(), a sam algorytm w metodzie getFromUniformDistribution():

#include <math.h>
#include <values.h>
#include <time.h>
#include "RandomNumberGenerator.h"

using namespace std;

RandomNumberGenerator* RandomNumberGenerator::instance = 0;

RandomNumberGenerator* RandomNumberGenerator::Instance(){
if (instance == 0){
       instance = new RandomNumberGenerator;
       instance->init();
}
return instance;
}

RandomNumberGenerator::RandomNumberGenerator(){}

void RandomNumberGenerator::init() {
a = time(NULL);
b = time(NULL);
c = time(NULL);
}

double RandomNumberGenerator::getFromUniformDistribution(){

static int n;
static double d,x;

d = (a + b + c) * 8192;
x = fmod (d,4294967291.0);
a = b; b = c; c= x;

if (x < (float) 2147483648.0) {
       n = (int) x;
} else {
       n = (int) (x - 4294967296.0);
}

return (n * 2.3283064365e-10 + 0.5);
}

CDN...

[1] http://pl.wikipedia.org/wiki/Generator_liczb_pseudolosowych
[2] "Komputerowe generatory liczb losowych" R. Wieczorkowski R Zieliński, WNT 1997