(root)/src/tools/randomgenerator.h - Rev 706
Rev 705 |
Rev 707 |
Go to most recent revision |
Blame |
Compare with Previous |
Last modification |
View Log
| RSS feed
#ifndef RANDOMWELL_H
#define RANDOMWELL_H
#include <cstdlib>
#include <math.h>
#include <time.h>
// see http://www.lomont.org/Math/Papers/2008/Lomont_PRNG_2008.pdf
// for the
class MTRand
{
public:
MTRand() {seed();}
/// get a random value from [0., 1.]
inline double rand() { return double(random_function()) * (1.0/4294967295.0); }
/// get a random value from [0., max_value]
inline double rand(const double max_value) { return max_value * rand(); }
/// get a random integer in [0,2^32-1]
inline unsigned long randInt(){ return random_function(); }
inline unsigned long randInt(const int max_value) { return max_value>0?randInt()%max_value:0; }
/// Access to nonuniform random number distributions
inline double randNorm( const double mean, const double stddev );
void seed();
void seed(unsigned int oneSeed);
private:
inline unsigned int random_function() { /* return WELLRNG512(); */
/* return xorshf96(); */
return fastrand(); }
/* initialize state to random bits */
unsigned long state[16];
/* init should also reset this to 0 */
unsigned int index;
/* return 32 bit random number */
inline unsigned long WELLRNG512(void)
{
unsigned long a, b, c, d;
a = state[index];
c = state[(index+13)&15];
b = a^c^(a<<16)^(c<<15);
c = state[(index+9)&15];
c ^= (c>>11);
a = state[index] = b^c;
d = a^((a<<5)&0xDA442D24UL);
index = (index + 15)&15;
a = state[index];
state[index] = a^b^d^(a<<2)^(b<<18)^(c<<28);
return state[index];
}
// The Marsaglia's xorshf generator:
// see: http://stackoverflow.com/questions/1640258/need-a-fast-random-generator-for-c and
// http://www.cse.yorku.ca/~oz/marsaglia-rng.html
//unsigned long x=123456789, y=362436069, z=521288629;
unsigned long x,y,z;
unsigned long xorshf96(void) { //period 2^96-1
unsigned long t;
x ^= x << 16;
x ^= x >> 5;
x ^= x << 1;
t = x;
x = y;
y = z;
z = t ^ x ^ y;
return z;
}
int g_seed;
inline unsigned int fastrand()
{
g_seed = (214013*g_seed+2531011);
return g_seed;
}
};
inline void MTRand::seed()
{
srand ( time(NULL) );
for (int i=0;i<16;i++)
state[i] = std::rand();
index = 0;
// inits for the xorshift algorithm...
x=123456789, y=362436069, z=521288629;
// inits for the fast rand....
g_seed = std::rand();
}
inline void MTRand::seed(unsigned int oneSeed)
{
srand ( oneSeed );
for (int i=0;i<16;i++)
state[i] = std::rand();
index = 0;
x=123456789, y=362436069, z=521288629;
g_seed = oneSeed;
}
inline double MTRand::randNorm( const double mean, const double stddev )
{
// Return a real number from a normal (Gaussian) distribution with given
// mean and standard deviation by polar form of Box-Muller transformation
double x, y, r;
do
{
x = 2.0 * rand() - 1.0;
y = 2.0 * rand() - 1.0;
r = x * x + y * y;
}
while ( r >= 1.0 || r == 0.0 );
double s = sqrt( -2.0 * log(r) / r );
return mean + x * s * stddev;
}
#endif // RANDOMWELL_H