39
\$\begingroup\$

Many people seed their Mersenne Twister engines like this:

std::mt19937 rng(std::random_device{}());

However, this only provides a single unsigned int, i.e. 32 bits on most systems, of seed randomness, which seems quite tiny when compared to the 19937 bit state space we want to seed. Indeed, if I find out the first number generated, my PC (Intel i7-4790K) only needs about 10 minutes to search through all 32 bit numbers and find the used seed. (I know that MT is not a cryptographic RNG, but I just did that to get a feel for how small 32 bit really is in these days.)

I am trying to build a function to properly seed a mt19937 and came up with this:

#include <algorithm>
#include <iostream>
#include <random>

auto RandomlySeededMersenneTwister () {
    // Magic number 624: The number of unsigned ints the MT uses as state
    std::vector<unsigned int> random_data(624);
    std::random_device source;
    std::generate(begin(random_data), end(random_data), [&](){return source();});
    std::seed_seq seeds(begin(random_data), end(random_data));
    std::mt19937 seededEngine (seeds);
    return seededEngine;
}


int main() {
    auto rng = RandomlySeededMersenneTwister();
    for (int i = 0; i < 10; ++i)
        std::cout << rng() << "\n";    
}

This does look like a safe solution to me; however, I have learned that problems with RNG are often times quite subtle.

Provided std::random_device produces good, random data on my system, does the code give me a correctly seeded std::mt19937?

\$\endgroup\$
  • \$\begingroup\$ As a side note: currently, std::random_device does not work correctly on every platform (I'm looking at you MinGW). If you're going the Boost way, then you could use boost::random_device instead. \$\endgroup\$ – Morwenn Nov 2 '15 at 11:14
32
\$\begingroup\$
  1. Well, first off, why do you use a std::vector for a comparatively small sequence of known length? A raw array or std::array suffice and avoids any dynamic allocation.

  2. Next, avoid needless magic numbers. Use std::mt19937::state_size instead of manually specifying 624.

  3. Why do you use a lambda? A simple std::ref(source) suffices.

The seeding itself looks perfectly fine and there's no actual error anywhere in your code.

template<class T = std::mt19937, std::size_t N = T::state_size * T::word_size>
auto ProperlySeededRandomEngine () -> typename std::enable_if<N, T>::type {
    unsigned random_data[(N - 1) / std::min(sizeof *random_data, 31 / CHAR_BIT + 1) + 1];
    std::random_device source;
    std::generate(std::begin(random_data), std::end(random_data), std::ref(source));
    std::seed_seq seeds(std::begin(random_data), std::end(random_data));
    return T(seeds);
}

You could avoid the need for random_data by using counting and transforming iterators as detailed in "Sequence iterator? Isn't there one in boost?".

This is not simpler, but maybe more efficient:

template<class T = std::mt19937, std::size_t N = T::state_size * T::word_size>
T ProperlySeededRandomEngine () {
    std::random_device source;
    auto make_iter = [&](std::size_t n) {
        return boost::make_transform_iterator(
            boost::counting_iterator<std::size_t>(n), [&](size_t){return source();});
    };
    std::seed_seq seeds(make_iter(0), make_iter((N - 1) / std::min(sizeof(unsigned), 31 / CHAR_BIT + 1) + 1));
    return T(seeds);
}

On coliru

\$\endgroup\$
13
\$\begingroup\$

I wrote more or less exactly the same function for my own use so of course I think it is pretty awesome. ;-)

Two things that I would do differently (only style, not security):

  • Don't hard-code the magic number 624. The std::mersenne_twister_engine template class has a static constexpr member word_size that you can use instead. Likewise, instead of unsigned, prefer using result_type.

  • Consider making the function a template so it can be used for std::mt19937_64 (and maybe other compatible engines) as well.

\$\endgroup\$
  • \$\begingroup\$ Does using monospace for template, class, static, and constexpr improve the clarity of what you're saying? If so, you missed for, do, and using. \$\endgroup\$ – bcrist Oct 31 '15 at 6:47
  • 3
    \$\begingroup\$ @bcrist No, because in the places I used these words, I am using them with their natural language meaning, not to refer to a C++ keyword. If I were talking about a for loop or a using directive, I would have made that clear by using the appropriate font. Likewise, I would write int if I refer to the C++ data type and say “integer” if I refer to the mathematical concept of numbers. \$\endgroup\$ – 5gon12eder Oct 31 '15 at 10:43
  • 4
    \$\begingroup\$ What I was trying to humorously point out is that template class, static, and constexpr are not just references to C++ keywords, they are natural language concepts in the vocabulary of any competent C++ dev. It is the concept that's important, not the token that represents it. Furthermore, in the case of describing mersenne_twister_engine, using class instead of class is actually incorrect. The standard does not mandate that it (or any other std type) use the class keyword. An implementation could just as easily use struct if it so desires. \$\endgroup\$ – bcrist Oct 31 '15 at 16:25
  • 1
    \$\begingroup\$ @bcrist: Kudos for a way to be more precise by getting more colloquial. \$\endgroup\$ – Deduplicator Nov 1 '15 at 18:16
