Pi day challange - http://www.fluentcpp.com/2017/03/02/the-pi-day-challenge-for-expressive-code-in-c/

estimate_pi.cpp

#include <cmath>
#include <functional>
#include <iomanip>
#include <iostream>
#include <random>
using namespace std;
int main() {
    double pi_ref=3.14159265359, r_max=1000000000.0, n_max=10000000.0;
    auto axis = bind(uniform_real_distribution<double>(0.0, 2*r_max), mt19937(random_device()()));
    auto estimate_pi = [&](double r, double n) -> double {
        double quantity_in_circle=0.0;
        for(double i=1.0;i<=n;++i) quantity_in_circle+=(hypot(fmod(axis(),2*r)-r,fmod(axis(),2*r)-r)<r);
        return 4.0*quantity_in_circle/n;
    };
    for(double r=1.0;r<=r_max;r*=10.0) for(double n=1;n<=n_max;n*=10.0) {
        cout<<setw(5)<<fabs(pi_ref-estimate_pi(r,n))<<setw(5)<<'\t'<<(n==n_max?"\n":"");
    }
    return EXIT_SUCCESS;
}


/*
axis() - randomly generate [0,2*r_max)
fmod(axis(),2*r) - randomly generate [0,2*r)
hypot(x-r, y-r) - distance of point (x,y) from centre of square (and circle).
hypot(x-r, y-r)<r - resolves to 1.0 if inside circle and 0.0 if not.


A recursive alternative to estimate_pi that would have worked except that
n_max = 10^7 is far too large to avoid stack overflow without manually
adjusting the stack space in the OS:

function<double(double,int)> estimate_pi = [&](double r, double n) -> double { return fabs(n)<0.5 ? 0 : ((n-1)/n)*estimate_pi(r,n-1)+(4/n)*(hypot(fmod(axis(),2*r)-r,fmod(axis(),2*r)-r)<r); };

*/