diff --git a/examples/bal_problem.cc b/examples/bal_problem.cc
index ceac89a..ccf7449 100644
--- a/examples/bal_problem.cc
+++ b/examples/bal_problem.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -30,22 +30,22 @@
 
 #include "bal_problem.h"
 
+#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 {
-namespace examples {
+namespace ceres::examples {
 namespace {
-typedef Eigen::Map<Eigen::VectorXd> VectorRef;
-typedef Eigen::Map<const Eigen::VectorXd> ConstVectorRef;
+using VectorRef = Eigen::Map<Eigen::VectorXd>;
+using ConstVectorRef = Eigen::Map<const Eigen::VectorXd>;
 
 template <typename T>
 void FscanfOrDie(FILE* fptr, const char* format, T* value) {
@@ -55,15 +55,14 @@
   }
 }
 
-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();
   }
 }
 
 double Median(std::vector<double>* data) {
-  int n = data->size();
-  std::vector<double>::iterator mid_point = data->begin() + n / 2;
+  auto mid_point = data->begin() + data->size() / 2;
   std::nth_element(data->begin(), mid_point, data->end());
   return *mid_point;
 }
@@ -73,12 +72,12 @@
 BALProblem::BALProblem(const std::string& filename, bool use_quaternions) {
   FILE* fptr = fopen(filename.c_str(), "r");
 
-  if (fptr == NULL) {
+  if (fptr == nullptr) {
     LOG(FATAL) << "Error: unable to open file " << filename;
     return;
   };
 
-  // This wil die horribly on invalid files. Them's the breaks.
+  // This will die horribly on invalid files. Them's the breaks.
   FscanfOrDie(fptr, "%d", &num_cameras_);
   FscanfOrDie(fptr, "%d", &num_points_);
   FscanfOrDie(fptr, "%d", &num_observations_);
@@ -111,7 +110,7 @@
   if (use_quaternions) {
     // Switch the angle-axis rotations to quaternions.
     num_parameters_ = 10 * num_cameras_ + 3 * num_points_;
-    double* quaternion_parameters = new double[num_parameters_];
+    auto* quaternion_parameters = new double[num_parameters_];
     double* original_cursor = parameters_;
     double* quaternion_cursor = quaternion_parameters;
     for (int i = 0; i < num_cameras_; ++i) {
@@ -137,7 +136,7 @@
 void BALProblem::WriteToFile(const std::string& filename) const {
   FILE* fptr = fopen(filename.c_str(), "w");
 
-  if (fptr == NULL) {
+  if (fptr == nullptr) {
     LOG(FATAL) << "Error: unable to open file " << filename;
     return;
   };
@@ -161,8 +160,8 @@
     } else {
       memcpy(angleaxis, parameters_ + 9 * i, 9 * sizeof(double));
     }
-    for (int j = 0; j < 9; ++j) {
-      fprintf(fptr, "%.16g\n", angleaxis[j]);
+    for (double coeff : angleaxis) {
+      fprintf(fptr, "%.16g\n", coeff);
     }
   }
 
@@ -297,14 +296,20 @@
   CHECK_GE(point_sigma, 0.0);
   CHECK_GE(rotation_sigma, 0.0);
   CHECK_GE(translation_sigma, 0.0);
-
+  std::mt19937 prng;
+  std::normal_distribution<double> 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<double> rotation_noise_distribution(0.0,
+                                                               point_sigma);
+  std::normal_distribution<double> 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);
     }
   }
 }
@@ -331,5 +338,4 @@
   delete[] parameters_;
 }
 
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples