diff --git a/src/distribution.cpp b/src/distribution.cpp
index 7f5b498add..cbdd131273 100644
--- a/src/distribution.cpp
+++ b/src/distribution.cpp
@@ -552,14 +552,9 @@ double Tabular::sample_unbiased(uint64_t* seed) const
   double c = prn(seed);
 
   // Find first CDF bin which is above the sampled value
-  double c_i = c_[0];
-  int i;
-  std::size_t n = c_.size();
-  for (i = 0; i < n - 1; ++i) {
-    if (c <= c_[i + 1])
-      break;
-    c_i = c_[i + 1];
-  }
+  auto c_iter = std::lower_bound(c_.begin() + 1, c_.end(), c);
+  int i = std::distance(c_.begin(), c_iter) - 1;
+  double c_i = c_[i];
 
   // Determine bounding PDF values
   double x_i = x_[i];
diff --git a/tests/cpp_unit_tests/test_distribution.cpp b/tests/cpp_unit_tests/test_distribution.cpp
index e46088a630..479e7c8a7a 100644
--- a/tests/cpp_unit_tests/test_distribution.cpp
+++ b/tests/cpp_unit_tests/test_distribution.cpp
@@ -82,6 +82,31 @@ TEST_CASE("Test alias sampling method for pugixml constructor")
   }
 }
 
+TEST_CASE("Test sampling a large linear-linear tabular distribution")
+{
+  constexpr int n_points = 10001;
+  constexpr int n_samples = 200000;
+  openmc::vector<double> x(n_points);
+  openmc::vector<double> p(n_points);
+  for (int i = 0; i < n_points; ++i) {
+    x[i] = static_cast<double>(i) / (n_points - 1);
+    p[i] = 2.0 * x[i];
+  }
+
+  openmc::Tabular dist(
+    x.data(), p.data(), n_points, openmc::Interpolation::lin_lin);
+  uint64_t seed = openmc::init_seed(0, 0);
+
+  double mean = 0.0;
+  for (int i = 0; i < n_samples; ++i) {
+    mean += dist.sample(&seed).first;
+  }
+  mean /= n_samples;
+
+  // The normalized PDF is 2x on [0, 1], which has a mean of 2/3.
+  REQUIRE_THAT(mean, Catch::Matchers::WithinAbs(2.0 / 3.0, 0.003));
+}
+
 TEST_CASE("Test construction of SpatialBox with parameters")
 {
   openmc::Position ll {-1, -2, -3};