Best way to seed mt19937_64 for Monte Carlo simulations

Mathhead200 picture Mathhead200 · Jun 20, 2014 · Viewed 7.6k times · Source

I'm working on a program that runs Monte Carlo simulation; specifically, I'm using a Metropolis algorithm. The program needs to generate possibly billions of "random" numbers. I know that the Mersenne twister is very popular for Monte Carlo simulation, but I would like to make sure that I am seeding the generator in the best way possible.

Currently I'm computing a 32-bit seed using the following method:

mt19937_64 prng; //pseudo random number generator
unsigned long seed; //store seed so that every run can follow the same sequence
unsigned char seed_count; //to help keep seeds from repeating because of temporal proximity

unsigned long genSeed() {
    return (  static_cast<unsigned long>(time(NULL))      << 16 )
         | ( (static_cast<unsigned long>(clock()) & 0xFF) << 8  )
         | ( (static_cast<unsigned long>(seed_count++) & 0xFF) );
}

//...

seed = genSeed();
prng.seed(seed);

I have a feeling there are much better ways to assure non-repeating new seeds, and I'm quite sure mt19937_64 can be seeded with more then 32-bits. Does anyone have any suggestions?

Answer

Praetorian picture Praetorian · Jun 20, 2014

Use std::random_device to generate the seed. It'll provide non-deterministic random numbers, provided your implementation supports it. Otherwise it's allowed to use some other random number engine.

std::mt19937_64 prng;
seed = std::random_device{}();
prng.seed(seed);

operator() of std::random_device returns an unsigned int, so if your platform has 32-bit ints, and you want a 64-bit seed, you'll need to call it twice.

std::mt19937_64 prng;
std::random_device device;
seed = (static_cast<uint64_t>(device()) << 32) | device();
prng.seed(seed);

Another available option is using std::seed_seq to seed the PRNG. This allows the PRNG to call seed_seq::generate, which produces a non-biased sequence over the range [0 ≤ i < 232), with an output range large enough to fill its entire state.

std::mt19937_64 prng;
std::random_device device;
std::seed_seq seq{device(), device(), device(), device()};
prng.seed(seq);

I'm calling the random_device 4 times to create a 4 element initial sequence for seed_seq. However, I'm not sure what the best practice for this is, as far as length or source of elements in the initial sequence is concerned.