12
\$\begingroup\$

I can offer a different possibility of correctly initializing a mt19937 random number generator:

auto RandomlySeededMersenneTwister () {
    std::mt19937 rng(std::random_device{}());
    rng.discard(700000);
    return rng;
}

According to this paper ("Improved long-period generators based on linear recurrences modulo 2", F. Panneton, P. L'Ecuyer, M. Matsumoto in AVM TOMS Volume 32 Issue 1, March 2006 Pages 1-16), especially figure 4 in chapter 7, one can obtain a mersenne twister state of high quality even when its initialization variables have only a few bits set to a non-zero value. You need to perform about 700000 random number generations (or twists) during/after the initialization phase.

This certainly extends to your case and is much easier to implement. It is, of course, much slower, but random number generator initialization should always only be triggered once in every binary, so this is probably negligible in many cases.

Also, this solution makes it much easier to reproduce results, as it is not needed to save 624 numbers but only one.

Another advantage is the fact, that you can be sure that you get a nice equal distribution. The original solution does not guarantee this property and relies instead on the quality of the source random number generator and its interaction with mt19937.

\$\endgroup\$
  • \$\begingroup\$ Thanks for writing that review. This was my first thought as well and could be trivially implemented using std::discard_block_engine::discard(). \$\endgroup\$ – Edward Nov 2 '15 at 12:46
  • \$\begingroup\$ Fortunately, you understood what I meant (std::mersenne_twister::discard()) rather than what I wrote! That takes less than 5 milliseconds on my 64-bit Linux box running a dual-core processor at 2.3GHz. \$\endgroup\$ – Edward Nov 2 '15 at 15:49
  • 5
    \$\begingroup\$ Burn-in is good, but it does not change the fact that you seed a 624 byte state with just a single integer. \$\endgroup\$ – Johan Lundberg Apr 20 '16 at 19:46
  • 3
    \$\begingroup\$ It does matter. After you seed with a 32 bit int and skip 700k iterations, there is still only 2^32 possible internal states. The paper you are referring to is talking about recovering from a bad internal state, this is not the problem at hand - if you pass a single number to the Mt constructor it uses a seed sequence to fill all its internal state, so you don't start with an almost all zeroes state. But this process (just like discarding) is deterministic, which is the problem the original question is talking about, and discarding does not fix it. \$\endgroup\$ – etarion Jun 30 '16 at 13:01
  • 1
    \$\begingroup\$ @etarion His remark " (I know that MT is not a cryptographic RNG, I just did that to get a feel for how small 32 bit really is in these days.)" also shows that he is not really keen on the safety aspect of his random numbers. And if you are not looking for a safe RNG my solution works just fine (and you can't use a MT anyways...). \$\endgroup\$ – Nils_M Jul 1 '16 at 16:05
3
\$\begingroup\$

Yes, it does technically provide you with a 'correctly' seeded std::mt19937 PRNG. However, the approach is comparatively clumsy - though you're only offering it as example code, to be fair. You can do without explicit magic numbers like (624) - for the required 19937 bits.


header:

#include <random>

class X
{
    ...
    typedef std::mt19937 random_engine;
    mutable random_engine rng; // mutable is a design decision.
    ...
};

implementation:

#include <algorithm>
#include <functional>

X::X ()
{
    std::random_device rdev;
    std::seed_seq::result_type data[random_engine::state_size];
    std::generate_n(data, random_engine::state_size, std::ref(rdev));

    std::seed_seq prng_seed (data, data + random_engine::state_size);
    rng.seed(prng_seed);
}

The member X::rng has a full random 19937-bit state ready for use. Of course this will depend on the policy of std::random_device. Not only do you have an essentially inexhaustible period of (2^19937 - 1) when using rng itself - you also have a random state vector of 19937 bits. The likelihood of re-creating the same initial PRNG state is now ~ (2^-19937), as opposed to (2^-32). Furthermore, all mt19937 states are theoretically reachable with this code.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.