Skip to content

Instantly share code, notes, and snippets.

@KayEss
Last active March 5, 2017 06:12
Show Gist options
  • Save KayEss/e250cec9d5d4bc423087bcc6c11cb7bf to your computer and use it in GitHub Desktop.
Save KayEss/e250cec9d5d4bc423087bcc6c11cb7bf to your computer and use it in GitHub Desktop.
Estimate pi using monte-carlo method varying radius and number of samples.
// 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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment