diff --git a/demos/demo_imi_flat.cpp b/demos/demo_imi_flat.cpp
index b03781732..72ef16a27 100644
--- a/demos/demo_imi_flat.cpp
+++ b/demos/demo_imi_flat.cpp
@@ -10,6 +10,7 @@
 #include <cmath>
 #include <cstdio>
 #include <cstdlib>
+#include <random>
 
 #include <sys/time.h>
 
@@ -77,6 +78,8 @@ int main ()
     // Use 4-1024 depending on the trade-off speed accuracy that you want
     index.nprobe = 2048;
 
+    std::mt19937 rng;
+    std::uniform_real_distribution<> distrib;
 
     { // training
         printf ("[%.3f s] Generating %ld vectors in %dD for training\n",
@@ -84,7 +87,7 @@ int main ()
 
         std::vector <float> trainvecs (nt * d);
         for (size_t i = 0; i < nt * d; i++) {
-            trainvecs[i] = drand48();
+            trainvecs[i] = distrib(rng);
         }
 
         printf ("[%.3f s] Training the index\n", elapsed() - t0);
@@ -101,7 +104,7 @@ int main ()
 
         std::vector <float> database (nb * d);
         for (size_t i = 0; i < nb * d; i++) {
-            database[i] = drand48();
+            database[i] = distrib(rng);
         }
 
         printf ("[%.3f s] Adding the vectors to the index\n", elapsed() - t0);
diff --git a/demos/demo_imi_pq.cpp b/demos/demo_imi_pq.cpp
index ea6f998c6..9ca11612e 100644
--- a/demos/demo_imi_pq.cpp
+++ b/demos/demo_imi_pq.cpp
@@ -10,6 +10,7 @@
 #include <cmath>
 #include <cstdio>
 #include <cstdlib>
+#include <random>
 
 #include <sys/time.h>
 
@@ -87,6 +88,9 @@ int main ()
     index.nprobe = 2048;
 
 
+    std::mt19937 rng;
+    std::uniform_real_distribution<> distrib;
+
     { // training.
 
         // The distribution of the training vectors should be the same
@@ -100,7 +104,7 @@ int main ()
         std::vector <float> trainvecs (nt * d);
         for (size_t i = 0; i < nt; i++) {
             for (size_t j = 0; j < d; j++) {
-                trainvecs[i * d + j] = drand48();
+                trainvecs[i * d + j] = distrib(rng);
             }
         }
 
@@ -124,7 +128,7 @@ int main ()
         std::vector <long> ids (nb);
         for (size_t i = 0; i < nb; i++) {
             for (size_t j = 0; j < d; j++) {
-                database[i * d + j] = drand48();
+                database[i * d + j] = distrib(rng);
             }
             ids[i] = 8760000000L + i;
         }
diff --git a/demos/demo_ivfpq_indexing.cpp b/demos/demo_ivfpq_indexing.cpp
index 743395ec2..ab2378b27 100644
--- a/demos/demo_ivfpq_indexing.cpp
+++ b/demos/demo_ivfpq_indexing.cpp
@@ -10,6 +10,7 @@
 #include <cmath>
 #include <cstdio>
 #include <cstdlib>
+#include <random>
 
 #include <sys/time.h>
 
@@ -54,13 +55,16 @@ int main ()
                              ncentroids, 4, 8);
 
 
+    std::mt19937 rng;
+
     { // training
         printf ("[%.3f s] Generating %ld vectors in %dD for training\n",
                 elapsed() - t0, nt, d);
 
         std::vector <float> trainvecs (nt * d);
+        std::uniform_real_distribution<> distrib;
         for (size_t i = 0; i < nt * d; i++) {
-            trainvecs[i] = drand48();
+            trainvecs[i] = distrib(rng);
         }
 
         printf ("[%.3f s] Training the index\n",
@@ -86,8 +90,9 @@ int main ()
                 elapsed() - t0, nb);
 
         std::vector <float> database (nb * d);
+        std::uniform_real_distribution<> distrib;
         for (size_t i = 0; i < nb * d; i++) {
-            database[i] = drand48();
+            database[i] = distrib(rng);
         }
 
         printf ("[%.3f s] Adding the vectors to the index\n",
diff --git a/faiss/gpu/test/TestUtils.cpp b/faiss/gpu/test/TestUtils.cpp
index 423d58b87..e3d4bb585 100644
--- a/faiss/gpu/test/TestUtils.cpp
+++ b/faiss/gpu/test/TestUtils.cpp
@@ -24,6 +24,8 @@ inline float relativeError(float a, float b) {
 // This seed is also used for the faiss float_rand API; in a test it
 // is all within a single thread, so it is ok
 long s_seed = 1;
+std::mt19937 rng(1);
+std::uniform_int_distribution<> distrib;
 
 void newTestSeed() {
   struct timespec t;
@@ -35,7 +37,7 @@ void newTestSeed() {
 void setTestSeed(long seed) {
   printf("testing with random seed %ld\n", seed);
 
-  srand48(seed);
+  rng = std::mt19937(seed);
   s_seed = seed;
 }
 
@@ -43,7 +45,7 @@ int randVal(int a, int b) {
   EXPECT_GE(a, 0);
   EXPECT_LE(a, b);
 
-  return a + (lrand48() % (b + 1 - a));
+  return a + (distrib(rng) % (b + 1 - a));
 }
 
 bool randBool() {
diff --git a/faiss/gpu/test/demo_ivfpq_indexing_gpu.cpp b/faiss/gpu/test/demo_ivfpq_indexing_gpu.cpp
index 0f9e6daec..e828ed357 100644
--- a/faiss/gpu/test/demo_ivfpq_indexing_gpu.cpp
+++ b/faiss/gpu/test/demo_ivfpq_indexing_gpu.cpp
@@ -8,6 +8,7 @@
 #include <cmath>
 #include <cstdio>
 #include <cstdlib>
+#include <random>
 
 #include <sys/time.h>
 
@@ -61,13 +62,16 @@ int main ()
     faiss::gpu::GpuIndexIVFPQ index (
       &resources, d, ncentroids, 4, 8, faiss::METRIC_L2, config);
 
+    std::mt19937 rng;
+
     { // training
         printf ("[%.3f s] Generating %ld vectors in %dD for training\n",
                 elapsed() - t0, nt, d);
 
         std::vector <float> trainvecs (nt * d);
+        std::uniform_real_distribution<> distrib;
         for (size_t i = 0; i < nt * d; i++) {
-            trainvecs[i] = drand48();
+            trainvecs[i] = distrib(rng);
         }
 
         printf ("[%.3f s] Training the index\n",
@@ -97,8 +101,9 @@ int main ()
                 elapsed() - t0, nb);
 
         std::vector <float> database (nb * d);
+        std::uniform_real_distribution<> distrib;
         for (size_t i = 0; i < nb * d; i++) {
-            database[i] = drand48();
+            database[i] = distrib(rng);
         }
 
         printf ("[%.3f s] Adding the vectors to the index\n",
diff --git a/misc/test_blas.cpp b/misc/test_blas.cpp
index be2536497..919636ac8 100644
--- a/misc/test_blas.cpp
+++ b/misc/test_blas.cpp
@@ -7,6 +7,7 @@
 
 #include <cstdio>
 #include <cstdlib>
+#include <random>
 
 #undef FINTEGER
 #define FINTEGER long
@@ -31,8 +32,10 @@ int sgeqrf_ (FINTEGER *m, FINTEGER *n, float *a, FINTEGER *lda,
 float *new_random_vec(int size)
 {
     float *x = new float[size];
+    std::mt19937 rng;
+    std::uniform_real_distribution<> distrib;
     for (int i = 0; i < size; i++)
-        x[i] = drand48();
+        x[i] = distrib(rng);
     return x;
 }
 
diff --git a/tests/test_binary_flat.cpp b/tests/test_binary_flat.cpp
index eb20cee87..a0efd5c96 100644
--- a/tests/test_binary_flat.cpp
+++ b/tests/test_binary_flat.cpp
@@ -23,8 +23,6 @@ TEST(BinaryFlat, accuracy) {
   // make the index object and train it
   faiss::IndexBinaryFlat index(d);
 
-  srand(35);
-
   std::vector<uint8_t> database(nb * (d / 8));
   for (size_t i = 0; i < nb * (d / 8); i++) {
     database[i] = rand() % 0x100;
diff --git a/tests/test_dealloc_invlists.cpp b/tests/test_dealloc_invlists.cpp
index d77cd242a..f1cf7f62f 100644
--- a/tests/test_dealloc_invlists.cpp
+++ b/tests/test_dealloc_invlists.cpp
@@ -10,6 +10,7 @@
 
 #include <memory>
 #include <vector>
+#include <random>
 
 #include <gtest/gtest.h>
 
@@ -19,6 +20,7 @@
 #include <faiss/index_io.h>
 #include <faiss/IVFlib.h>
 
+
 using namespace faiss;
 
 namespace {
@@ -38,12 +40,15 @@ size_t nb = 1000;
 // nb of queries
 size_t nq = 200;
 
+std::mt19937 rng;
 
 std::vector<float> make_data(size_t n)
 {
     std::vector <float> database (n * d);
+    std::uniform_real_distribution<> distrib;
+
     for (size_t i = 0; i < n * d; i++) {
-        database[i] = drand48();
+        database[i] = distrib(rng);
     }
     return database;
 }
diff --git a/tests/test_ivfpq_codec.cpp b/tests/test_ivfpq_codec.cpp
index 8d18ac0ad..2732b6787 100644
--- a/tests/test_ivfpq_codec.cpp
+++ b/tests/test_ivfpq_codec.cpp
@@ -7,6 +7,7 @@
 
 #include <cstdio>
 #include <cstdlib>
+#include <random>
 
 #include <gtest/gtest.h>
 
@@ -51,8 +52,10 @@ double eval_codec_error (long ncentroids, long m, const std::vector<float> &v)
 TEST(IVFPQ, codec) {
 
     std::vector <float> database (nb * d);
+    std::mt19937 rng;
+    std::uniform_real_distribution<> distrib;
     for (size_t i = 0; i < nb * d; i++) {
-        database[i] = drand48();
+        database[i] = distrib(rng);
     }
 
     double err0 = eval_codec_error(16, 8, database);
diff --git a/tests/test_ivfpq_indexing.cpp b/tests/test_ivfpq_indexing.cpp
index 9f4bbcd2c..383007831 100644
--- a/tests/test_ivfpq_indexing.cpp
+++ b/tests/test_ivfpq_indexing.cpp
@@ -8,6 +8,7 @@
 
 #include <cstdio>
 #include <cstdlib>
+#include <random>
 
 #include <gtest/gtest.h>
 
@@ -39,13 +40,14 @@ TEST(IVFPQ, accuracy) {
     // index that gives the ground-truth
     faiss::IndexFlatL2 index_gt (d);
 
-    srand48 (35);
+    std::mt19937 rng;
+    std::uniform_real_distribution<> distrib;
 
     { // training
 
         std::vector <float> trainvecs (nt * d);
         for (size_t i = 0; i < nt * d; i++) {
-            trainvecs[i] = drand48();
+            trainvecs[i] = distrib(rng);
         }
         index.verbose = true;
         index.train (nt, trainvecs.data());
@@ -55,7 +57,7 @@ TEST(IVFPQ, accuracy) {
 
         std::vector <float> database (nb * d);
         for (size_t i = 0; i < nb * d; i++) {
-            database[i] = drand48();
+            database[i] = distrib(rng);
         }
 
         index.add (nb, database.data());
@@ -69,7 +71,7 @@ TEST(IVFPQ, accuracy) {
 
         std::vector <float> queries (nq * d);
         for (size_t i = 0; i < nq * d; i++) {
-            queries[i] = drand48();
+            queries[i] = distrib(rng);
         }
 
         std::vector<faiss::Index::idx_t> gt_nns (nq);
diff --git a/tests/test_lowlevel_ivf.cpp b/tests/test_lowlevel_ivf.cpp
index 7fd62eefc..2a3c8f63b 100644
--- a/tests/test_lowlevel_ivf.cpp
+++ b/tests/test_lowlevel_ivf.cpp
@@ -12,6 +12,7 @@
 #include <memory>
 #include <vector>
 #include <thread>
+#include <random>
 
 #include <gtest/gtest.h>
 
@@ -24,6 +25,7 @@
 #include <faiss/IVFlib.h>
 #include <faiss/VectorTransform.h>
 
+
 using namespace faiss;
 
 namespace {
@@ -45,12 +47,15 @@ size_t nq = 200;
 
 int k = 10;
 
+std::mt19937 rng;
+
 
 std::vector<float> make_data(size_t n)
 {
     std::vector <float> database (n * d);
+    std::uniform_real_distribution<> distrib;
     for (size_t i = 0; i < n * d; i++) {
-        database[i] = drand48();
+        database[i] = distrib(rng);
     }
     return database;
 }
@@ -274,8 +279,9 @@ std::vector<uint8_t> make_data_binary(size_t n)
 {
 
     std::vector <uint8_t> database (n * nbit / 8);
+    std::uniform_int_distribution<> distrib;
     for (size_t i = 0; i < n * d; i++) {
-        database[i] = lrand48();
+        database[i] = distrib(rng);
     }
     return database;
 }
diff --git a/tests/test_merge.cpp b/tests/test_merge.cpp
index 47af10614..edd89a9ab 100644
--- a/tests/test_merge.cpp
+++ b/tests/test_merge.cpp
@@ -7,6 +7,7 @@
 
 #include <cstdio>
 #include <cstdlib>
+#include <random>
 
 #include <gtest/gtest.h>
 
@@ -70,11 +71,13 @@ struct CommonData {
 
     CommonData(): database (nb * d), queries (nq * d), ids(nb), quantizer (d) {
 
+        std::mt19937 rng;
+        std::uniform_real_distribution<> distrib;
         for (size_t i = 0; i < nb * d; i++) {
-            database[i] = drand48();
+            database[i] = distrib(rng);
         }
         for (size_t i = 0; i < nq * d; i++) {
-            queries[i] = drand48();
+            queries[i] = distrib(rng);
         }
         for (int i = 0; i < nb; i++) {
             ids[i] = 123 + 456 * i;
diff --git a/tests/test_ondisk_ivf.cpp b/tests/test_ondisk_ivf.cpp
index c7f717faf..cf5f97e77 100644
--- a/tests/test_ondisk_ivf.cpp
+++ b/tests/test_ondisk_ivf.cpp
@@ -7,6 +7,7 @@
 
 #include <cstdio>
 #include <cstdlib>
+#include <random>
 
 #include <omp.h>
 
@@ -69,8 +70,10 @@ TEST(ONDISK, make_invlists) {
 
     {
         std::vector<uint8_t> code(32);
+        std::mt19937 rng;
+        std::uniform_real_distribution<> distrib;
         for (int i = 0; i < nadd; i++) {
-            double d = drand48();
+            double d = distrib(rng);
             int list_no = int(nlist * d * d); // skewed distribution
             int * ar = (int*)code.data();
             ar[0] = i;
@@ -183,8 +186,10 @@ TEST(ONDISK, make_invlists_threaded) {
 
     std::vector<int> list_nos (nadd);
 
+    std::mt19937 rng;
+    std::uniform_real_distribution<> distrib;
     for (int i = 0; i < nadd; i++) {
-        double d = drand48();
+        double d = distrib(rng);
         list_nos[i] = int(nlist * d * d); // skewed distribution
     }
 
diff --git a/tests/test_pairs_decoding.cpp b/tests/test_pairs_decoding.cpp
index 7857d0fb5..2318b3562 100644
--- a/tests/test_pairs_decoding.cpp
+++ b/tests/test_pairs_decoding.cpp
@@ -10,6 +10,7 @@
 
 #include <memory>
 #include <vector>
+#include <random>
 
 #include <gtest/gtest.h>
 
@@ -37,11 +38,14 @@ size_t nb = 8000;
 // nb of queries
 size_t nq = 200;
 
+std::mt19937 rng;
+
 std::vector<float> make_data(size_t n)
 {
     std::vector <float> database (n * d);
+    std::uniform_real_distribution<> distrib;
     for (size_t i = 0; i < n * d; i++) {
-        database[i] = drand48();
+        database[i] = distrib(rng);
     }
     return database;
 }
diff --git a/tests/test_params_override.cpp b/tests/test_params_override.cpp
index d6df2a4ef..352189d9c 100644
--- a/tests/test_params_override.cpp
+++ b/tests/test_params_override.cpp
@@ -10,6 +10,7 @@
 
 #include <memory>
 #include <vector>
+#include <random>
 
 #include <gtest/gtest.h>
 
@@ -19,6 +20,7 @@
 #include <faiss/AutoTune.h>
 #include <faiss/IVFlib.h>
 
+
 using namespace faiss;
 
 namespace {
@@ -35,13 +37,15 @@ size_t nb = 1000;
 // nb of queries
 size_t nq = 200;
 
+std::mt19937 rng;
 
 
 std::vector<float> make_data(size_t n)
 {
     std::vector <float> database (n * d);
+    std::uniform_real_distribution<> distrib;
     for (size_t i = 0; i < n * d; i++) {
-        database[i] = drand48();
+        database[i] = distrib(rng);
     }
     return database;
 }
@@ -155,8 +159,9 @@ TEST(TPO, IVFFlatPP) {
 
 std::vector<uint8_t> make_data_binary(size_t n) {
     std::vector <uint8_t> database (n * d / 8);
+    std::uniform_int_distribution<> distrib;
     for (size_t i = 0; i < n * d / 8; i++) {
-        database[i] = lrand48();
+        database[i] = distrib(rng);
     }
     return database;
 }
diff --git a/tests/test_sliding_ivf.cpp b/tests/test_sliding_ivf.cpp
index ea77dd9b4..bf1b70ade 100644
--- a/tests/test_sliding_ivf.cpp
+++ b/tests/test_sliding_ivf.cpp
@@ -9,6 +9,7 @@
 #include <cstdlib>
 
 #include <memory>
+#include <random>
 #include <vector>
 
 #include <gtest/gtest.h>
@@ -19,6 +20,7 @@
 #include <faiss/clone_index.h>
 #include <faiss/IVFlib.h>
 
+
 using namespace faiss;
 
 typedef Index::idx_t idx_t;
@@ -47,8 +49,10 @@ int window_size = 10;
 std::vector<float> make_data(size_t n)
 {
     std::vector <float> database (n * d);
+    std::mt19937 rng;
+    std::uniform_real_distribution<> distrib;
     for (size_t i = 0; i < n * d; i++) {
-        database[i] = drand48();
+        database[i] = distrib(rng);
     }
     return database;
 }
@@ -92,8 +96,10 @@ void make_index_slices (const Index* trained_index,
 
         auto xb = make_data(nb * d);
         std::vector<faiss::Index::idx_t> ids (nb);
+        std::mt19937 rng;
+        std::uniform_int_distribution<> distrib;
         for (int j = 0; j < nb; j++) {
-            ids[j] = lrand48();
+            ids[j] = distrib(rng);
         }
         index->add_with_ids (nb, xb.data(), ids.data());
     }
diff --git a/tutorial/cpp/1-Flat.cpp b/tutorial/cpp/1-Flat.cpp
index 2f25436ee..e4a388ba4 100644
--- a/tutorial/cpp/1-Flat.cpp
+++ b/tutorial/cpp/1-Flat.cpp
@@ -7,6 +7,7 @@
 
 #include <cstdio>
 #include <cstdlib>
+#include <random>
 
 #include <faiss/IndexFlat.h>
 
@@ -18,18 +19,21 @@ int main() {
     int nb = 100000;                       // database size
     int nq = 10000;                        // nb of queries
 
+    std::mt19937 rng;
+    std::uniform_real_distribution<> distrib;
+
     float *xb = new float[d * nb];
     float *xq = new float[d * nq];
 
     for(int i = 0; i < nb; i++) {
         for(int j = 0; j < d; j++)
-            xb[d * i + j] = drand48();
+            xb[d * i + j] = distrib(rng);
         xb[d * i] += i / 1000.;
     }
 
     for(int i = 0; i < nq; i++) {
         for(int j = 0; j < d; j++)
-            xq[d * i + j] = drand48();
+            xq[d * i + j] = distrib(rng);
         xq[d * i] += i / 1000.;
     }
 
diff --git a/tutorial/cpp/2-IVFFlat.cpp b/tutorial/cpp/2-IVFFlat.cpp
index d0bf4ff92..c778eaf0f 100644
--- a/tutorial/cpp/2-IVFFlat.cpp
+++ b/tutorial/cpp/2-IVFFlat.cpp
@@ -8,6 +8,7 @@
 #include <cstdio>
 #include <cstdlib>
 #include <cassert>
+#include <random>
 
 #include <faiss/IndexFlat.h>
 #include <faiss/IndexIVFFlat.h>
@@ -19,18 +20,21 @@ int main() {
     int nb = 100000;                       // database size
     int nq = 10000;                        // nb of queries
 
+    std::mt19937 rng;
+    std::uniform_real_distribution<> distrib;
+
     float *xb = new float[d * nb];
     float *xq = new float[d * nq];
 
     for(int i = 0; i < nb; i++) {
         for(int j = 0; j < d; j++)
-            xb[d * i + j] = drand48();
+            xb[d * i + j] = distrib(rng);
         xb[d * i] += i / 1000.;
     }
 
     for(int i = 0; i < nq; i++) {
         for(int j = 0; j < d; j++)
-            xq[d * i + j] = drand48();
+            xq[d * i + j] = distrib(rng);
         xq[d * i] += i / 1000.;
     }
 
diff --git a/tutorial/cpp/3-IVFPQ.cpp b/tutorial/cpp/3-IVFPQ.cpp
index 7569f899d..5165a9ef4 100644
--- a/tutorial/cpp/3-IVFPQ.cpp
+++ b/tutorial/cpp/3-IVFPQ.cpp
@@ -7,6 +7,7 @@
 
 #include <cstdio>
 #include <cstdlib>
+#include <random>
 
 #include <faiss/IndexFlat.h>
 #include <faiss/IndexIVFPQ.h>
@@ -18,18 +19,21 @@ int main() {
     int nb = 100000;                       // database size
     int nq = 10000;                        // nb of queries
 
+    std::mt19937 rng;
+    std::uniform_real_distribution<> distrib;
+
     float *xb = new float[d * nb];
     float *xq = new float[d * nq];
 
     for(int i = 0; i < nb; i++) {
         for(int j = 0; j < d; j++)
-            xb[d * i + j] = drand48();
+            xb[d * i + j] = distrib(rng);
         xb[d * i] += i / 1000.;
     }
 
     for(int i = 0; i < nq; i++) {
         for(int j = 0; j < d; j++)
-            xq[d * i + j] = drand48();
+            xq[d * i + j] = distrib(rng);
         xq[d * i] += i / 1000.;
     }
 
diff --git a/tutorial/cpp/4-GPU.cpp b/tutorial/cpp/4-GPU.cpp
index 9737971ae..ec34c54cd 100644
--- a/tutorial/cpp/4-GPU.cpp
+++ b/tutorial/cpp/4-GPU.cpp
@@ -8,6 +8,7 @@
 #include <cstdio>
 #include <cstdlib>
 #include <cassert>
+#include <random>
 
 #include <faiss/IndexFlat.h>
 #include <faiss/gpu/GpuIndexFlat.h>
@@ -20,18 +21,21 @@ int main() {
     int nb = 100000;                       // database size
     int nq = 10000;                        // nb of queries
 
+    std::mt19937 rng;
+    std::uniform_real_distribution<> distrib;
+
     float *xb = new float[d * nb];
     float *xq = new float[d * nq];
 
     for(int i = 0; i < nb; i++) {
         for(int j = 0; j < d; j++)
-            xb[d * i + j] = drand48();
+            xb[d * i + j] = distrib(rng);
         xb[d * i] += i / 1000.;
     }
 
     for(int i = 0; i < nq; i++) {
         for(int j = 0; j < d; j++)
-            xq[d * i + j] = drand48();
+            xq[d * i + j] = distrib(rng);
         xq[d * i] += i / 1000.;
     }
 
diff --git a/tutorial/cpp/5-Multiple-GPUs.cpp b/tutorial/cpp/5-Multiple-GPUs.cpp
index 3152b731a..852f59191 100644
--- a/tutorial/cpp/5-Multiple-GPUs.cpp
+++ b/tutorial/cpp/5-Multiple-GPUs.cpp
@@ -7,6 +7,7 @@
 
 #include <cstdio>
 #include <cstdlib>
+#include <random>
 
 #include <faiss/IndexFlat.h>
 #include <faiss/gpu/GpuAutoTune.h>
@@ -21,18 +22,21 @@ int main() {
     int nb = 100000;                       // database size
     int nq = 10000;                        // nb of queries
 
+    std::mt19937 rng;
+    std::uniform_real_distribution<> distrib;
+
     float *xb = new float[d * nb];
     float *xq = new float[d * nq];
 
     for(int i = 0; i < nb; i++) {
         for(int j = 0; j < d; j++)
-            xb[d * i + j] = drand48();
+            xb[d * i + j] = distrib(rng);
         xb[d * i] += i / 1000.;
     }
 
     for(int i = 0; i < nq; i++) {
         for(int j = 0; j < d; j++)
-            xq[d * i + j] = drand48();
+            xq[d * i + j] = distrib(rng);
         xq[d * i] += i / 1000.;
     }