Thread safety of a static random number generator

Question

I have a bunch of threads, each one needs a thread safe random number. Since in my real program threads are spawned and joined repeatedly, I wouldn't like to create random_device and mt19937 each time I enter a new parallel region which calls the same function, so I put them as static:

#include <iostream>
#include <random>
#include <omp.h>

void test(void) {
    static std::random_device rd;
    static std::mt19937 rng(rd());
    static std::uniform_int_distribution<int> uni(1, 1000);

    int x = uni(rng);
#   pragma omp critical
    std::cout << "thread " << omp_get_thread_num() << " | x = " << x << std::endl;
}

int main() {
#   pragma omp parallel num_threads(4)
    test();
}

I cannot place them as threadprivate because of Error C3057: dynamic initialization of 'threadprivate' symbols is not currently supported. Some sources say random_device and mt19937 are thread safe, but I haven't managed to find any docs which would prove it.

1. Is this randomization thread safe?
2. If no, which of the static objects can be left as static to preserve thread safety?

Answer 1

Here is a different approach. I keep a global seeding value so that the random_device is only used once. Since using it can be very slow, I think it is prudent to only use it as rarely as possible.

Instead we increment the seeding value per thread and also per use. That way we avoid the birthday paradox and we minimize the thread-local state to a single integer.

#include <omp.h>

#include <algorithm>
#include <array>
#include <random>


using seed_type = std::array<std::mt19937::result_type, std::mt19937::state_size>;


namespace {

  seed_type init_seed()
  {
    seed_type rtrn;
    std::random_device rdev;
    std::generate(rtrn.begin(), rtrn.end(), std::ref(rdev));
    return rtrn;
  }
  
}
/**
 * Provides a process-global random seeding value
 *
 * Thread-safe (assuming the C++ compiler if standard-conforming.
 * Seed is initialized on first call
 */
seed_type global_seed()
{
  static seed_type rtrn = init_seed();
  return rtrn;
}
/**
 * Creates a new random number generator
 *
 * Operation is thread-safe, Each thread will get its own RNG with a different
 * seed. Repeated calls within a thread will create different RNGs, too.
 */
std::mt19937 make_rng()
{
  static std::mt19937::result_type sequence_number = 0;
# pragma omp threadprivate(sequence_number)
  seed_type seed = global_seed();
  static_assert(seed.size() >= 3);
  seed[0] += sequence_number++;
  seed[1] += static_cast<std::mt19937::result_type>(omp_get_thread_num());
  seed[2] += static_cast<std::mt19937::result_type>(omp_get_level());
  std::seed_seq sseq(seed.begin(), seed.end());
  return std::mt19937(sseq);
}

See also this: How to make this code thread safe with openMP? Monte Carlo two-dimensional integration

For the approach of just increment the seeding value, see this: https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/

Answer 2

I think threadprivate is the right approach still, and you can obviate the initialization problem by doing a parallel assignment later.

static random_device rd;
static mt19937 rng;
#pragma omp threadprivate(rd)
#pragma omp threadprivate(rng)

int main() {

#pragma omp parallel
  rng = mt19937(rd());

#pragma omp parallel
  {
    stringstream res;
    uniform_int_distribution<int> uni(1, 100);
    res << "Thread " << omp_get_thread_num() << ": " << uni(rng) << "\n";
    cout << res.str();
  }

  return 0;
}

Btw, note the stringstream: OpenMP has a tendency to split output lines at the << operators.