When I was in grad school I needed to do some heavy-duty Monte Carlo calculations. These require high-quality RNGs. Not cryptographic-grade, but each data point ended using more than 10[sup]8[/sup] random numbers. Obviously, system RNGs with periods of only 2[sup]32[/sup]~10[sup]9[/sup] would not be appropriate.
I did some literature research, and ended up with a paper by Marsaglia and Zaman (Comp. in Phys. 8, 117 (1994), sorry couldn’t find an online version). I coded it, and it worked well. It combines a linear congruential and a Fibonacci sequence for an overall period of around 2[sup]192[/sup]~10[sup]57[/sup]. BTW, all the better RNGs carefully combine two or more classes of simpler generators. This avoids the weaknesses of any single generator.
Here is a copy of the core of the C++ code. It’s written by me and is public domain, so feel free to use or abuse it (a version is published in my dissertation, but you’d never find it). The kernel is from Marsaglia and Zaman (recoded into C++). Each object of Random_t defines a separate instance of the RNG. The state can be initialized in its entirety, with a single integer, or with bytes pulled from /dev/urandom. An object can generate a uniformly distributed integer (0 to 2[sup]32[/sup]-1), a uniformly distributed double (0 <= X < 1.0), or a normally distributed double (mean = 0.0, variance = 1.0).
Class header file Random.h:
// -- Random.h
#ifndef __Random_h__
#define __Random_h__
#include <vector>
class Random_t
{
public:
typedef unsigned int Int_t; // MUST be 32 bits
typedef std::vector<Int_t> State_t;
private:
Int_t
n, // multiplicative modulo seed
v, w, x, y, z, // Fibonacci modulo seeds
c; // Fibonacci carry bit
public:
Random_t(void); // initialize from /dev/urandom
Random_t(Int_t S); // initialize from one integer
Random_t(const State_t& S); // initialize full state
Int_t operator()(void) {return next();} // random integer
double Uniform(void); // uniform distribution, 0.0 <= X < 1.0
double Normal(void); // normal distribution, mean = 0.0, variance = 1.0
State_t State(void) const; // returns current state vector
private:
Int_t next(void); // internal state updater
};
#endif // __Random_h__
// -- Random.h
Class code file Random.cpp:
// -- Random.cpp
#include <cmath>
#include <fstream>
#include "Random.h"
namespace // private namespace
{
template<typename Type_t>
void ReadBinary(std::istream& S, Type_t& X)
{
S.read((char*)(&X),sizeof(Type_t));
}
} // private namespace
Random_t::Random_t(void)
{
std::ifstream fin("/dev/urandom");
ReadBinary(fin,n);
ReadBinary(fin,v); ReadBinary(fin,w); ReadBinary(fin,x);
ReadBinary(fin,y); ReadBinary(fin,z); ReadBinary(fin,c);
}
Random_t::Random_t(Int_t S): n(S), v(n+2), w(v+3), x(w+5), y(x+7), z(y+11),
c(z+13)
{
for(size_t i=0; i<16; ++i) z=next(); // mix things for a better state
}
Random_t::Random_t(const State_t& S)
{
State_t s=S;
while(s.size()<7) s.push_back(s.size()); // fill out if state is too small
n=s[0]; v=s[1]; w=s[2]; x=s[3]; y=s[4]; z=s[5]; c=s[6];
}
double Random_t::Uniform(void)
{
static const double base=1.0/(Int_t(-1)+1.0);
return (next()*base+next())*base; // we need more than 32 bits
}
double Random_t::Normal(void)
{
static bool isCache=false;
static double cache=0.0;
if(isCache)
{
isCache=false;
return cache;
}
else
{
double v1, v2, vv;
do
{
v1=Uniform()*2.0-1.0;
v2=Uniform()*2.0-1.0;
vv=v1*v1+v2*v2;
}
while(vv>=1.0 || vv==0.0);
vv=sqrt(-2.0*log(vv)/vv);
isCache=true; cache=v1*vv;
return v2*vv;
}
}
Random_t::State_t Random_t::State(void) const
{
State_t s(7);
s[0]=n; s[1]=v; s[2]=w; s[3]=x; s[4]=y; s[5]=z; s[6]=c;
return s;
}
Random_t::Int_t Random_t::next(void)
// from G. Marsaglia, A. Zaman, Comp. in Phys. 8, 117 (1994)
// combines two independent pseudo-random number generators
// for a net period of about 2^192 (about 10^58)
// DO NOT TINKER with this function!!! (it breaks easily)
{
Int_t s=v+c; c=(y<s); s=y-s; if(c) s-=10;
v=w; w=x; x=y; y=z; z=s;
n=n*69069+1013904243;
return s+n;
}
// -- Random.cpp
And some test code TestRandom.cpp:
// -- TestRandom.cpp
#include <cstdlib>
#include <iostream>
#include "Random.h"
int main(int ArgNum, const char* Args[])
{
std::size_t num=10;
if(ArgNum>1) num=std::atol(Args[1]);
Random_t rnd1, rnd2, rnd3;
std::cout.precision(21);
for(; num; --num)
std::cout<<rnd1()<<" "<<rnd2.Uniform()<<" "<<rnd3.Normal()<<std::endl;
return 0;
}
// -- TestRandom.cpp
) You’ll want to use true random numbers from an outside source. Here’s a list of facilities provided in the OpenBSD operating system, for example.