diff --git a/internal/ceres/autodiff_test.cc b/internal/ceres/autodiff_test.cc
index 629eea8..0916a45 100644
--- a/internal/ceres/autodiff_test.cc
+++ b/internal/ceres/autodiff_test.cc
@@ -30,6 +30,8 @@
 
 #include "ceres/internal/autodiff.h"
 
+#include <algorithm>
+#include <iterator>
 #include <random>
 
 #include "gtest/gtest.h"
@@ -173,9 +175,9 @@
 
   // Make random P and X, in a single vector.
   double PX[12 + 4];
-  for (double& PX_i : PX) {
-    PX_i = uniform01(prng);
-  }
+  std::generate(std::begin(PX), std::end(PX), [&prng, &uniform01] {
+    return uniform01(prng);
+  });
 
   // Handy names for the P and X parts.
   double* P = PX + 0;
@@ -294,7 +296,10 @@
   double qcX[4 + 3 + 3];
   std::mt19937 prng;
   std::uniform_real_distribution<double> uniform01(0.0, 1.0);
-  for (double& qcX_i : qcX) qcX_i = uniform01(prng);
+
+  std::generate(std::begin(qcX), std::end(qcX), [&prng, &uniform01] {
+    return uniform01(prng);
+  });
 
   // Handy names.
   double* q = qcX;
diff --git a/internal/ceres/block_sparse_matrix.cc b/internal/ceres/block_sparse_matrix.cc
index ae6bd3a..90bd926 100644
--- a/internal/ceres/block_sparse_matrix.cc
+++ b/internal/ceres/block_sparse_matrix.cc
@@ -375,7 +375,7 @@
 }
 
 std::unique_ptr<BlockSparseMatrix> BlockSparseMatrix::CreateRandomMatrix(
-    const BlockSparseMatrix::RandomMatrixOptions& options) {
+    const BlockSparseMatrix::RandomMatrixOptions& options, std::mt19937& prng) {
   CHECK_GT(options.num_row_blocks, 0);
   CHECK_GT(options.min_row_block_size, 0);
   CHECK_GT(options.max_row_block_size, 0);
@@ -383,7 +383,6 @@
   CHECK_GT(options.block_density, 0.0);
   CHECK_LE(options.block_density, 1.0);
 
-  std::mt19937 prng;
   std::uniform_int_distribution<int> col_distribution(
       options.min_col_block_size, options.max_col_block_size);
   std::uniform_int_distribution<int> row_distribution(
@@ -436,9 +435,10 @@
   auto matrix = std::make_unique<BlockSparseMatrix>(bs);
   double* values = matrix->mutable_values();
   std::normal_distribution<double> standard_normal_distribution;
-  for (int i = 0; i < matrix->num_nonzeros(); ++i) {
-    values[i] = standard_normal_distribution(prng);
-  }
+  std::generate_n(
+      values, matrix->num_nonzeros(), [&standard_normal_distribution, &prng] {
+        return standard_normal_distribution(prng);
+      });
 
   return matrix;
 }
diff --git a/internal/ceres/block_sparse_matrix.h b/internal/ceres/block_sparse_matrix.h
index 7cef18d..44cb5e5 100644
--- a/internal/ceres/block_sparse_matrix.h
+++ b/internal/ceres/block_sparse_matrix.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -35,6 +35,7 @@
 #define CERES_INTERNAL_BLOCK_SPARSE_MATRIX_H_
 
 #include <memory>
+#include <random>
 
 #include "ceres/block_structure.h"
 #include "ceres/crs_matrix.h"
@@ -123,7 +124,7 @@
   // distributed and whose structure is determined by
   // RandomMatrixOptions.
   static std::unique_ptr<BlockSparseMatrix> CreateRandomMatrix(
-      const RandomMatrixOptions& options);
+      const RandomMatrixOptions& options, std::mt19937& prng);
 
  private:
   int num_rows_;
diff --git a/internal/ceres/compressed_row_sparse_matrix.cc b/internal/ceres/compressed_row_sparse_matrix.cc
index 574b0c0..ce587da 100644
--- a/internal/ceres/compressed_row_sparse_matrix.cc
+++ b/internal/ceres/compressed_row_sparse_matrix.cc
@@ -119,11 +119,12 @@
   transpose_rows[0] = 0;
 }
 
+template <class RandomNormalFunctor>
 void AddRandomBlock(const int num_rows,
                     const int num_cols,
                     const int row_block_begin,
                     const int col_block_begin,
-                    std::function<double()> dist,
+                    RandomNormalFunctor&& randn,
                     std::vector<int>* rows,
                     std::vector<int>* cols,
                     std::vector<double>* values) {
@@ -131,20 +132,21 @@
     for (int c = 0; c < num_cols; ++c) {
       rows->push_back(row_block_begin + r);
       cols->push_back(col_block_begin + c);
-      values->push_back(dist());
+      values->push_back(randn());
     }
   }
 }
 
