diff --git a/examples/BUILD b/examples/BUILD
index 7bacf22..00b799e 100644
--- a/examples/BUILD
+++ b/examples/BUILD
@@ -47,7 +47,6 @@
         "bal_problem.cc",
         "bal_problem.h",
         "bundle_adjuster.cc",
-        "random.h",
         "snavely_reprojection_error.h",
     ],
     copts = EXAMPLE_COPTS,
@@ -69,7 +68,6 @@
 cc_binary(
     name = "robot_pose_mle",
     srcs = [
-        "random.h",
         "robot_pose_mle.cc",
     ],
     copts = EXAMPLE_COPTS,
diff --git a/examples/bal_problem.cc b/examples/bal_problem.cc
index 99bc075..2423fb6 100644
--- a/examples/bal_problem.cc
+++ b/examples/bal_problem.cc
@@ -32,15 +32,15 @@
 
 #include <algorithm>
 #include <cstdio>
-#include <cstdlib>
 #include <fstream>
+#include <functional>
+#include <random>
 #include <string>
 #include <vector>
 
 #include "Eigen/Core"
 #include "ceres/rotation.h"
 #include "glog/logging.h"
-#include "random.h"
 
 namespace ceres::examples {
 namespace {
@@ -55,9 +55,9 @@
   }
 }
 
-void PerturbPoint3(const double sigma, double* point) {
+void PerturbPoint3(std::function<double()> dist, double* point) {
   for (int i = 0; i < 3; ++i) {
-    point[i] += RandNormal() * sigma;
+    point[i] += dist();
   }
 }
 
@@ -297,14 +297,19 @@
   CHECK_GE(point_sigma, 0.0);
   CHECK_GE(rotation_sigma, 0.0);
   CHECK_GE(translation_sigma, 0.0);
-
+  std::mt19937 prng;
+  std::normal_distribution point_noise_distribution(0.0, point_sigma);
   double* points = mutable_points();
   if (point_sigma > 0) {
     for (int i = 0; i < num_points_; ++i) {
-      PerturbPoint3(point_sigma, points + 3 * i);
+      PerturbPoint3(std::bind(point_noise_distribution, std::ref(prng)),
+                    points + 3 * i);
     }
   }
 
+  std::normal_distribution rotation_noise_distribution(0.0, point_sigma);
+  std::normal_distribution translation_noise_distribution(0.0,
+                                                          translation_sigma);
   for (int i = 0; i < num_cameras_; ++i) {
     double* camera = mutable_cameras() + camera_block_size() * i;
 
@@ -314,12 +319,14 @@
     // representation.
     CameraToAngleAxisAndCenter(camera, angle_axis, center);
     if (rotation_sigma > 0.0) {
-      PerturbPoint3(rotation_sigma, angle_axis);
+      PerturbPoint3(std::bind(rotation_noise_distribution, std::ref(prng)),
+                    angle_axis);
     }
     AngleAxisAndCenterToCamera(angle_axis, center, camera);
 
     if (translation_sigma > 0.0) {
-      PerturbPoint3(translation_sigma, camera + camera_block_size() - 6);
+      PerturbPoint3(std::bind(translation_noise_distribution, std::ref(prng)),
+                    camera + camera_block_size() - 6);
     }
   }
 }
diff --git a/examples/robot_pose_mle.cc b/examples/robot_pose_mle.cc
index c199b89..f924bc3 100644
--- a/examples/robot_pose_mle.cc
+++ b/examples/robot_pose_mle.cc
@@ -127,13 +127,13 @@
 
 #include <cmath>
 #include <cstdio>
+#include <random>
 #include <vector>
 
 #include "ceres/ceres.h"
 #include "ceres/dynamic_autodiff_cost_function.h"
 #include "gflags/gflags.h"
 #include "glog/logging.h"
-#include "random.h"
 
 using ceres::AutoDiffCostFunction;
 using ceres::CauchyLoss;
@@ -143,7 +143,6 @@
 using ceres::Problem;
 using ceres::Solve;
 using ceres::Solver;
-using ceres::examples::RandNormal;
 using std::min;
 using std::vector;
 
@@ -247,6 +246,10 @@
   const int num_steps =
       static_cast<int>(ceil(CERES_GET_FLAG(FLAGS_corridor_length) /
                             CERES_GET_FLAG(FLAGS_pose_separation)));
+  std::mt19937 prng;
+  std::normal_distribution odometry_noise(
+      0.0, CERES_GET_FLAG(FLAGS_odometry_stddev));
+  std::normal_distribution range_noise(0.0, CERES_GET_FLAG(FLAGS_range_stddev));
 
   // The robot starts out at the origin.
   double robot_location = 0.0;
@@ -258,10 +261,8 @@
     const double actual_range =
         CERES_GET_FLAG(FLAGS_corridor_length) - robot_location;
     const double observed_odometry =
-        RandNormal() * CERES_GET_FLAG(FLAGS_odometry_stddev) +
-        actual_odometry_value;
-    const double observed_range =
-        RandNormal() * CERES_GET_FLAG(FLAGS_range_stddev) + actual_range;
+        actual_odometry_value + odometry_noise(prng);
+    const double observed_range = actual_range + range_noise(prng);
     odometry_values->push_back(observed_odometry);
     range_readings->push_back(observed_range);
   }
