#include <iostream>
#include <iomanip>
#ifdef USE_OLD_RAND
#include <stdlib.h>
inline double getRandDart() {return rand() * 1.0 / RAND_MAX;}
#else
#include <random>
std::default_random_engine generator;
std::uniform_real_distribution<double> distribution(0,1);
inline double getRandDart() {return distribution(generator);}
#endif
// Monte Carlo Simulator to estimate the value of PI.
//
// If we have a circle with a radius of 1.
// Then the smallest square that encloses the circle as sides of length 2.
//
// Area of circle pi r^2 = pi
// Area of square 2.r.2.r = 4
//
// Ratio of overlapping area: pi/4
//
// If we throw darts randomly at a dart board (with an even distribution) and always hit the square.
// Then the ratio of darts falling into the circle should be pi/4 of the total number of darts thrown.
//
// pi/4 * countInSquare = countInCircle
//
// pi = 4 . countInCircle / countInSquare
//
// To simplify the maths.
// We will set the center point as 0,0 and only use the top right quadrant of the circle.
// We have 1/4 the size of the square and circle but the same maths still apply.
//
// A dart thrown has a random x/y value in the range 0->1 (top right quadrant).
// A dart is outside the circle if x^2 + y^2 > 1 (note 1^2 is 1)
//
int main()
{
long countInSquare = 0;
long countInCircle = 0;
for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {
double x = getRandDart();
double y = getRandDart();
double d = (x * x) + (y * y);
countInSquare += 1;
countInCircle += (d >= 1.0) ? 0 : 1;
if (iteration % 10'000'000 == 0) {
std::cout << iteration << " " << (4.0 * countInCircle / countInSquare) << "\n";
}
}
std::cout << "\n\n" << std::setprecision(9) << (4.0 * countInCircle / countInSquare) << "\n";
}
Output:
> ./a.out
9990000000 3.14158
10000000000 3.14158
3.14158355