Is a `std::mt19937` static function variable thread-safe?

Question

I have a function rand_double() which generates a random floating-point value using a static function variable with type std::mt19937. Several threads may use this function concurrently.

#include <random>

auto random_double(double min, double max) {
    static std::mt19937 generator{std::random_device{}()}; // Forgot to seed it, thanks @user4581301
    return std::uniform_real_distribution<>{min, max}(generator);
}

My first question is, is the above code as-is thread safe? I know that since C++11 the initialization of static function variables is always thread-safe, but how about the concurrent use of generator?

Answer 1

To generate a random number, the std::uniform_real_distribution calls operator() of the argument generator.

In case of mersenne twister the operator has the side effect:

"The state of the engine is advanced by one position"

which is not (by any indicator) considered a thread-safe operation.

You can use it as thread_local which would have the intended effect i.e. static in each specific thread:

auto random_double(double min, double max) {
    thread_local std::mt19937 generator(std::random_device{}());
    return std::uniform_real_distribution<>{min, max}(generator);
}

Answer 2

thread_local is the way to go

The OP asks:

Is a std::mt19937 static function variable thread-safe?

In a word, no. As was identified in the comments, calling operator() on a std::mt19937 object is not thread safe.

For applications that make many calls to std::mt19937, declaring the engine to be thread_local is best solution. That way, each thread will get it's own instance of the random number engine.

Using param_type with a random number distribution

By making std::uniform_real_distribution thread_local, as well, you can speed things up a tiny bit . That avoids an unnecessary constructor call every time you generate a random number.

To make this work, you must supply a param_type object to uniform_real_distribution, together with the std::mt19937 object, when you generate a random number.

In the function below, param_type{ min, max } constructs a temporary param object, and initializes it with min and max.

auto random_double(double const min, double const max)
{
    thread_local std::mt19937 mt{ std::random_device{}() };
    thread_local std::uniform_real_distribution<double> dist;
    using param_type = std::uniform_real_distribution<double>::param_type;
    return dist(mt, param_type{ min, max });
}

This is the most direct solution to the question from the OP.

As an alternative, you can speed things up even more by constructing a param object before you call function random_double, and reuse it for multiple invocations.

using param_type = std::uniform_real_distribution<double>::param_type;

auto random_double(param_type const& p)
{
    thread_local std::mt19937 mt{ std::random_device{}() };
    thread_local std::uniform_real_distribution<double> dist;
    return dist(mt, p);  // p is passed by const& straight through to distribution function.
}

void some_function()
{
    // Construct the `param` object once, ...
    param_type p{ -100.0, +100.0 };

    // and use it many times.
    for (int n{ 10 }; n--; )
        std::cout << "   " << random_double(p);
    std::cout.put('\n');
}

What if min and max are constant?

If min and max never change, you make things even more efficient by hard-coding them. This version of function random_double avoids the overhead of passing in the arguments min and max. For the example below, I have set them, respectively, to 1.0 and 10.0.

auto random_double()
{
    thread_local std::mt19937 mt{ std::random_device{}() };
    thread_local std::uniform_real_distribution<double> dist{ 1.0, 10.0 };
    return dist(mt);
}

Use std::next_after if you need a closed interval

The C++ Standard requires that the random numbers returned by a std::uniform_real_distribution object are uniformly distributed over a half-open interval, say [a,b). Thus, they are greater than or equal to a, and less than (but never equal to) b.

To create a distribution over the closed interval [a,b], std::nextafter(b, std::numeric_limits<RealType>::max()) can be used as the second argument of a distribution object's constructor. The same is true for calls to function param and the operator() call function.

Best way to seed std::mt19937 using std::random_device

std::mt19937 and std::19937_64 both have 19,968 bits of state, so trying to seed them with a single seed of only 32 or 64 bits is borderline foolish.

To make it easy to seed all 19,968 bits of state, I wrote class seed_seq_rd. seed_seq_rd stands for seed-sequence using random-device. This class mimics the interface of std::seed_seq, but uses std::random_device to generate seeds. Objects of this type are seed sequences that can be used as arguments to member function seed in a random number engine.

// Example: Seed mt19937 with random seeds from std::random_device.
tbx::seed_seq_rd s;
std::mt19937 mt;
mt.seed( s );

A seed_seq_rd object can also be used as an argument to the constructor of a random number engine, and anywhere else that a std:seed_seq object can be used.

// Example: Seed mt19937_64 at the time of its construction.
tbx::seed_seq_rd s;
std::mt19937_64 mt64{ s };

The souce code for seed_seq_rd is available in a single-file header on GitHub. It works with any random number engine that can be seeded with a seed-sequence, including all of the engines in the Standard Library.

Class seed_seq_rd is defined in header tbx.cpp14.seed_randomly.h. (It uses only the features of C++14, nothing later. Hence the name.) Using it to seed std::mt19937, would change your function to something like the following.

// main.cpp
#include <iostream>
#include <random>
#include "tbx.cpp14.seed_randomly.h"

auto random_double(double const min, double const max)
{
    thread_local tbx::seed_seq_rd s;  // uses only 4 bytes of storage!
    thread_local std::mt19937 mt{ s };
    thread_local std::uniform_real_distribution<double> dist;
    using param_type = std::uniform_real_distribution<double>::param_type;
    return dist(mt, param_type{ min, max });
}

int main()
{
    for (int n{ 10 }; n--; )
        std::cout << "   " << random_double(-5.0, +5.0);
    std::cout.put('\n');
    return 0;
}
// end file: main.cpp

Output:

   -1.95065   4.29782   2.23657   -1.3178   4.73166   -2.22954   -3.14094   4.15405   0.847722   2.52301