Benchmark two different versions of MC-based π approximation. A result looks like this: Loop: (7.431e+07 ± 140019) ns Vector: (8.50605e+07 ± 159916) ns.

pi_benchmark.cc

#include <algorithm>
#include <chrono>
#include <cmath>
#include <functional>
#include <iostream>
#include <random>
#include <utility>
#include <vector>


typedef double Number;
typedef std::pair<Number, Number> Point;
typedef std::mt19937 RandomGenerator;


Number ApproximatePiVector(const size_t N, RandomGenerator &rng) noexcept {
  static std::uniform_real_distribution<Number> dist(0.0, 1.0);
  std::vector<Point> v(N);
  std::generate(std::begin(v), std::end(v), [&]() {
    return Point{dist(rng), dist(rng)};
  });
  const size_t n = std::count_if(std::begin(v), std::end(v),
    [](const Point &p) { return p.first * p.first + p.second * p.second <= 1; }
  );
  return 4.0 * n / N;
}


Number ApproximatePiLoop(const size_t N, std::mt19937 &rng) noexcept {
  static std::uniform_real_distribution<Number> dist(0.0, 1.0);
  size_t i,
         n = 0;
  Number x,
         y;
  for (i = 0; i < N; i++) {
    x = dist(rng);
    y = dist(rng);
    if (x * x + y * y <= 1) {
      ++n;
    }
  }
  return 4.0 * n / N;
}


std::pair<double, double> Benchmark(
    const size_t num_runs,
    std::function<Number(const size_t, std::mt19937&)> f,
    const size_t num_iterations,
    std::mt19937 &rng) {
  size_t i;
  std::vector<double> timings(num_runs);
  for (i = 0; i < num_runs; i++) {
    std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();
    f(num_iterations, rng);
    std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
    auto time_span = std::chrono::duration_cast<std::chrono::nanoseconds>(end - start);
    timings[i] = time_span.count();
  }
  // Calc mean and sem
  double mean = std::accumulate(std::begin(timings), std::end(timings), 0.0)
              / num_runs;
  double sem = sqrt(std::accumulate(std::begin(timings), std::end(timings),
      0.0, [&](const double acc, const double x) {
        return acc + (x - mean) * (x - mean);
      }
  ) / (num_runs * (num_runs - 1)));
  return std::make_pair(mean, sem);
}


int main() {
  std::random_device rd;
  std::mt19937 rng(rd());
  auto bench_loop = Benchmark(100, ApproximatePiLoop, 1000000, rng);
  auto bench_vector = Benchmark(100, ApproximatePiVector, 1000000, rng);

  std::cout << "Loop: (" << bench_loop.first << " ± " << bench_loop.second
            << ") ns" << std::endl;
  std::cout << "Vector: (" << bench_vector.first << " ± "
            << bench_vector.second << ") ns" << std::endl;
}