+template <class RandomNormalFunctor>
 void AddSymmetricRandomBlock(const int num_rows,
                              const int row_block_begin,
-                             std::function<double()> dist,
+                             RandomNormalFunctor&& randn,
                              std::vector<int>* rows,
                              std::vector<int>* cols,
                              std::vector<double>* values) {
   for (int r = 0; r < num_rows; ++r) {
     for (int c = r; c < num_rows; ++c) {
-      const double v = dist();
+      const double v = randn();
       rows->push_back(row_block_begin + r);
       cols->push_back(row_block_begin + c);
       values->push_back(v);
@@ -625,7 +627,8 @@
 
 std::unique_ptr<CompressedRowSparseMatrix>
 CompressedRowSparseMatrix::CreateRandomMatrix(
-    CompressedRowSparseMatrix::RandomMatrixOptions options) {
+    CompressedRowSparseMatrix::RandomMatrixOptions options,
+    std::mt19937& prng) {
   CHECK_GT(options.num_row_blocks, 0);
   CHECK_GT(options.min_row_block_size, 0);
   CHECK_GT(options.max_row_block_size, 0);
@@ -651,7 +654,6 @@
   vector<int> col_blocks;
   col_blocks.reserve(options.num_col_blocks);
 
-  std::mt19937 prng;
   std::uniform_int_distribution<int> col_distribution(
       options.min_col_block_size, options.max_col_block_size);
   std::uniform_int_distribution<int> row_distribution(
@@ -708,6 +710,9 @@
 
         // Randomly determine if this block is present or not.
         if (uniform01(prng) <= options.block_density) {
+          auto randn = [&standard_normal_distribution, &prng] {
+            return standard_normal_distribution(prng);
+          };
           // If the matrix is symmetric, then we take care to generate
           // symmetric diagonal blocks.
           if (options.storage_type == StorageType::UNSYMMETRIC || r != c) {
@@ -715,14 +720,14 @@
                            col_blocks[c],
                            row_block_begin,
                            col_block_begin,
-                           values_dist,
+                           randn,
                            &tsm_rows,
                            &tsm_cols,
                            &tsm_values);
           } else {
             AddSymmetricRandomBlock(row_blocks[r],
                                     row_block_begin,
-                                    values_dist,
+                                    randn,
                                     &tsm_rows,
                                     &tsm_cols,
                                     &tsm_values);
diff --git a/internal/ceres/compressed_row_sparse_matrix.h b/internal/ceres/compressed_row_sparse_matrix.h
index 2580045..d164bef 100644
--- a/internal/ceres/compressed_row_sparse_matrix.h
+++ b/internal/ceres/compressed_row_sparse_matrix.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -32,6 +32,7 @@
 #define CERES_INTERNAL_COMPRESSED_ROW_SPARSE_MATRIX_H_
 
 #include <memory>
+#include <random>
 #include <vector>
 
 #include "ceres/internal/disable_warnings.h"
@@ -195,7 +196,7 @@
   // normally distributed and whose structure is determined by
   // RandomMatrixOptions.
   static std::unique_ptr<CompressedRowSparseMatrix> CreateRandomMatrix(
-      RandomMatrixOptions options);
+      RandomMatrixOptions options, std::mt19937& prng);
 
  private:
   static std::unique_ptr<CompressedRowSparseMatrix> FromTripletSparseMatrix(
diff --git a/internal/ceres/compressed_row_sparse_matrix_test.cc b/internal/ceres/compressed_row_sparse_matrix_test.cc
index 42f5498..718a659 100644
--- a/internal/ceres/compressed_row_sparse_matrix_test.cc
+++ b/internal/ceres/compressed_row_sparse_matrix_test.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -328,6 +328,7 @@
 }
 
 TEST(CompressedRowSparseMatrix, FromTripletSparseMatrix) {
+  std::mt19937 prng;
   TripletSparseMatrix::RandomMatrixOptions options;
   options.num_rows = 5;
   options.num_cols = 7;
@@ -336,7 +337,7 @@
   const int kNumTrials = 10;
   for (int i = 0; i < kNumTrials; ++i) {
     std::unique_ptr<TripletSparseMatrix> tsm =
-        TripletSparseMatrix::CreateRandomMatrix(options);
+        TripletSparseMatrix::CreateRandomMatrix(options, prng);
     std::unique_ptr<CompressedRowSparseMatrix> crsm =
         CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm);
 
@@ -354,6 +355,7 @@
 }
 
 TEST(CompressedRowSparseMatrix, FromTripletSparseMatrixTransposed) {
+  std::mt19937 prng;
   TripletSparseMatrix::RandomMatrixOptions options;
   options.num_rows = 5;
   options.num_cols = 7;
@@ -362,7 +364,7 @@
   const int kNumTrials = 10;
   for (int i = 0; i < kNumTrials; ++i) {
     std::unique_ptr<TripletSparseMatrix> tsm =
-        TripletSparseMatrix::CreateRandomMatrix(options);
+        TripletSparseMatrix::CreateRandomMatrix(options, prng);
     std::unique_ptr<CompressedRowSparseMatrix> crsm =
         CompressedRowSparseMatrix::FromTripletSparseMatrixTransposed(*tsm);
 
@@ -421,7 +423,7 @@
       options.block_density = uniform(prng);
       options.storage_type = ::testing::get<0>(param);
       std::unique_ptr<CompressedRowSparseMatrix> matrix =
-          CompressedRowSparseMatrix::CreateRandomMatrix(options);
+          CompressedRowSparseMatrix::CreateRandomMatrix(options, prng);
       const int num_rows = matrix->num_rows();
       const int num_cols = matrix->num_cols();
 
@@ -491,7 +493,7 @@
       options.block_density = uniform(prng);
       options.storage_type = ::testing::get<0>(param);
       std::unique_ptr<CompressedRowSparseMatrix> matrix =
-          CompressedRowSparseMatrix::CreateRandomMatrix(options);
+          CompressedRowSparseMatrix::CreateRandomMatrix(options, prng);
       const int num_rows = matrix->num_rows();
       const int num_cols = matrix->num_cols();
 
@@ -561,7 +563,7 @@
       options.block_density = uniform(prng);
       options.storage_type = ::testing::get<0>(param);
       std::unique_ptr<CompressedRowSparseMatrix> matrix =
-          CompressedRowSparseMatrix::CreateRandomMatrix(options);
+          CompressedRowSparseMatrix::CreateRandomMatrix(options, prng);
       const int num_cols = matrix->num_cols();
 
       Vector actual(num_cols);
diff --git a/internal/ceres/corrector_test.cc b/internal/ceres/corrector_test.cc
index 7f58efa..a5dd29c 100644
--- a/internal/ceres/corrector_test.cc
+++ b/internal/ceres/corrector_test.cc
@@ -32,7 +32,6 @@
 
 #include <algorithm>
 #include <cmath>
-#include <cstdlib>
 #include <cstring>
 #include <random>
 
diff --git a/internal/ceres/gradient_checker_test.cc b/internal/ceres/gradient_checker_test.cc
index 3c5f5da..8d7d10f 100644
--- a/internal/ceres/gradient_checker_test.cc
+++ b/internal/ceres/gradient_checker_test.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2021 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -33,7 +33,6 @@
 #include "ceres/gradient_checker.h"
 
 #include <cmath>
-#include <cstdlib>
 #include <random>
 #include <utility>
 #include <vector>
@@ -59,16 +58,16 @@
 // version, they are both block vectors, of course.
 class GoodTestTerm : public CostFunction {
  public:
-  GoodTestTerm(int arity, int const* dim) : arity_(arity), return_value_(true) {
-    std::mt19937 prng;
-    std::uniform_real_distribution<double> distribution(-1.0, 1.0);
-
+  template <class UniformRandomFunctor>
+  GoodTestTerm(int arity, int const* dim, UniformRandomFunctor&& randu)
+      : arity_(arity), return_value_(true) {
+    std::uniform_real_distribution distribution(-1.0, 1.0);
     // Make 'arity' random vectors.
     a_.resize(arity_);
     for (int j = 0; j < arity_; ++j) {
       a_[j].resize(dim[j]);
       for (int u = 0; u < dim[j]; ++u) {
-        a_[j][u] = distribution(prng);
+        a_[j][u] = randu();
       }
     }
 
@@ -121,15 +120,15 @@
 
 class BadTestTerm : public CostFunction {
  public:
-  BadTestTerm(int arity, int const* dim) : arity_(arity) {
-    std::mt19937 prng;
-    std::uniform_real_distribution<double> distribution(-1.0, 1.0);
+  template <class UniformRandomFunctor>
+  BadTestTerm(int arity, int const* dim, UniformRandomFunctor&& randu)
+      : arity_(arity) {
     // Make 'arity' random vectors.
     a_.resize(arity_);
     for (int j = 0; j < arity_; ++j) {
       a_[j].resize(dim[j]);
       for (int u = 0; u < dim[j]; ++u) {
-        a_[j][u] = distribution(prng);
+        a_[j][u] = randu();
       }
     }
 
@@ -210,10 +209,11 @@
   FixedArray<double*> parameters(num_parameters);
   std::mt19937 prng;
   std::uniform_real_distribution<double> distribution(-1.0, 1.0);
+  auto randu = [&prng, &distribution] { return distribution(prng); };
   for (int j = 0; j < num_parameters; ++j) {
     parameters[j] = new double[parameter_sizes[j]];
     for (int u = 0; u < parameter_sizes[j]; ++u) {
-      parameters[j][u] = distribution(prng);
+      parameters[j][u] = randu();
     }
   }
 
@@ -221,7 +221,7 @@
   GradientChecker::ProbeResults results;
 
   // Test that Probe returns true for correct Jacobians.
-  GoodTestTerm good_term(num_parameters, parameter_sizes.data());
+  GoodTestTerm good_term(num_parameters, parameter_sizes.data(), randu);
   std::vector<const Manifold*>* manifolds = nullptr;
   GradientChecker good_gradient_checker(
       &good_term, manifolds, numeric_diff_options);
@@ -258,7 +258,7 @@
   EXPECT_FALSE(results.error_log.empty());
 
   // Test that Probe returns false for incorrect Jacobians.
-  BadTestTerm bad_term(num_parameters, parameter_sizes.data());
+  BadTestTerm bad_term(num_parameters, parameter_sizes.data(), randu);
   GradientChecker bad_gradient_checker(
       &bad_term, manifolds, numeric_diff_options);
   EXPECT_FALSE(
diff --git a/internal/ceres/gradient_checking_cost_function_test.cc b/internal/ceres/gradient_checking_cost_function_test.cc
index caf4f25..87c3278 100644
--- a/internal/ceres/gradient_checking_cost_function_test.cc
+++ b/internal/ceres/gradient_checking_cost_function_test.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -69,15 +69,15 @@
  public:
   // The constructor of this function needs to know the number
   // of blocks desired, and the size of each block.
-  TestTerm(int arity, int const* dim) : arity_(arity) {
-    std::mt19937 prng;
-    std::uniform_real_distribution<double> uniform01(-1.0, 1.0);
+  template <class UniformRandomFunctor>
+  TestTerm(int arity, int const* dim, UniformRandomFunctor&& randu)
+      : arity_(arity) {
     // Make 'arity' random vectors.
     a_.resize(arity_);
     for (int j = 0; j < arity_; ++j) {
       a_[j].resize(dim[j]);
       for (int u = 0; u < dim[j]; ++u) {
-        a_[j][u] = uniform01(prng);
+        a_[j][u] = randu();
       }
     }
 
@@ -139,11 +139,12 @@
   // Make a random set of blocks.
   vector<double*> parameters(arity);
   std::mt19937 prng;
-  std::uniform_real_distribution<double> uniform01(-1.0, 1.0);
+  std::uniform_real_distribution<double> distribution(-1.0, 1.0);
+  auto randu = [&prng, &distribution] { return distribution(prng); };
   for (int j = 0; j < arity; ++j) {
     parameters[j] = new double[dim[j]];
     for (int u = 0; u < dim[j]; ++u) {
-      parameters[j][u] = uniform01(prng);
+      parameters[j][u] = randu();
     }
   }
 
@@ -162,7 +163,7 @@
   const double kRelativeStepSize = 1e-6;
   const double kRelativePrecision = 1e-4;
 
-  TestTerm<-1, -1> term(arity, dim);
+  TestTerm<-1, -1> term(arity, dim, randu);
   GradientCheckingIterationCallback callback;
   auto gradient_checking_cost_function =
       CreateGradientCheckingCostFunction(&term,
@@ -196,11 +197,12 @@
   // Make a random set of blocks.
   vector<double*> parameters(arity);
   std::mt19937 prng;
-  std::uniform_real_distribution<double> uniform01(-1.0, 1.0);
+  std::uniform_real_distribution<double> distribution(-1.0, 1.0);
+  auto randu = [&prng, &distribution] { return distribution(prng); };
   for (int j = 0; j < arity; ++j) {
     parameters[j] = new double[dim[j]];
     for (int u = 0; u < dim[j]; ++u) {
-      parameters[j][u] = uniform01(prng);
+      parameters[j][u] = randu();
     }
   }
 
@@ -218,7 +220,7 @@
   // Should have one term that's bad, causing everything to get dumped.
   LOG(INFO) << "Bad gradient";
   {
-    TestTerm<1, 2> term(arity, dim);
+    TestTerm<1, 2> term(arity, dim, randu);
     GradientCheckingIterationCallback callback;
     auto gradient_checking_cost_function =
         CreateGradientCheckingCostFunction(&term,
@@ -238,7 +240,7 @@
   // The gradient is correct, so no errors are reported.
   LOG(INFO) << "Good gradient";
   {
-    TestTerm<-1, -1> term(arity, dim);
+    TestTerm<-1, -1> term(arity, dim, randu);
     GradientCheckingIterationCallback callback;
     auto gradient_checking_cost_function =
         CreateGradientCheckingCostFunction(&term,
diff --git a/internal/ceres/inner_product_computer_test.cc b/internal/ceres/inner_product_computer_test.cc
index 5bd4c50..4573b6b 100644
--- a/internal/ceres/inner_product_computer_test.cc
+++ b/internal/ceres/inner_product_computer_test.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2017 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -109,7 +109,7 @@
         VLOG(2) << "block density: " << options.block_density;
 
         std::unique_ptr<BlockSparseMatrix> random_matrix(
-            BlockSparseMatrix::CreateRandomMatrix(options));
+            BlockSparseMatrix::CreateRandomMatrix(options, prng));
 
         TripletSparseMatrix tsm(random_matrix->num_rows(),
                                 random_matrix->num_cols(),
@@ -169,7 +169,7 @@
     VLOG(2) << "block density: " << options.block_density;
 
     std::unique_ptr<BlockSparseMatrix> random_matrix(
-        BlockSparseMatrix::CreateRandomMatrix(options));
+        BlockSparseMatrix::CreateRandomMatrix(options, prng));
 
     const std::vector<CompressedRow>& row_blocks =
         random_matrix->block_structure()->rows;
diff --git a/internal/ceres/normal_prior_test.cc b/internal/ceres/normal_prior_test.cc
index 25e6322..881038f 100644
--- a/internal/ceres/normal_prior_test.cc
+++ b/internal/ceres/normal_prior_test.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -42,6 +42,7 @@
 TEST(NormalPriorTest, ResidualAtRandomPosition) {
   std::mt19937 prng;
   std::uniform_real_distribution<double> distribution(-1.0, 1.0);
+  auto randu = [&distribution, &prng] { return distribution(prng); };
   for (int num_rows = 1; num_rows < 5; ++num_rows) {
     for (int num_cols = 1; num_cols < 5; ++num_cols) {
       Vector b(num_cols);
@@ -50,7 +51,7 @@
       A.setRandom();
 
       auto* x = new double[num_cols];
-      for (int i = 0; i < num_cols; ++i) x[i] = distribution(prng);
+      std::generate_n(x, num_cols, randu);
 
       auto* jacobian = new double[num_rows * num_cols];
       Vector residuals(num_rows);
@@ -77,6 +78,7 @@
 TEST(NormalPriorTest, ResidualAtRandomPositionNullJacobians) {
   std::mt19937 prng;
   std::uniform_real_distribution<double> distribution(-1.0, 1.0);
+  auto randu = [&distribution, &prng] { return distribution(prng); };
   for (int num_rows = 1; num_rows < 5; ++num_rows) {
     for (int num_cols = 1; num_cols < 5; ++num_cols) {
       Vector b(num_cols);
@@ -85,7 +87,7 @@
       A.setRandom();
 
       auto* x = new double[num_cols];
-      for (int i = 0; i < num_cols; ++i) x[i] = distribution(prng);
+      std::generate_n(x, num_cols, randu);
 
       double* jacobians[1];
       jacobians[0] = nullptr;
diff --git a/internal/ceres/numeric_diff_cost_function_test.cc b/internal/ceres/numeric_diff_cost_function_test.cc
index b9a0beb..dcffcff 100644
--- a/internal/ceres/numeric_diff_cost_function_test.cc
+++ b/internal/ceres/numeric_diff_cost_function_test.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -35,6 +35,7 @@
 #include <array>
 #include <cmath>
 #include <memory>
+#include <random>
 #include <string>
 #include <vector>
 
@@ -310,6 +311,7 @@
 }
 
 TEST(NumericDiffCostFunction, RandomizedFunctorRidders) {
+  std::mt19937 prng;
   NumericDiffOptions options;
   // Larger initial step size is chosen to produce robust results in the
   // presence of random noise.
@@ -321,15 +323,16 @@
                                                1,  // number of residuals
                                                1   // size of x1
                                                >>(
-          new RandomizedFunctor(kNoiseFactor, kRandomSeed),
+          new RandomizedFunctor(kNoiseFactor, prng),
           TAKE_OWNERSHIP,
           1,
           options);
-  RandomizedFunctor functor(kNoiseFactor, kRandomSeed);
+  RandomizedFunctor functor(kNoiseFactor, prng);
   functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function);
 }
 
 TEST(NumericDiffCostFunction, RandomizedCostFunctionRidders) {
+  std::mt19937 prng;
   NumericDiffOptions options;
   // Larger initial step size is chosen to produce robust results in the
   // presence of random noise.
@@ -341,12 +344,12 @@
                                                1,  // number of residuals
                                                1   // size of x1
                                                >>(
-          new RandomizedCostFunction(kNoiseFactor, kRandomSeed),
+          new RandomizedCostFunction(kNoiseFactor, prng),
           TAKE_OWNERSHIP,
           1,
           options);
 
-  RandomizedFunctor functor(kNoiseFactor, kRandomSeed);
+  RandomizedFunctor functor(kNoiseFactor, prng);
   functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function);
 }
 
diff --git a/internal/ceres/numeric_diff_test_utils.cc b/internal/ceres/numeric_diff_test_utils.cc
index 3d78e36..9da514e 100644
--- a/internal/ceres/numeric_diff_test_utils.cc
+++ b/internal/ceres/numeric_diff_test_utils.cc
@@ -226,14 +226,7 @@
 }
 
 bool RandomizedFunctor::operator()(const double* x1, double* residuals) const {
-  double random_value =
-      static_cast<double>(rand()) / static_cast<double>(RAND_MAX);
-
-  // Normalize noise to [-factor, factor].
-  random_value *= 2.0;
-  random_value -= 1.0;
-  random_value *= noise_factor_;
-
+  double random_value = uniform_distribution_(*prng_);
   residuals[0] = x1[0] * x1[0] + random_value;
   return true;
 }
@@ -244,9 +237,6 @@
 
   const double kTolerance = 2e-4;
 
-  // Initialize random number generator with given seed.
-  srand(random_seed_);
-
   for (double& test : kTests) {
     double* parameters[] = {&test};
     double dydx;
diff --git a/internal/ceres/numeric_diff_test_utils.h b/internal/ceres/numeric_diff_test_utils.h
index 68ae5ed..8d3f14b 100644
--- a/internal/ceres/numeric_diff_test_utils.h
+++ b/internal/ceres/numeric_diff_test_utils.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -31,6 +31,8 @@
 #ifndef CERES_INTERNAL_NUMERIC_DIFF_TEST_UTILS_H_
 #define CERES_INTERNAL_NUMERIC_DIFF_TEST_UTILS_H_
 
+#include <random>
+
 #include "ceres/cost_function.h"
 #include "ceres/internal/export.h"
 #include "ceres/sized_cost_function.h"
@@ -116,8 +118,10 @@
 // y = x^2 + [random noise], dy/dx ~ 2x
 class CERES_NO_EXPORT RandomizedFunctor {
  public:
-  RandomizedFunctor(double noise_factor, unsigned int random_seed)
-      : noise_factor_(noise_factor), random_seed_(random_seed) {}
+  RandomizedFunctor(double noise_factor, std::mt19937& prng)
+      : noise_factor_(noise_factor),
+        prng_(&prng),
+        uniform_distribution_{-noise_factor, noise_factor} {}
 
   bool operator()(const double* x1, double* residuals) const;
   void ExpectCostFunctionEvaluationIsNearlyCorrect(
@@ -125,13 +129,16 @@
 
  private:
   double noise_factor_;
-  unsigned int random_seed_;
+  // Store the generator as a pointer to be able to modify the instance the
+  // pointer is pointing to.
+  std::mt19937* prng_;
+  mutable std::uniform_real_distribution<> uniform_distribution_;
 };
 
 class CERES_NO_EXPORT RandomizedCostFunction : public SizedCostFunction<1, 1> {
  public:
-  RandomizedCostFunction(double noise_factor, unsigned int random_seed)
-      : functor_(noise_factor, random_seed) {}
+  RandomizedCostFunction(double noise_factor, std::mt19937& prng)
+      : functor_(noise_factor, prng) {}
 
   bool Evaluate(double const* const* parameters,
                 double* residuals,
diff --git a/internal/ceres/partitioned_matrix_view_test.cc b/internal/ceres/partitioned_matrix_view_test.cc
index cb3dd14..e4484dd 100644
--- a/internal/ceres/partitioned_matrix_view_test.cc
+++ b/internal/ceres/partitioned_matrix_view_test.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
diff --git a/internal/ceres/reorder_program_test.cc b/internal/ceres/reorder_program_test.cc
index 872f381..d6a6577 100644
--- a/internal/ceres/reorder_program_test.cc
+++ b/internal/ceres/reorder_program_test.cc
@@ -265,7 +265,7 @@
   problem.GetResidualBlocks(&residual_block_ids);
   std::vector<ResidualBlock*> residual_blocks =
       problem.program().residual_blocks();
-  auto rng = std::default_random_engine{};
+  auto rng = std::mt19937{};
   for (int i = 1; i < 6; ++i) {
     std::shuffle(
         std::begin(residual_block_ids), std::end(residual_block_ids), rng);
diff --git a/internal/ceres/rotation_test.cc b/internal/ceres/rotation_test.cc
index d6de650..97ea514 100644
--- a/internal/ceres/rotation_test.cc
+++ b/internal/ceres/rotation_test.cc
@@ -32,6 +32,7 @@
 
 #include <cmath>
 #include <limits>
+#include <random>
 #include <string>
 
 #include "ceres/internal/eigen.h"
@@ -56,11 +57,6 @@
 const double kPi = 3.14159265358979323846;
 const double kHalfSqrt2 = 0.707106781186547524401;
 
-static double RandDouble() {
-  double r = rand();
-  return r / RAND_MAX;
-}
-
 // A tolerance value for floating-point comparisons.
 static double const kTolerance = numeric_limits<double>::epsilon() * 10;
 
@@ -322,20 +318,22 @@
 // Takes a bunch of random axis/angle values, converts them to quaternions,
 // and back again.
 TEST(Rotation, AngleAxisToQuaterionAndBack) {
-  srand(5);
+  std::mt19937 prng;
+  std::uniform_real_distribution<double> uniform_distribution{-1.0, 1.0};
   for (int i = 0; i < kNumTrials; i++) {
     double axis_angle[3];
     // Make an axis by choosing three random numbers in [-1, 1) and
     // normalizing.
     double norm = 0;
     for (double& coeff : axis_angle) {
-      coeff = RandDouble() * 2 - 1;
+      coeff = uniform_distribution(prng);
       norm += coeff * coeff;
     }
     norm = sqrt(norm);
 
     // Angle in [-pi, pi).
-    double theta = kPi * 2 * RandDouble() - kPi;
+    double theta = uniform_distribution(
+        prng, std::uniform_real_distribution<double>::param_type{-kPi, kPi});
     for (double& coeff : axis_angle) {
       coeff = coeff * theta / norm;
     }
@@ -355,13 +353,14 @@
 // Takes a bunch of random quaternions, converts them to axis/angle,
 // and back again.
 TEST(Rotation, QuaterionToAngleAxisAndBack) {
-  srand(5);
+  std::mt19937 prng;
+  std::uniform_real_distribution<double> uniform_distribution{-1.0, 1.0};
   for (int i = 0; i < kNumTrials; i++) {
     double quaternion[4];
     // Choose four random numbers in [-1, 1) and normalize.
     double norm = 0;
     for (double& coeff : quaternion) {
-      coeff = RandDouble() * 2 - 1;
+      coeff = uniform_distribution(prng);
       norm += coeff * coeff;
     }
     norm = sqrt(norm);
@@ -432,20 +431,24 @@
   double matrix[9];
   double out_axis_angle[3];
 
-  srand(5);
+  std::mt19937 prng;
+  std::uniform_real_distribution<double> uniform_distribution{-1.0, 1.0};
   for (int i = 0; i < kNumTrials; i++) {
     // Make an axis by choosing three random numbers in [-1, 1) and
     // normalizing.
     double norm = 0;
     for (double& coeff : in_axis_angle) {
-      coeff = RandDouble() * 2 - 1;
+      coeff = uniform_distribution(prng);
       norm += coeff * coeff;
     }
     norm = sqrt(norm);
 
     // Angle in [pi - kMaxSmallAngle, pi).
-    const double kMaxSmallAngle = 1e-8;
-    double theta = kPi - kMaxSmallAngle * RandDouble();
+    constexpr double kMaxSmallAngle = 1e-8;
+    double theta =
+        uniform_distribution(prng,
+                             std::uniform_real_distribution<double>::param_type{
+                                 kPi - kMaxSmallAngle, kPi});
 
     for (double& coeff : in_axis_angle) {
       coeff *= (theta / norm);
@@ -531,20 +534,22 @@
 // Takes a bunch of random axis/angle values, converts them to rotation
 // matrices, and back again.
 TEST(Rotation, AngleAxisToRotationMatrixAndBack) {
-  srand(5);
+  std::mt19937 prng;
+  std::uniform_real_distribution<double> uniform_distribution{-1.0, 1.0};
   for (int i = 0; i < kNumTrials; i++) {
     double axis_angle[3];
     // Make an axis by choosing three random numbers in [-1, 1) and
     // normalizing.
     double norm = 0;
     for (double& i : axis_angle) {
-      i = RandDouble() * 2 - 1;
+      i = uniform_distribution(prng);
       norm += i * i;
     }
     norm = sqrt(norm);
 
     // Angle in [-pi, pi).
-    double theta = kPi * 2 * RandDouble() - kPi;
+    double theta = uniform_distribution(
+        prng, std::uniform_real_distribution<double>::param_type{-kPi, kPi});
     for (double& i : axis_angle) {
       i = i * theta / norm;
     }
@@ -564,20 +569,25 @@
 // Takes a bunch of random axis/angle values near zero, converts them
 // to rotation matrices, and back again.
 TEST(Rotation, AngleAxisToRotationMatrixAndBackNearZero) {
-  srand(5);
+  std::mt19937 prng;
+  std::uniform_real_distribution<double> uniform_distribution{-1.0, 1.0};
   for (int i = 0; i < kNumTrials; i++) {
     double axis_angle[3];
     // Make an axis by choosing three random numbers in [-1, 1) and
     // normalizing.
     double norm = 0;
     for (double& i : axis_angle) {
-      i = RandDouble() * 2 - 1;
+      i = uniform_distribution(prng);
       norm += i * i;
     }
     norm = sqrt(norm);
 
     // Tiny theta.
-    double theta = 1e-16 * (kPi * 2 * RandDouble() - kPi);
+    constexpr double kScale = 1e-16;
+    double theta =
+        uniform_distribution(prng,
+                             std::uniform_real_distribution<double>::param_type{
+                                 -kScale * kPi, kScale * kPi});
     for (double& i : axis_angle) {
       i = i * theta / norm;
     }
@@ -647,11 +657,12 @@
 // Test that a random rotation produces an orthonormal rotation
 // matrix.
 TEST(EulerAnglesToRotationMatrix, IsOrthonormal) {
-  srand(5);
+  std::mt19937 prng;
+  std::uniform_real_distribution<double> uniform_distribution{-180.0, 180.0};
   for (int trial = 0; trial < kNumTrials; ++trial) {
     double euler_angles_degrees[3];
     for (double& euler_angles_degree : euler_angles_degrees) {
-      euler_angles_degree = RandDouble() * 360.0 - 180.0;
+      euler_angles_degree = uniform_distribution(prng);
     }
     double rotation_matrix[9];
     EulerAnglesToRotationMatrix(euler_angles_degrees, 3, rotation_matrix);
@@ -906,13 +917,15 @@
 
 // Verify that (a * b) * c == a * (b * c).
 TEST(Quaternion, MultiplicationIsAssociative) {
+  std::mt19937 prng;
+  std::uniform_real_distribution<double> uniform_distribution{-1.0, 1.0};
   double a[4];
   double b[4];
   double c[4];
   for (int i = 0; i < 4; ++i) {
-    a[i] = 2 * RandDouble() - 1;
-    b[i] = 2 * RandDouble() - 1;
-    c[i] = 2 * RandDouble() - 1;
+    a[i] = uniform_distribution(prng);
+    b[i] = uniform_distribution(prng);
+    c[i] = uniform_distribution(prng);
   }
 
   double ab[4];
@@ -932,6 +945,8 @@
 }
 
 TEST(AngleAxis, RotatePointGivesSameAnswerAsRotationMatrix) {
+  std::mt19937 prng;
+  std::uniform_real_distribution<double> uniform_distribution{-1.0, 1.0};
   double angle_axis[3];
   double R[9];
   double p[3];
@@ -943,8 +958,8 @@
     for (int j = 0; j < 50; ++j) {
       double norm2 = 0.0;
       for (int k = 0; k < 3; ++k) {
-        angle_axis[k] = 2.0 * RandDouble() - 1.0;
-        p[k] = 2.0 * RandDouble() - 1.0;
+        angle_axis[k] = uniform_distribution(prng);
+        p[k] = uniform_distribution(prng);
         norm2 = angle_axis[k] * angle_axis[k];
       }
 
@@ -976,6 +991,8 @@
 }
 
 TEST(AngleAxis, NearZeroRotatePointGivesSameAnswerAsRotationMatrix) {
+  std::mt19937 prng;
+  std::uniform_real_distribution<double> uniform_distribution{-1.0, 1.0};
   double angle_axis[3];
   double R[9];
   double p[3];
@@ -985,8 +1002,8 @@
   for (int i = 0; i < 10000; ++i) {
     double norm2 = 0.0;
     for (int k = 0; k < 3; ++k) {
-      angle_axis[k] = 2.0 * RandDouble() - 1.0;
-      p[k] = 2.0 * RandDouble() - 1.0;
+      angle_axis[k] = uniform_distribution(prng);
+      p[k] = uniform_distribution(prng);
       norm2 = angle_axis[k] * angle_axis[k];
     }
 
@@ -1111,7 +1128,12 @@
 }
 
 TEST(RotationMatrixToAngleAxis, ExhaustiveRoundTrip) {
-  const double kMaxSmallAngle = 1e-8;
+  constexpr double kMaxSmallAngle = 1e-8;
+  std::mt19937 prng;
+  std::uniform_real_distribution<double> uniform_distribution1{
+      kPi - kMaxSmallAngle, kPi};
+  std::uniform_real_distribution<double> uniform_distribution2{
+      -1.0, 2.0 * kMaxSmallAngle - 1.0};
   const int kNumSteps = 1000;
   for (int i = 0; i < kNumSteps; ++i) {
     const double theta = static_cast<double>(i) / kNumSteps * 2.0 * kPi;
@@ -1121,10 +1143,10 @@
       CheckRotationMatrixToAngleAxisRoundTrip(theta, phi, kPi);
       // Rotation of angle approximately Pi.
       CheckRotationMatrixToAngleAxisRoundTrip(
-          theta, phi, kPi - kMaxSmallAngle * RandDouble());
+          theta, phi, uniform_distribution1(prng));
       // Rotations of angle approximately zero.
       CheckRotationMatrixToAngleAxisRoundTrip(
-          theta, phi, kMaxSmallAngle * 2.0 * RandDouble() - 1.0);
+          theta, phi, uniform_distribution2(prng));
     }
   }
 }
diff --git a/internal/ceres/schur_eliminator_benchmark.cc b/internal/ceres/schur_eliminator_benchmark.cc
index 0973e73..c6425ce 100644
--- a/internal/ceres/schur_eliminator_benchmark.cc
+++ b/internal/ceres/schur_eliminator_benchmark.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2019 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -91,11 +91,9 @@
 
     matrix_ = std::make_unique<BlockSparseMatrix>(bs);
     double* values = matrix_->mutable_values();
-    std::mt19937 prng;
-    std::normal_distribution<double> standard_normal;
-    for (int i = 0; i < matrix_->num_nonzeros(); ++i) {
-      values[i] = standard_normal(prng);
-    }
+    std::generate_n(values, matrix_->num_nonzeros(), [this] {
+      return standard_normal(prng);
+    });
 
     b_.resize(matrix_->num_rows());
     b_.setRandom();
@@ -128,6 +126,8 @@
   Vector diagonal_;
   Vector z_;
   Vector y_;
+  std::mt19937 prng;
+  std::normal_distribution<> standard_normal;
 };
 
 static void BM_SchurEliminatorEliminate(benchmark::State& state) {
diff --git a/internal/ceres/schur_eliminator_test.cc b/internal/ceres/schur_eliminator_test.cc
index da27019..17ce67d 100644
--- a/internal/ceres/schur_eliminator_test.cc
+++ b/internal/ceres/schur_eliminator_test.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2019 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -284,10 +284,10 @@
   BlockSparseMatrix matrix(bs);
   double* values = matrix.mutable_values();
   std::mt19937 prng;
-  std::normal_distribution<double> standard_normal;
-  for (int i = 0; i < matrix.num_nonzeros(); ++i) {
-    values[i] = standard_normal(prng);
-  }
+  std::normal_distribution<> standard_normal;
+  std::generate_n(values, matrix.num_nonzeros(), [&prng, &standard_normal] {
+    return standard_normal(prng);
+  });
 
   Vector b(matrix.num_rows());
   b.setRandom();
diff --git a/internal/ceres/sparse_cholesky_test.cc b/internal/ceres/sparse_cholesky_test.cc
index 73e1455..2397fdd 100644
--- a/internal/ceres/sparse_cholesky_test.cc
+++ b/internal/ceres/sparse_cholesky_test.cc
@@ -55,7 +55,8 @@
     const int num_col_blocks,
     const int min_col_block_size,
     const int max_col_block_size,
-    const double block_density) {
+    const double block_density,
+    std::mt19937& prng) {
   // Create a random matrix
   BlockSparseMatrix::RandomMatrixOptions options;
   options.num_col_blocks = num_col_blocks;
@@ -66,7 +67,7 @@
   options.min_row_block_size = 1;
   options.max_row_block_size = max_col_block_size;
   options.block_density = block_density;
-  auto random_matrix = BlockSparseMatrix::CreateRandomMatrix(options);
+  auto random_matrix = BlockSparseMatrix::CreateRandomMatrix(options, prng);
 
   // Add a diagonal block sparse matrix to make it full rank.
   Vector diagonal = Vector::Ones(random_matrix->num_cols());
@@ -110,7 +111,8 @@
     const int num_blocks,
     const int min_block_size,
     const int max_block_size,
-    const double block_density) {
+    const double block_density,
+    std::mt19937& prng) {
   LinearSolver::Options sparse_cholesky_options;
   sparse_cholesky_options.sparse_linear_algebra_library_type =
       sparse_linear_algebra_library_type;
@@ -120,7 +122,7 @@
       sparse_cholesky->StorageType();
 
   auto m = CreateRandomFullRankMatrix(
-      num_blocks, min_block_size, max_block_size, block_density);
+      num_blocks, min_block_size, max_block_size, block_density, prng);
   auto inner_product_computer = InnerProductComputer::Create(*m, storage_type);
   inner_product_computer->Compute();
   CompressedRowSparseMatrix* lhs = inner_product_computer->mutable_result();
@@ -184,7 +186,8 @@
                                    num_blocks,
                                    kMinBlockSize,
                                    kMaxBlockSize,
-                                   block_density);
+                                   block_density,
+                                   prng);
     }
   }
 }
diff --git a/internal/ceres/subset_preconditioner_test.cc b/internal/ceres/subset_preconditioner_test.cc
index cd86695..5fb7f26 100644
--- a/internal/ceres/subset_preconditioner_test.cc
+++ b/internal/ceres/subset_preconditioner_test.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2017 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -31,6 +31,7 @@
 #include "ceres/subset_preconditioner.h"
 
 #include <memory>
+#include <random>
 
 #include "Eigen/Dense"
 #include "Eigen/SparseCore"
@@ -100,13 +101,13 @@
     options.max_row_block_size = 4;
     options.block_density = 0.9;
 
-    m_ = BlockSparseMatrix::CreateRandomMatrix(options);
+    m_ = BlockSparseMatrix::CreateRandomMatrix(options, prng);
     start_row_block_ = m_->block_structure()->rows.size();
 
     // Ensure that the bottom part of the matrix has the same column
     // block structure.
     options.col_blocks = m_->block_structure()->cols;
-    b_ = BlockSparseMatrix::CreateRandomMatrix(options);
+    b_ = BlockSparseMatrix::CreateRandomMatrix(options, prng);
     m_->AppendRows(*b_);
 
     // Create a Identity block diagonal matrix with the same column
@@ -132,6 +133,7 @@
   std::unique_ptr<Preconditioner> preconditioner_;
   Vector diagonal_;
   int start_row_block_;
+  std::mt19937 prng;
 };
 
 TEST_P(SubsetPreconditionerTest, foo) {
diff --git a/internal/ceres/triplet_sparse_matrix.cc b/internal/ceres/triplet_sparse_matrix.cc
index 49d367a..37a779f 100644
--- a/internal/ceres/triplet_sparse_matrix.cc
+++ b/internal/ceres/triplet_sparse_matrix.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -310,7 +310,8 @@
 }
 
 std::unique_ptr<TripletSparseMatrix> TripletSparseMatrix::CreateRandomMatrix(
-    const TripletSparseMatrix::RandomMatrixOptions& options) {
+    const TripletSparseMatrix::RandomMatrixOptions& options,
+    std::mt19937& prng) {
   CHECK_GT(options.num_rows, 0);
   CHECK_GT(options.num_cols, 0);
   CHECK_GT(options.density, 0.0);
@@ -319,7 +320,6 @@
   std::vector<int> rows;
   std::vector<int> cols;
   std::vector<double> values;
-  std::mt19937 prng;
   std::uniform_real_distribution<double> uniform01(0.0, 1.0);
   std::normal_distribution<double> standard_normal;
   while (rows.empty()) {
diff --git a/internal/ceres/triplet_sparse_matrix.h b/internal/ceres/triplet_sparse_matrix.h
index c962455..549966a 100644
--- a/internal/ceres/triplet_sparse_matrix.h
+++ b/internal/ceres/triplet_sparse_matrix.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -32,6 +32,7 @@
 #define CERES_INTERNAL_TRIPLET_SPARSE_MATRIX_H_
 
 #include <memory>
+#include <random>
 #include <vector>
 
 #include "ceres/crs_matrix.h"
@@ -135,7 +136,8 @@
   // normally distributed and whose structure is determined by
   // RandomMatrixOptions.
   static std::unique_ptr<TripletSparseMatrix> CreateRandomMatrix(
-      const TripletSparseMatrix::RandomMatrixOptions& options);
+      const TripletSparseMatrix::RandomMatrixOptions& options,
+      std::mt19937& prng);
 
   // Load a triplet sparse matrix from a text file.
   static std::unique_ptr<TripletSparseMatrix> CreateFromTextFile(FILE* file);