Mersenne Twister
Author |
Message |
md
|
Posted: Mon Apr 24, 2006 10:28 pm Post subject: Mersenne Twister |
|
|
I've grown fedup with the very poor random number generators in C/C++; so I wrote this simple routine to get much better random numbers. Implements the Mersenne Twister
c++: |
/*
The Mersenne Twister
The Mersenne Twister is a random number generator, invented/discovered in 1996 by Matsumora and Nishimura.
MT is a twisted GFSR(624,397), similar in spirit to R250 and in speed to R250. MT has an period of 2^19937-1.
It is a very good random number generator.
*/
#include <cstdlib>
#include <ctime>
#include <fstream>
const int MT_LEN = 624;
int mt_index;
unsigned long mt_buffer[MT_LEN];
void mt_init()
{
#ifdef __linux__
std::ifstream entropy("/dev/urandom", std::ios::in | std::ios::binary);
int i;
srand(time(0));
for (i = 0; i < MT_LEN; i++)
entropy.read((char*)(& mt_buffer[i]), sizeof(unsigned long));
mt_index = 0;
entropy.close();
#else
// default seeding
int i;
srand(time(0));
for (i = 0; i < MT_LEN; i++)
mt_buffer[i] = rand();
mt_index = 0;
#endif
}
const int MT_IA = 397;
const int MT_IB = MT_LEN - MT_IA;
const unsigned long UPPER_MASK = 0x80000000;
const unsigned long LOWER_MASK = 0x7FFFFFFF;
const unsigned long MATRIX_A = 0x9908B0DF;
unsigned long mt_random()
{
unsigned long * b = mt_buffer;
int idx = mt_index;
unsigned long s;
int i;
if (idx == MT_LEN*sizeof(unsigned long))
{
idx = 0;
i = 0;
for (; i < MT_IB; i++)
{
s = ((b)[i] & UPPER_MASK) | ((b)[i+1] & LOWER_MASK);
b[i] = b[i + MT_IA] ^ (s >> 1) ^ (((s)&1)*MATRIX_A);
}
for (; i < MT_LEN-1; i++)
{
s = ((b)[i] & UPPER_MASK) | ((b)[i+1] & LOWER_MASK);
b[i] = b[i - MT_IB] ^ (s >> 1) ^ (((s)&1)*MATRIX_A);
}
s = ((b)[MT_LEN-1] & UPPER_MASK) | ((b)[0] & LOWER_MASK);
b[MT_LEN-1] = b[MT_IA-1] ^ (s >> 1) ^ (((s)&1)*MATRIX_A);
}
mt_index = idx + sizeof(unsigned long);
return *(unsigned long *)((unsigned char *)b + idx);
}
|
Code is also on svn://svn.nxor.org/public/mt_rng |
|
|
|
|
|
Sponsor Sponsor
|
|
|
wtd
|
Posted: Mon Apr 24, 2006 11:59 pm Post subject: (No subject) |
|
|
What about refactoring this such that a random number generator is an object that automatically inits when its constructor is called? |
|
|
|
|
|
md
|
Posted: Tue Apr 25, 2006 9:26 am Post subject: (No subject) |
|
|
Good idea! That'll be for tonight |
|
|
|
|
|
wtd
|
Posted: Tue Apr 25, 2006 10:17 am Post subject: (No subject) |
|
|
code: | class mt_rang
{
// ...
template <typename t>
mt_rng& operator>>(t& dest)
{
// get next random number, cast it to t,
// and store in dest
return *this;
}
}; |
Then:
code: | generator >> i >> j; |
|
|
|
|
|
|
md
|
Posted: Tue Apr 25, 2006 7:27 pm Post subject: (No subject) |
|
|
New code just for you wtd
cMTRandom.h
c++: |
/*
The Mersenne Twister
The Mersenne Twister is a random number generator, invented/discovered in 1996 by Matsumora and Nishimura.
MT is a twisted GFSR(624,397), similar in spirit to R250 and in speed to R250. MT has an period of 2^19937-1.
It is a very good random number generator.
*/
#ifndef cmtrandom_include
#define cmtrandom_include
class cMTRandom
{
public:
cMTRandom(void); // no seed
cMTRandom(unsigned long seed); // seeded (but why?!)
unsigned long Rand(void);
template <typename T> inline T& operator=(cMTRandom& boogie) { return static_cast<T> (boogie.Rand()); }
template <typename T> inline cMTRandom& operator>>(T& boogie) { boogie = static_cast<T>(Rand()); return *this; }
protected:
int index;
unsigned long MT[623]; //623 is a semi-magic number; it's part of the algorithm and
// not chosen at random.
void GenerateNumbers(void);
};
#endif
|
cMTRandom.cpp
c++: |
/*
The Mersenne Twister
The Mersenne Twister is a random number generator, invented/discovered in 1996 by Matsumora and Nishimura.
MT is a twisted GFSR(624,397), similar in spirit to R250 and in speed to R250. MT has an period of 2^19937-1.
It is a very good random number generator.
*/
#include <cstdlib>
#include <ctime>
#include <fstream>
#include "cMTRandom.h"
/*
Constants
*/
const unsigned long UPPER_MASK = 0x80000000;
const unsigned long LOWER_MASK = 0x7FFFFFFF;
const unsigned long MATRIX = 0x9908B0DF;
const unsigned long MANGLE_A = 0x9D2C5680;
const unsigned long MANGLE_B = 0xEFC60000;
/*
Variables
*/
cMTRandom::cMTRandom(void)
{
/*
On linux we can use /dev/urandom to get our starting seed, as it's
less predictable then time() or srand()/rand() seeded from a specific
time.
On other systems just use time() since it's just as random as srand()/
rand() seeded at time(), and a quicker.
*/
#ifdef __linux__
std::ifstream entropy("/dev/urandom", std::ios::in | std::ios::binary);
entropy.read((char*)(&MT[0]), sizeof(MT[0]));
entropy.close();
#else
MT[0] = time(0);
#endif
// now initialize the generator
for (int i = 1; i < 624; i++)
MT[i] = (69069 * MT[i-1]);
// generate the initial set of numbers
GenerateNumbers();
index = 0;
}
cMTRandom::cMTRandom(unsigned long seed)
{
// now initialize the generator
MT[0] = seed;
for (int i = 1; i < 624; i++)
MT[i] = (69069 * MT[i-1]);
// generate the initial set of numbers
GenerateNumbers();
index = 0;
}
void cMTRandom::GenerateNumbers(void)
{
int y;
for (int i = 0; i < 623; i++)
{
y = (MT[i] & UPPER_MASK) + (MT[i+1] & LOWER_MASK);
if (y & 1 == 0)
MT[i] = MT[(i + 397) % 624] ^ (y >> 1);
else
MT[i] = MT[(i + 397) % 624] ^ (y >> 1) ^ MATRIX;
}
y = (MT[623] & UPPER_MASK) + (MT[0] & LOWER_MASK);
if ( y & 1 == 0)
MT[623] = MT[396] ^ (y >> 1);
else
MT[623] = MT[396] ^ (y >> 1) ^ MATRIX;
}
unsigned long cMTRandom::Rand(void)
{
if (index == 624)
{
index = 0;
GenerateNumbers();
}
unsigned long x = ((MT[index] & 0xFFFF0000) & MANGLE_A) ^ (MT[index] << 16);
unsigned long y = ((MT[index] & 0x0000FFFF) & MANGLE_B) ^ (MT[index] >> 16);
index++;
return (x ^ y);
}
|
|
|
|
|
|
|
|
|