Estimate pi using monte-carlo method varying radius and number of samples.

pi.cpp

// Compile with:
//    clang++ --std=c++14 -O3 -o pi pi.cpp

#include <algorithm>
#include <cmath>
#include <experimental/iterator>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <type_traits>
#include <random>
#include <vector>


// Reference value for pi
const auto pi{3.14159265359};
// Numeric type derived from pi
using numeric = std::remove_const_t<decltype(pi)>;


// A 2d point within the square
using point_type = std::pair<numeric, numeric>;
// A vector of samples used to calculate pi
using points_type = std::vector<point_type>;


// Create a vector of samples
auto random_points(std::size_t samples, numeric radius) {
    // We need to be able to generate some random points
    std::random_device random_device;
    std::mt19937 engine{random_device()};
    std::uniform_real_distribution<numeric> location{-radius, radius};
    auto make_point = [&]() {
        return point_type{location(engine), location(engine)};
    };
    // Create our samples
    points_type points;
    points.reserve(samples); // ~30% slower without this line
    std::generate_n(std::back_inserter(points), samples, make_point);
    return points;
}


// Calculate the error from pi
auto calculate_error(std::size_t samples, numeric radius) {
    // Get our points
    auto points = random_points(samples, radius);
    // Number of points inside the circle gives our estimate
    auto inside = std::count_if(points.begin(), points.end(), [&](auto p) {
        return p.first * p.first + p.second * p.second < radius * radius;
    });
    auto estimate = inside * numeric{4} / samples;
    // Return error from pi
    return std::abs(pi - estimate);
}


int main() {
    std::cout << std::setprecision(6) << std::fixed;
    // Calculate the errors and print as we go
    std::size_t n{};
    auto newline = std::ostream_iterator<const char *>(std::cout, "\n");
    std::generate_n(newline, 10, [&]() {
        std::size_t m{};
        auto line = std::experimental::make_ostream_joiner(std::cout, "  ");
        std::generate_n(line, 8, [&]() {
            return calculate_error(std::pow(10, m++), std::pow(10, n++));
        });
        return "";
    });

    return 0;
}