Mersenne Twister

Tässä:

MersenneTwisterInt32.h
------------------------------------
#ifndef MERSENNETWISTERINT32_H
#define MERSENNETWISTERINT32_H

#pragma once

class MersenneTwisterInt32
{
public:
MersenneTwisterInt32();
MersenneTwisterInt32(unsigned long seed);
MersenneTwisterInt32(const unsigned long* pArray, int size);
virtual ~MersenneTwisterInt32();

void Seed(unsigned long seed);
void Seed(const unsigned long* pArray, int size);

unsigned long operator()();

protected:
unsigned long GetRandomInt32();

private:
unsigned long Twiddle(unsigned long u, unsigned long v);
void GenerateState();

private:
static const int n = 624;
static const int m = 397;
static unsigned long state[n];
static int p;
static bool init;
};

inline unsigned long MersenneTwisterInt32::GetRandomInt32()
{
if (p == n)
{
GenerateState();
}

unsigned long x = state[p++];

x ^= (x >> 11);
x ^= (x 18);
}

inline unsigned long MersenneTwisterInt32::Twiddle(unsigned long u, unsigned long v)
{
return (((u & 0x80000000UL) | (v & 0x7FFFFFFFUL)) >> 1) ^ ((v & 1UL) ? 0x9908B0DFUL : 0x0UL);
}

#endif
----------------------------------------

MersenneTwisterInt32.cpp
----------------------------------------
#include "MersenneTwisterInt32.h"

unsigned long MersenneTwisterInt32::state[n] = {0x0UL};
int MersenneTwisterInt32::p = 0;
bool MersenneTwisterInt32::init = false;

MersenneTwisterInt32::MersenneTwisterInt32()
{
if (!init)
{
Seed(5489UL);
}

init = true;
}

MersenneTwisterInt32::MersenneTwisterInt32(unsigned long s)
{
Seed(s);

init = true;
}

MersenneTwisterInt32::MersenneTwisterInt32(const unsigned long* pArray, int size)
{
Seed(pArray, size);
}

MersenneTwisterInt32::~MersenneTwisterInt32()
{
}

void MersenneTwisterInt32::Seed(unsigned long s)
{
state[0] = s & 0xFFFFFFFFUL;

for (int i = 1; i < n; ++i)
{
state[i] = 1812433253UL * (state[i - 1] ^ (state[i - 1] >> 30)) + i;
state[i] &= 0xFFFFFFFFUL;
}

p = n;
}

void MersenneTwisterInt32::Seed(const unsigned long* pArray, int size)
{
int i = 1, j = 0;

Seed(19650218UL);

for (int k = ((n > size) ? n : size); k; --k)
{
state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1664525UL)) + pArray[j] + j;
state[i] &= 0xFFFFFFFFUL;

++j;

j %= size;

if ((++i) == n)
{
state[0] = state[n - 1];
i = 1;
}
}

for (int k = n - 1; k; --k)
{
state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1566083941UL)) - i;
state[i] &= 0xFFFFFFFFUL;

if ((++i) == n)
{
state[0] = state[n - 1];
i = 1;
}
}

state[0] = 0x80000000UL;
p = n;
}

unsigned long MersenneTwisterInt32::operator()()
{
return GetRandomInt32();
}

void MersenneTwisterInt32::GenerateState()
{
for (int i = 0; i < (n - m); ++i)
{
state[i] = state[i + m] ^ Twiddle(state[i], state[i + 1]);
}

for (int i = n - m; i < (n - 1); ++i)
{
state[i] = state[i + m - n] ^ Twiddle(state[i], state[i + 1]);
}

state[n - 1] = state[m - 1] ^ Twiddle(state[n - 1], state[0]);

p = 0;
}
-----------------------------------

Ja jos haluat esim. gaussiaaneja:
-----------------------------------
#ifndef MERSENNETWISTERGAUSSIAN_H
#define MERSENNETWISTERGAUSSIAN_H

#pragma once

#include "MersenneTwisterInt32.h"
#include

class MersenneTwisterGaussian: public MersenneTwisterInt32
{
public:
MersenneTwisterGaussian(): MersenneTwisterInt32() {};
MersenneTwisterGaussian(unsigned long seed): MersenneTwisterInt32(seed) {};
MersenneTwisterGaussian(const unsigned long* pArray, int size): MersenneTwisterInt32(pArray, size) {};
virtual ~MersenneTwisterGaussian() {};

double operator()();
};

inline double MersenneTwisterGaussian::operator()()
{
return std::sqrt(-2. * std::log(1. - static_cast(GetRandomInt32())
* (1. / 4294967296.))) * std::cos(static_cast(GetRandomInt32()) *
(2. * 3.14159265358979323846 / 4294967296.));
}


#endif