Pasted

diff --git a/utils/rand_generator.h b/utils/rand_generator.h
index 45e09d0..eb7f41d 100644
--- a/utils/rand_generator.h
+++ b/utils/rand_generator.h
@@ -2,6 +2,7 @@
 #define _RAND_GENERATOR_H
 
 #include <stdlib.h>
+#include <random>
 #include <vector>
 
 #ifdef OPENMP
@@ -17,8 +18,23 @@ class rand_generator
   public:
   rand_generator(size_t rsize) : _rands(rsize), _ridx(0)
   {
+    /*
+     * We use a single seed value that is incremented per thread.
+     * See here for motivation:
+     * https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/
+     */
+    static unsigned global_seed = std::random_device()();
+    unsigned thread_seed;
+#   pragma omp atomic relaxed capture
+    thread_seed = global_seed++;
+    /*
+     * Consider using a better random number generator. This one should resemble
+     * drand48(), which was used previously
+     */
+    std::minstd_rand rng(thread_seed);
+    std::uniform_real_distribution<float> distr;
     for(size_t i=0; i<_rands.size(); i++)
-      _rands[i] = drand48();
+      _rands[i] = distr(rng);
   }
 
   rand_generator() : rand_generator(4*1024*1024)
@@ -32,7 +48,11 @@ class rand_generator
 };
 
 class rand_generator_state {
-  mt_vector<rand_generator> global_rands;
+  /*
+   * use std::vector instead of mt_vector because we have to avoid calling
+   * new rand_generator[N] which would invoke the costly default constructor
+   */
+  std::vector<rand_generator> global_rands;
 
   public:
   rand_generator_state() {
@@ -41,8 +61,13 @@ class rand_generator_state {
     {
       int tid      = omp_get_thread_num();
       int nthreads = omp_get_num_threads();
-      if (tid == 0)
-        global_rands.resize(nthreads);
+      /*
+       * Create a vector of empty generators, then let each thread
+       * initialize its own
+       */
+#     pragma omp single
+      global_rands.resize(nthreads, rand_generator(0));
+      global_rands[tid] = rand_generator();
     }
 #else
     global_rands.resize(1);