Seed std::mt19937 from std::random_device

Question

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?

Answer 1

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 * sizeof(typename T::result_type)>
auto ProperlySeededRandomEngine () -> typename std::enable_if<N, T>::type {
    std::random_device source;
    std::random_device::result_type random_data[(N - 1) / sizeof(source()) + 1];
    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 * sizeof(typename T::result_type)>
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) / sizeof(source()) + 1));
    return T(seeds);
}

On coliru

If you can upgrade to C++20, use ranges and views (godbolt):

template<class T = std::mt19937, std::size_t N = T::state_size * sizeof(typename T::result_type)>
T ProperlySeededRandomEngine () {
    std::random_device source;
    auto random_data = std::views::iota(std::size_t(), (N - 1) / sizeof(source()) + 1)
        | std::views::transform([&](auto){ return source(); });
    std::seed_seq seeds(std::begin(random_data), std::end(random_data));
    return T(seeds);
}

Answer 2

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.

Answer 3

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.

Answer 4

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
{
    ...
    mutable std::mt19937 rng; // mutable is a design decision.
    ...

    // thanks to Konrad Rudolph for correcting me on potentially misleading
    // typedefs, that may be misleading with other RNGs.

};

implementation:

#include <algorithm>
#include <functional>

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

    std::seed_seq prng_seed (data, data + std::mt19937::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.