diff --git a/internal/ceres/autodiff_test.cc b/internal/ceres/autodiff_test.cc
index 5f769a0..629eea8 100644
--- a/internal/ceres/autodiff_test.cc
+++ b/internal/ceres/autodiff_test.cc
@@ -30,7 +30,8 @@
 
 #include "ceres/internal/autodiff.h"
 
-#include "ceres/random.h"
+#include <random>
+
 #include "gtest/gtest.h"
 
 namespace ceres::internal {
@@ -162,17 +163,18 @@
 
 // Test projective camera model projector.
 TEST(AutoDiff, ProjectiveCameraModel) {
-  srand(5);
   double const tol = 1e-10;  // floating-point tolerance.
   double const del = 1e-4;   // finite-difference step.
   double const err = 1e-6;   // finite-difference tolerance.
 
   Projective b;
+  std::mt19937 prng;
+  std::uniform_real_distribution<double> uniform01(0.0, 1.0);
 
   // Make random P and X, in a single vector.
   double PX[12 + 4];
   for (double& PX_i : PX) {
-    PX_i = RandDouble();
+    PX_i = uniform01(prng);
   }
 
   // Handy names for the P and X parts.
@@ -282,16 +284,17 @@
 
 // This test is similar in structure to the previous one.
 TEST(AutoDiff, Metric) {
-  srand(5);
   double const tol = 1e-10;  // floating-point tolerance.
   double const del = 1e-4;   // finite-difference step.
-  double const err = 1e-5;   // finite-difference tolerance.
+  double const err = 2e-5;   // finite-difference tolerance.
 
   Metric b;
 
   // Make random parameter vector.
   double qcX[4 + 3 + 3];
-  for (double& qcX_i : qcX) qcX_i = RandDouble();
+  std::mt19937 prng;
+  std::uniform_real_distribution<double> uniform01(0.0, 1.0);
+  for (double& qcX_i : qcX) qcX_i = uniform01(prng);
 
   // Handy names.
   double* q = qcX;
diff --git a/internal/ceres/block_sparse_matrix.cc b/internal/ceres/block_sparse_matrix.cc
index 22224a3..9ae4384 100644
--- a/internal/ceres/block_sparse_matrix.cc
+++ b/internal/ceres/block_sparse_matrix.cc
@@ -33,12 +33,12 @@
 #include <algorithm>
 #include <cstddef>
 #include <memory>
+#include <random>
 #include <vector>
 
 #include "ceres/block_structure.h"
 #include "ceres/crs_matrix.h"
 #include "ceres/internal/eigen.h"
-#include "ceres/random.h"
 #include "ceres/small_blas.h"
 #include "ceres/triplet_sparse_matrix.h"
 #include "glog/logging.h"
@@ -382,6 +382,11 @@
   CHECK_GT(options.block_density, 0.0);
   CHECK_LE(options.block_density, 1.0);
 
+  std::mt19937 prng;
+  std::uniform_int_distribution col_distribution(options.min_col_block_size,
+                                                 options.max_col_block_size);
+  std::uniform_int_distribution row_distribution(options.min_row_block_size,
+                                                 options.max_row_block_size);
   auto* bs = new CompressedRowBlockStructure();
   if (options.col_blocks.empty()) {
     CHECK_GT(options.num_col_blocks, 0);
@@ -392,10 +397,7 @@
     // Generate the col block structure.
     int col_block_position = 0;
     for (int i = 0; i < options.num_col_blocks; ++i) {
-      // Generate a random integer in [min_col_block_size, max_col_block_size]
-      const int delta_block_size =
-          Uniform(options.max_col_block_size - options.min_col_block_size);
-      const int col_block_size = options.min_col_block_size + delta_block_size;
+      const int col_block_size = col_distribution(prng);
       bs->cols.emplace_back(col_block_size, col_block_position);
       col_block_position += col_block_size;
     }
@@ -404,22 +406,21 @@
   }
 
   bool matrix_has_blocks = false;
+  std::uniform_real_distribution uniform01(0.0, 1.0);
   while (!matrix_has_blocks) {
     VLOG(1) << "Clearing";
     bs->rows.clear();
     int row_block_position = 0;
     int value_position = 0;
     for (int r = 0; r < options.num_row_blocks; ++r) {
-      const int delta_block_size =
-          Uniform(options.max_row_block_size - options.min_row_block_size);
-      const int row_block_size = options.min_row_block_size + delta_block_size;
+      const int row_block_size = row_distribution(prng);
       bs->rows.emplace_back();
       CompressedRow& row = bs->rows.back();
       row.block.size = row_block_size;
       row.block.position = row_block_position;
       row_block_position += row_block_size;
       for (int c = 0; c < bs->cols.size(); ++c) {
-        if (RandDouble() > options.block_density) continue;
+        if (uniform01(prng) > options.block_density) continue;
 
         row.cells.emplace_back();
         Cell& cell = row.cells.back();
@@ -433,8 +434,9 @@
 
   auto matrix = std::make_unique<BlockSparseMatrix>(bs);
   double* values = matrix->mutable_values();
+  std::normal_distribution standard_normal_distribution;
   for (int i = 0; i < matrix->num_nonzeros(); ++i) {
-    values[i] = RandNormal();
+    values[i] = standard_normal_distribution(prng);
   }
 
   return matrix;
diff --git a/internal/ceres/compressed_row_sparse_matrix.cc b/internal/ceres/compressed_row_sparse_matrix.cc
index d663f4a..bd94fdc 100644
--- a/internal/ceres/compressed_row_sparse_matrix.cc
+++ b/internal/ceres/compressed_row_sparse_matrix.cc
@@ -31,13 +31,14 @@
 #include "ceres/compressed_row_sparse_matrix.h"
 
 #include <algorithm>
+#include <functional>
 #include <memory>
 #include <numeric>
+#include <random>
 #include <vector>
 
 #include "ceres/crs_matrix.h"
 #include "ceres/internal/export.h"
-#include "ceres/random.h"
 #include "ceres/triplet_sparse_matrix.h"
 #include "glog/logging.h"
 
@@ -122,6 +123,7 @@
                     const int num_cols,
                     const int row_block_begin,
                     const int col_block_begin,
+                    std::function<double()> dist,
                     std::vector<int>* rows,
                     std::vector<int>* cols,
                     std::vector<double>* values) {
@@ -129,19 +131,20 @@
     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(RandNormal());
+      values->push_back(dist());
     }
   }
 }
 
 void AddSymmetricRandomBlock(const int num_rows,
                              const int row_block_begin,
+                             std::function<double()> dist,
                              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 = RandNormal();
+      const double v = dist();
       rows->push_back(row_block_begin + r);
       cols->push_back(row_block_begin + c);
       values->push_back(v);
@@ -644,21 +647,26 @@
   vector<int> row_blocks;
   vector<int> col_blocks;
 
+  std::mt19937 prng;
+  std::uniform_int_distribution col_distribution(options.min_col_block_size,
+                                                 options.max_col_block_size);
+  std::uniform_int_distribution row_distribution(options.min_row_block_size,
+                                                 options.max_row_block_size);
+  std::uniform_real_distribution uniform01(0.0, 1.0);
+  std::normal_distribution standard_normal_distribution;
+  auto values_dist = std::bind(standard_normal_distribution, std::ref(prng));
+
   // Generate the row block structure.
   for (int i = 0; i < options.num_row_blocks; ++i) {
     // Generate a random integer in [min_row_block_size, max_row_block_size]
-    const int delta_block_size =
-        Uniform(options.max_row_block_size - options.min_row_block_size);
-    row_blocks.push_back(options.min_row_block_size + delta_block_size);
+    row_blocks.push_back(row_distribution(prng));
   }
 
   if (options.storage_type == StorageType::UNSYMMETRIC) {
     // Generate the col block structure.
     for (int i = 0; i < options.num_col_blocks; ++i) {
       // Generate a random integer in [min_col_block_size, max_col_block_size]
-      const int delta_block_size =
-          Uniform(options.max_col_block_size - options.min_col_block_size);
-      col_blocks.push_back(options.min_col_block_size + delta_block_size);
+      col_blocks.push_back(col_distribution(prng));
     }
   } else {
     // Symmetric matrices (LOWER_TRIANGULAR or UPPER_TRIANGULAR);
@@ -695,7 +703,7 @@
         }
 
         // Randomly determine if this block is present or not.
-        if (RandDouble() <= options.block_density) {
+        if (uniform01(prng) <= options.block_density) {
           // If the matrix is symmetric, then we take care to generate
           // symmetric diagonal blocks.
           if (options.storage_type == StorageType::UNSYMMETRIC || r != c) {
@@ -703,12 +711,14 @@
                            col_blocks[c],
                            row_block_begin,
                            col_block_begin,
+                           values_dist,
                            &tsm_rows,
                            &tsm_cols,
                            &tsm_values);
           } else {
             AddSymmetricRandomBlock(row_blocks[r],
                                     row_block_begin,
+                                    values_dist,
                                     &tsm_rows,
                                     &tsm_cols,
                                     &tsm_values);
diff --git a/internal/ceres/compressed_row_sparse_matrix_test.cc b/internal/ceres/compressed_row_sparse_matrix_test.cc
index 4151845..db90d8b 100644
--- a/internal/ceres/compressed_row_sparse_matrix_test.cc
+++ b/internal/ceres/compressed_row_sparse_matrix_test.cc
@@ -33,6 +33,7 @@
 #include <algorithm>
 #include <memory>
 #include <numeric>
+#include <random>
 #include <string>
 
 #include "Eigen/SparseCore"
@@ -40,7 +41,6 @@
 #include "ceres/crs_matrix.h"
 #include "ceres/internal/eigen.h"
 #include "ceres/linear_least_squares_problems.h"
-#include "ceres/random.h"
 #include "ceres/triplet_sparse_matrix.h"
 #include "glog/logging.h"
 #include "gtest/gtest.h"
@@ -404,7 +404,8 @@
   const int kMinBlockSize = 1;
   const int kMaxBlockSize = 5;
   const int kNumTrials = 10;
-
+  std::mt19937 prng;
+  std::uniform_real_distribution uniform(0.5, 1.0);
   for (int num_blocks = kMinNumBlocks; num_blocks < kMaxNumBlocks;
        ++num_blocks) {
     for (int trial = 0; trial < kNumTrials; ++trial) {
@@ -416,7 +417,7 @@
       options.num_row_blocks = 2 * num_blocks;
       options.min_row_block_size = kMinBlockSize;
       options.max_row_block_size = kMaxBlockSize;
-      options.block_density = std::max(0.5, RandDouble());
+      options.block_density = uniform(prng);
       options.storage_type = ::testing::get<0>(param);
       std::unique_ptr<CompressedRowSparseMatrix> matrix =
           CompressedRowSparseMatrix::CreateRandomMatrix(options);
@@ -473,7 +474,8 @@
   const int kMinBlockSize = 1;
   const int kMaxBlockSize = 5;
   const int kNumTrials = 10;
-
+  std::mt19937 prng;
+  std::uniform_real_distribution uniform(0.5, 1.0);
   for (int num_blocks = kMinNumBlocks; num_blocks < kMaxNumBlocks;
        ++num_blocks) {
     for (int trial = 0; trial < kNumTrials; ++trial) {
@@ -485,7 +487,7 @@
       options.num_row_blocks = 2 * num_blocks;
       options.min_row_block_size = kMinBlockSize;
       options.max_row_block_size = kMaxBlockSize;
-      options.block_density = std::max(0.5, RandDouble());
+      options.block_density = uniform(prng);
       options.storage_type = ::testing::get<0>(param);
       std::unique_ptr<CompressedRowSparseMatrix> matrix =
           CompressedRowSparseMatrix::CreateRandomMatrix(options);
@@ -542,7 +544,8 @@
   const int kMinBlockSize = 1;
   const int kMaxBlockSize = 5;
   const int kNumTrials = 10;
-
+  std::mt19937 prng;
+  std::uniform_real_distribution uniform(0.5, 1.0);
   for (int num_blocks = kMinNumBlocks; num_blocks < kMaxNumBlocks;
        ++num_blocks) {
     for (int trial = 0; trial < kNumTrials; ++trial) {
@@ -554,7 +557,7 @@
       options.num_row_blocks = 2 * num_blocks;
       options.min_row_block_size = kMinBlockSize;
       options.max_row_block_size = kMaxBlockSize;
-      options.block_density = std::max(0.5, RandDouble());
+      options.block_density = uniform(prng);
       options.storage_type = ::testing::get<0>(param);
       std::unique_ptr<CompressedRowSparseMatrix> matrix =
           CompressedRowSparseMatrix::CreateRandomMatrix(options);
diff --git a/internal/ceres/corrector_test.cc b/internal/ceres/corrector_test.cc
index d2f7c84..b1e9bee 100644
--- a/internal/ceres/corrector_test.cc
+++ b/internal/ceres/corrector_test.cc
@@ -34,9 +34,9 @@
 #include <cmath>
 #include <cstdlib>
 #include <cstring>
+#include <random>
 
 #include "ceres/internal/eigen.h"
-#include "ceres/random.h"
 #include "gtest/gtest.h"
 
 namespace ceres {
@@ -160,18 +160,19 @@
   // and hessians.
   Matrix c_hess(2, 2);
   Vector c_grad(2);
-
-  srand(5);
+  std::mt19937 prng;
+  std::uniform_real_distribution uniform01(0.0, 1.0);
   for (int iter = 0; iter < 10000; ++iter) {
     // Initialize the jacobian and residual.
-    for (double& jacobian_entry : jacobian) jacobian_entry = RandDouble();
-    for (double& residual : residuals) residual = RandDouble();
+    for (double& jacobian_entry : jacobian) jacobian_entry = uniform01(prng);
+    for (double& residual : residuals) residual = uniform01(prng);
 
     const double sq_norm = res.dot(res);
 
     rho[0] = sq_norm;
-    rho[1] = RandDouble();
-    rho[2] = 2.0 * RandDouble() - 1.0;
+    rho[1] = uniform01(prng);
+    rho[2] = uniform01(
+        prng, std::uniform_real_distribution<double>::param_type(-1, 1));
 
     // If rho[2] > 0, then the curvature correction to the correction
     // and the gauss newton approximation will match. Otherwise, we
@@ -227,10 +228,11 @@
   Matrix c_hess(2, 2);
   Vector c_grad(2);
 
-  srand(5);
+  std::mt19937 prng;
+  std::uniform_real_distribution uniform01(0.0, 1.0);
   for (int iter = 0; iter < 10000; ++iter) {
     // Initialize the jacobian.
-    for (double& jacobian_entry : jacobian) jacobian_entry = RandDouble();
+    for (double& jacobian_entry : jacobian) jacobian_entry = uniform01(prng);
 
     // Zero residuals
     res.setZero();
@@ -238,8 +240,9 @@
     const double sq_norm = res.dot(res);
 
     rho[0] = sq_norm;
-    rho[1] = RandDouble();
-    rho[2] = 2 * RandDouble() - 1.0;
+    rho[1] = uniform01(prng);
+    rho[2] = uniform01(
+        prng, std::uniform_real_distribution<double>::param_type(-1, 1));
 
     // Ground truth values.
     g_res = sqrt(rho[1]) * res;
diff --git a/internal/ceres/dense_linear_solver_benchmark.cc b/internal/ceres/dense_linear_solver_benchmark.cc
index 53628ad..a01e0d6 100644
--- a/internal/ceres/dense_linear_solver_benchmark.cc
+++ b/internal/ceres/dense_linear_solver_benchmark.cc
@@ -34,7 +34,6 @@
 #include "ceres/dense_sparse_matrix.h"
 #include "ceres/internal/config.h"
 #include "ceres/linear_solver.h"
-#include "ceres/random.h"
 
 namespace ceres::internal {
 
diff --git a/internal/ceres/gradient_checker_test.cc b/internal/ceres/gradient_checker_test.cc
index aedb0b9..5de4778 100644
--- a/internal/ceres/gradient_checker_test.cc
+++ b/internal/ceres/gradient_checker_test.cc
@@ -34,12 +34,12 @@
 
 #include <cmath>
 #include <cstdlib>
+#include <random>
 #include <utility>
 #include <vector>
 
 #include "ceres/cost_function.h"
 #include "ceres/problem.h"
-#include "ceres/random.h"
 #include "ceres/solver.h"
 #include "ceres/test_util.h"
 #include "glog/logging.h"
@@ -60,12 +60,15 @@
 class GoodTestTerm : public CostFunction {
  public:
   GoodTestTerm(int arity, int const* dim) : arity_(arity), return_value_(true) {
+    std::mt19937 prng;
+    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] = 2.0 * RandDouble() - 1.0;
+        a_[j][u] = distribution(prng);
       }
     }
 
@@ -119,12 +122,14 @@
 class BadTestTerm : public CostFunction {
  public:
   BadTestTerm(int arity, int const* dim) : arity_(arity) {
+    std::mt19937 prng;
+    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] = 2.0 * RandDouble() - 1.0;
+        a_[j][u] = distribution(prng);
       }
     }
 
@@ -194,8 +199,6 @@
 }
 
 TEST(GradientChecker, SmokeTest) {
-  srand(5);
-
   // Test with 3 blocks of size 2, 3 and 4.
   int const num_parameters = 3;
   std::vector<int> parameter_sizes(3);
@@ -205,10 +208,12 @@
 
   // Make a random set of blocks.
   FixedArray<double*> parameters(num_parameters);
+  std::mt19937 prng;
+  std::uniform_real_distribution distribution(-1.0, 1.0);
   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] = 2.0 * RandDouble() - 1.0;
+      parameters[j][u] = distribution(prng);
     }
   }
 
diff --git a/internal/ceres/gradient_checking_cost_function_test.cc b/internal/ceres/gradient_checking_cost_function_test.cc
index 5c6046f..99e0727 100644
--- a/internal/ceres/gradient_checking_cost_function_test.cc
+++ b/internal/ceres/gradient_checking_cost_function_test.cc
@@ -33,6 +33,7 @@
 #include <cmath>
 #include <cstdint>
 #include <memory>
+#include <random>
 #include <vector>
 
 #include "ceres/cost_function.h"
@@ -41,7 +42,6 @@
 #include "ceres/parameter_block.h"
 #include "ceres/problem_impl.h"
 #include "ceres/program.h"
-#include "ceres/random.h"
 #include "ceres/residual_block.h"
 #include "ceres/sized_cost_function.h"
 #include "ceres/types.h"
@@ -70,12 +70,14 @@
   // 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 uniform01(-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] = 2.0 * RandDouble() - 1.0;
+        a_[j][u] = uniform01(prng);
       }
     }
 
@@ -130,18 +132,18 @@
 };
 
 TEST(GradientCheckingCostFunction, ResidualsAndJacobiansArePreservedTest) {
-  srand(5);
-
   // Test with 3 blocks of size 2, 3 and 4.
   int const arity = 3;
   int const dim[arity] = {2, 3, 4};
 
   // Make a random set of blocks.
   vector<double*> parameters(arity);
+  std::mt19937 prng;
+  std::uniform_real_distribution uniform01(-1.0, 1.0);
   for (int j = 0; j < arity; ++j) {
     parameters[j] = new double[dim[j]];
     for (int u = 0; u < dim[j]; ++u) {
-      parameters[j][u] = 2.0 * RandDouble() - 1.0;
+      parameters[j][u] = uniform01(prng);
     }
   }
 
@@ -187,18 +189,18 @@
 }
 
 TEST(GradientCheckingCostFunction, SmokeTest) {
-  srand(5);
-
   // Test with 3 blocks of size 2, 3 and 4.
   int const arity = 3;
   int const dim[arity] = {2, 3, 4};
 
   // Make a random set of blocks.
   vector<double*> parameters(arity);
+  std::mt19937 prng;
+  std::uniform_real_distribution uniform01(-1.0, 1.0);
   for (int j = 0; j < arity; ++j) {
     parameters[j] = new double[dim[j]];
     for (int u = 0; u < dim[j]; ++u) {
-      parameters[j][u] = 2.0 * RandDouble() - 1.0;
+      parameters[j][u] = uniform01(prng);
     }
   }
 
diff --git a/internal/ceres/inner_product_computer_test.cc b/internal/ceres/inner_product_computer_test.cc
index a1f5c49..9c0382b 100644
--- a/internal/ceres/inner_product_computer_test.cc
+++ b/internal/ceres/inner_product_computer_test.cc
@@ -32,11 +32,11 @@
 
 #include <memory>
 #include <numeric>
+#include <random>
 
 #include "Eigen/SparseCore"
 #include "ceres/block_sparse_matrix.h"
 #include "ceres/internal/eigen.h"
-#include "ceres/random.h"
 #include "ceres/triplet_sparse_matrix.h"
 #include "glog/logging.h"
 #include "gtest/gtest.h"
@@ -77,11 +77,11 @@
   }
 
 TEST(InnerProductComputer, NormalOperation) {
-  // "Randomly generated seed."
-  SetRandomState(29823);
   const int kMaxNumRowBlocks = 10;
   const int kMaxNumColBlocks = 10;
   const int kNumTrials = 10;
+  std::mt19937 prng;
+  std::uniform_real_distribution distribution(0.01, 1.0);
 
   // Create a random matrix, compute its outer product using Eigen and
   // ComputeOuterProduct. Convert both matrices to dense matrices and
@@ -98,7 +98,8 @@
         options.max_row_block_size = 5;
         options.min_col_block_size = 1;
         options.max_col_block_size = 10;
-        options.block_density = std::max(0.1, RandDouble());
+        options.block_density = distribution(prng);
+
         VLOG(2) << "num row blocks: " << options.num_row_blocks;
         VLOG(2) << "num col blocks: " << options.num_col_blocks;
         VLOG(2) << "min row block size: " << options.min_row_block_size;
@@ -143,6 +144,8 @@
   const int kNumRowBlocks = 10;
   const int kNumColBlocks = 20;
   const int kNumTrials = 5;
+  std::mt19937 prng;
+  std::uniform_real_distribution distribution(0.01, 1.0);
 
   // Create a random matrix, compute its outer product using Eigen and
   // ComputeInnerProductComputer. Convert both matrices to dense matrices and
@@ -155,7 +158,7 @@
     options.max_row_block_size = 5;
     options.min_col_block_size = 1;
     options.max_col_block_size = 10;
-    options.block_density = std::min(0.01, RandDouble());
+    options.block_density = distribution(prng);
 
     VLOG(2) << "num row blocks: " << options.num_row_blocks;
     VLOG(2) << "num col blocks: " << options.num_col_blocks;
diff --git a/internal/ceres/normal_prior_test.cc b/internal/ceres/normal_prior_test.cc
index c66f9d6..9c84555 100644
--- a/internal/ceres/normal_prior_test.cc
+++ b/internal/ceres/normal_prior_test.cc
@@ -31,43 +31,26 @@
 #include "ceres/normal_prior.h"
 
 #include <cstddef>
+#include <random>
 
 #include "ceres/internal/eigen.h"
-#include "ceres/random.h"
 #include "gtest/gtest.h"
 
 namespace ceres {
 namespace internal {
 
-namespace {
-
-void RandomVector(Vector* v) {
-  for (int r = 0; r < v->rows(); ++r) (*v)[r] = 2 * RandDouble() - 1;
-}
-
-void RandomMatrix(Matrix* m) {
-  for (int r = 0; r < m->rows(); ++r) {
-    for (int c = 0; c < m->cols(); ++c) {
-      (*m)(r, c) = 2 * RandDouble() - 1;
-    }
-  }
-}
-
-}  // namespace
-
 TEST(NormalPriorTest, ResidualAtRandomPosition) {
-  srand(5);
-
+  std::mt19937 prng;
+  std::uniform_real_distribution distribution(-1.0, 1.0);
   for (int num_rows = 1; num_rows < 5; ++num_rows) {
     for (int num_cols = 1; num_cols < 5; ++num_cols) {
       Vector b(num_cols);
-      RandomVector(&b);
-
+      b.setRandom();
       Matrix A(num_rows, num_cols);
-      RandomMatrix(&A);
+      A.setRandom();
 
       auto* x = new double[num_cols];
-      for (int i = 0; i < num_cols; ++i) x[i] = 2 * RandDouble() - 1;
+      for (int i = 0; i < num_cols; ++i) x[i] = distribution(prng);
 
       auto* jacobian = new double[num_rows * num_cols];
       Vector residuals(num_rows);
@@ -92,18 +75,17 @@
 }
 
 TEST(NormalPriorTest, ResidualAtRandomPositionNullJacobians) {
-  srand(5);
-
+  std::mt19937 prng;
+  std::uniform_real_distribution distribution(-1.0, 1.0);
   for (int num_rows = 1; num_rows < 5; ++num_rows) {
     for (int num_cols = 1; num_cols < 5; ++num_cols) {
       Vector b(num_cols);
-      RandomVector(&b);
-
+      b.setRandom();
       Matrix A(num_rows, num_cols);
-      RandomMatrix(&A);
+      A.setRandom();
 
       auto* x = new double[num_cols];
-      for (int i = 0; i < num_cols; ++i) x[i] = 2 * RandDouble() - 1;
+      for (int i = 0; i < num_cols; ++i) x[i] = distribution(prng);
 
       double* jacobians[1];
       jacobians[0] = nullptr;
diff --git a/internal/ceres/partitioned_matrix_view_test.cc b/internal/ceres/partitioned_matrix_view_test.cc
index d01dd3b..6e9dd92 100644
--- a/internal/ceres/partitioned_matrix_view_test.cc
+++ b/internal/ceres/partitioned_matrix_view_test.cc
@@ -31,13 +31,13 @@
 #include "ceres/partitioned_matrix_view.h"
 
 #include <memory>
+#include <random>
 #include <vector>
 
 #include "ceres/block_structure.h"
 #include "ceres/casts.h"
 #include "ceres/internal/eigen.h"
 #include "ceres/linear_least_squares_problems.h"
-#include "ceres/random.h"
 #include "ceres/sparse_matrix.h"
 #include "glog/logging.h"
 #include "gtest/gtest.h"
@@ -50,7 +50,6 @@
 class PartitionedMatrixViewTest : public ::testing::Test {
  protected:
   void SetUp() final {
-    srand(5);
     std::unique_ptr<LinearLeastSquaresProblem> problem =
         CreateLinearLeastSquaresProblemFromId(2);
     CHECK(problem != nullptr);
@@ -65,11 +64,16 @@
         options, *down_cast<BlockSparseMatrix*>(A_.get()));
   }
 
+  double RandDouble() { return distribution_(prng_); }
+
   int num_rows_;
   int num_cols_;
   int num_eliminate_blocks_;
   std::unique_ptr<SparseMatrix> A_;
   std::unique_ptr<PartitionedMatrixViewBase> pmv_;
+  std::mt19937 prng_;
+  std::uniform_real_distribution<double> distribution_ =
+      std::uniform_real_distribution(0.0, 1.0);
 };
 
 TEST_F(PartitionedMatrixViewTest, DimensionsTest) {
diff --git a/internal/ceres/random.h b/internal/ceres/random.h
deleted file mode 100644
index 0495d67..0000000
--- a/internal/ceres/random.h
+++ /dev/null
@@ -1,73 +0,0 @@
-// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
-// http://ceres-solver.org/
-//
-// Redistribution and use in source and binary forms, with or without
-// modification, are permitted provided that the following conditions are met:
-//
-// * Redistributions of source code must retain the above copyright notice,
-//   this list of conditions and the following disclaimer.
-// * Redistributions in binary form must reproduce the above copyright notice,
-//   this list of conditions and the following disclaimer in the documentation
-//   and/or other materials provided with the distribution.
-// * Neither the name of Google Inc. nor the names of its contributors may be
-//   used to endorse or promote products derived from this software without
-//   specific prior written permission.
-//
-// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
-// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-// POSSIBILITY OF SUCH DAMAGE.
-//
-// Author: keir@google.com (Keir Mierle)
-//         sameeragarwal@google.com (Sameer Agarwal)
-
-#ifndef CERES_INTERNAL_RANDOM_H_
-#define CERES_INTERNAL_RANDOM_H_
-
-#include <cmath>
-#include <cstdlib>
-
-#include "ceres/internal/export.h"
-
-namespace ceres {
-
-inline void SetRandomState(int state) { srand(state); }
-
-inline int Uniform(int n) {
-  if (n) {
-    return rand() % n;
-  } else {
-    return 0;
-  }
-}
-
-inline double RandDouble() {
-  auto r = static_cast<double>(rand());
-  return r / RAND_MAX;
-}
-
-// Box-Muller algorithm for normal random number generation.
-// http://en.wikipedia.org/wiki/Box-Muller_transform
-inline double RandNormal() {
-  double x1, x2, w;
-  do {
-    x1 = 2.0 * RandDouble() - 1.0;
-    x2 = 2.0 * RandDouble() - 1.0;
-    w = x1 * x1 + x2 * x2;
-  } while (w >= 1.0 || w == 0.0);
-
-  w = sqrt((-2.0 * log(w)) / w);
-  return x1 * w;
-}
-
-}  // namespace ceres
-
-#endif  // CERES_INTERNAL_RANDOM_H_
diff --git a/internal/ceres/schur_eliminator_benchmark.cc b/internal/ceres/schur_eliminator_benchmark.cc
index 2bb18ca..0973e73 100644
--- a/internal/ceres/schur_eliminator_benchmark.cc
+++ b/internal/ceres/schur_eliminator_benchmark.cc
@@ -29,13 +29,13 @@
 // Authors: sameeragarwal@google.com (Sameer Agarwal)
 
 #include <memory>
+#include <random>
 
 #include "Eigen/Dense"
 #include "benchmark/benchmark.h"
 #include "ceres/block_random_access_dense_matrix.h"
 #include "ceres/block_sparse_matrix.h"
 #include "ceres/block_structure.h"
-#include "ceres/random.h"
 #include "ceres/schur_eliminator.h"
 
 namespace ceres::internal {
@@ -91,8 +91,10 @@
 
     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] = RandNormal();
+      values[i] = standard_normal(prng);
     }
 
     b_.resize(matrix_->num_rows());
diff --git a/internal/ceres/schur_eliminator_test.cc b/internal/ceres/schur_eliminator_test.cc
index a8c5f71..da27019 100644
--- a/internal/ceres/schur_eliminator_test.cc
+++ b/internal/ceres/schur_eliminator_test.cc
@@ -31,6 +31,7 @@
 #include "ceres/schur_eliminator.h"
 
 #include <memory>
+#include <random>
 
 #include "Eigen/Dense"
 #include "ceres/block_random_access_dense_matrix.h"
@@ -41,7 +42,6 @@
 #include "ceres/detect_structure.h"
 #include "ceres/internal/eigen.h"
 #include "ceres/linear_least_squares_problems.h"
-#include "ceres/random.h"
 #include "ceres/test_util.h"
 #include "ceres/triplet_sparse_matrix.h"
 #include "ceres/types.h"
@@ -283,8 +283,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] = RandNormal();
+    values[i] = standard_normal(prng);
   }
 
   Vector b(matrix.num_rows());
diff --git a/internal/ceres/sparse_cholesky_test.cc b/internal/ceres/sparse_cholesky_test.cc
index 0522232..e99d59b 100644
--- a/internal/ceres/sparse_cholesky_test.cc
+++ b/internal/ceres/sparse_cholesky_test.cc
@@ -32,6 +32,7 @@
 
 #include <memory>
 #include <numeric>
+#include <random>
 #include <vector>
 
 #include "Eigen/Dense"
@@ -42,7 +43,6 @@
 #include "ceres/internal/config.h"
 #include "ceres/internal/eigen.h"
 #include "ceres/iterative_refiner.h"
-#include "ceres/random.h"
 #include "glog/logging.h"
 #include "gmock/gmock.h"
 #include "gtest/gtest.h"
@@ -165,17 +165,18 @@
 class SparseCholeskyTest : public ::testing::TestWithParam<Param> {};
 
 TEST_P(SparseCholeskyTest, FactorAndSolve) {
-  SetRandomState(2982);
   const int kMinNumBlocks = 1;
   const int kMaxNumBlocks = 10;
   const int kNumTrials = 10;
   const int kMinBlockSize = 1;
   const int kMaxBlockSize = 5;
+  std::mt19937 prng;
+  std::uniform_real_distribution distribution(0.1, 1.0);
 
   for (int num_blocks = kMinNumBlocks; num_blocks < kMaxNumBlocks;
        ++num_blocks) {
     for (int trial = 0; trial < kNumTrials; ++trial) {
-      const double block_density = std::max(0.1, RandDouble());
+      const double block_density = distribution(prng);
       Param param = GetParam();
       SparseCholeskySolverUnitTest(::testing::get<0>(param),
                                    ::testing::get<1>(param),
diff --git a/internal/ceres/triplet_sparse_matrix.cc b/internal/ceres/triplet_sparse_matrix.cc
index 86a444e..c5ff2e5 100644
--- a/internal/ceres/triplet_sparse_matrix.cc
+++ b/internal/ceres/triplet_sparse_matrix.cc
@@ -32,12 +32,12 @@
 
 #include <algorithm>
 #include <memory>
+#include <random>
 
 #include "ceres/compressed_row_sparse_matrix.h"
 #include "ceres/crs_matrix.h"
 #include "ceres/internal/eigen.h"
 #include "ceres/internal/export.h"
-#include "ceres/random.h"
 #include "ceres/types.h"
 #include "glog/logging.h"
 
@@ -282,8 +282,8 @@
   }
 }
 
-std::unique_ptr<TripletSparseMatrix>
-TripletSparseMatrix::CreateFromTextFile(FILE* file) {
+std::unique_ptr<TripletSparseMatrix> TripletSparseMatrix::CreateFromTextFile(
+    FILE* file) {
   CHECK(file != nullptr);
   int num_rows = 0;
   int num_cols = 0;
@@ -317,16 +317,19 @@
   std::vector<int> rows;
   std::vector<int> cols;
   std::vector<double> values;
+  std::mt19937 prng;
+  std::uniform_real_distribution uniform01(0.0, 1.0);
+  std::uniform_real_distribution standard_normal;
   while (rows.empty()) {
     rows.clear();
     cols.clear();
     values.clear();
     for (int r = 0; r < options.num_rows; ++r) {
       for (int c = 0; c < options.num_cols; ++c) {
-        if (RandDouble() <= options.density) {
+        if (uniform01(prng) <= options.density) {
           rows.push_back(r);
           cols.push_back(c);
-          values.push_back(RandNormal());
+          values.push_back(standard_normal(prng));
         }
       }
     }