Monte Carlo simulation to approximate the value of PI

Question

This is a C++ implementation of a simple Monte Carlo simulation to approximate the value of pi. The program uses the standard library Mersenne twister engine to generate two random numbers between -1.0 to 1.0, one for x and one for y, for each point generated.

Every time the Point class constructor is called it checks to see if the values passed in satisfy the equation:

$$x^2 + y^2 = 1$$

If it does, both the number of points and the number of points within the circle increase by one. If not, only the total number of points is increased.

The program operates on the property that a circle inscribed within a unit square has an area of pi/4 so I divide the number of points in the circle by the total number of points and then multiply by four to get the approximate value of pi.

/** Monte Carlo Simulation to estimate the value of PI
 *
 */

unsigned long long totalPoints = 0;
unsigned long long pointsInCircle = 0;

unsigned long long overallTotalPoints = 0;
unsigned long long overallPointsInCircle = 0;

#include "includes\Main.hpp"

const long double actualPI = 3.14159265359;

long double approximatePI(unsigned long long inCircle, unsigned long long total);
void printPercentError(const long double myPI, const long double actualPI);


int main()
{
    system("cls");

    int zeroCounter = 0;
    int currentIteration = 0;

    long double myPI = 0;

    // Test Parameters
    const int M = 50;                   // Simulations per program execution
    const int N = 20'000'000;           // Points per simulation

    for (int i = 0; i < M; ++i) {    

        totalPoints = 0;
        pointsInCircle = 0;

        srand(time(NULL));

        static ::std::random_device rd;
        static std::mt19937_64 rng(rand());

        static std::uniform_real_distribution<double> uid1(-1.0,1.0);   // x
        static std::uniform_real_distribution<double> uid2(-1.0,1.0);   // y


        for (int j = 0; j < N; ++j) {
            Point p([&]{
                return uid1(rng);
            }(), [&]{
                return uid2(rng);
            }());
        }

        std::cout << "Iteration: " << currentIteration + 1 << std::endl;
        std::cout << "Points in circle: " << pointsInCircle << std::endl;
        std::cout << "Total points: " << totalPoints << std::endl;


        printf("PI: %.10Le\n", approximatePI(pointsInCircle, totalPoints));
        myPI += approximatePI(overallPointsInCircle, overallTotalPoints);

        ++currentIteration;

        // Divide myPI every 5 iterations to put off an overflow for as long as possible
        if (zeroCounter % 5 == 0) {
            if (currentIteration == M) {
                // Do not zero out myPI
            } else {
                myPI /= 5;
                zeroCounter = 0;

                overallTotalPoints = 0;
                overallPointsInCircle = 0;
            }
        }

        printf("\n");
    }

    printf("\n\n");

    std::cout << "Points in circle overall: " << overallPointsInCircle << std::endl;
    std::cout << "Total points overall: " << overallTotalPoints << std::endl;

    myPI = approximatePI(overallPointsInCircle, overallTotalPoints);
    printf("PI: %.10Le\n", myPI);
    printPercentError(myPI, actualPI);

    return EXIT_SUCCESS;
}


long double approximatePI(const unsigned long long inCircle, const unsigned long long total)
{
    return (((long double) inCircle / (long double) total) * 4);
}


void printPercentError(const long double myPI, const long double actualPI)
{
    printf("\n\tPercent Error: %.10Le%%\n", 100 * ((myPI - actualPI) / actualPI));
}

Output:

Iteration: 1
Points in circle: 15708507
Total points: 20000000
PI: 3.1417014000e+000

Iteration: 2
Points in circle: 15708384
Total points: 20000000
PI: 3.1416768000e+000

Iteration: 3
Points in circle: 15705582
Total points: 20000000
PI: 3.1411164000e+000

Iteration: 4
Points in circle: 15710275
Total points: 20000000
PI: 3.1420550000e+000

Iteration: 5
Points in circle: 15706656
Total points: 20000000
PI: 3.1413312000e+000

Iteration: 6
Points in circle: 15705987
Total points: 20000000
PI: 3.1411974000e+000

Iteration: 7
Points in circle: 15709017
Total points: 20000000
PI: 3.1418034000e+000

Iteration: 8
Points in circle: 15706029
Total points: 20000000
PI: 3.1412058000e+000

Iteration: 9
Points in circle: 15707857
Total points: 20000000
PI: 3.1415714000e+000

Iteration: 10
Points in circle: 15709423
Total points: 20000000
PI: 3.1418846000e+000

Iteration: 11
Points in circle: 15708768
Total points: 20000000
PI: 3.1417536000e+000

Iteration: 12
Points in circle: 15708054
Total points: 20000000
PI: 3.1416108000e+000

Iteration: 13
Points in circle: 15706587
Total points: 20000000
PI: 3.1413174000e+000

Iteration: 14
Points in circle: 15706242
Total points: 20000000
PI: 3.1412484000e+000

Iteration: 15
Points in circle: 15705368
Total points: 20000000
PI: 3.1410736000e+000

Iteration: 16
Points in circle: 15706917
Total points: 20000000
PI: 3.1413834000e+000

Iteration: 17
Points in circle: 15708899
Total points: 20000000
PI: 3.1417798000e+000

Iteration: 18
Points in circle: 15706711
Total points: 20000000
PI: 3.1413422000e+000

Iteration: 19
Points in circle: 15709541
Total points: 20000000
PI: 3.1419082000e+000

Iteration: 20
Points in circle: 15710444
Total points: 20000000
PI: 3.1420888000e+000

Iteration: 21
Points in circle: 15706853
Total points: 20000000
PI: 3.1413706000e+000

Iteration: 22
Points in circle: 15712085
Total points: 20000000
PI: 3.1424170000e+000

Iteration: 23
Points in circle: 15710104
Total points: 20000000
PI: 3.1420208000e+000

Iteration: 24
Points in circle: 15709046
Total points: 20000000
PI: 3.1418092000e+000

Iteration: 25
Points in circle: 15707144
Total points: 20000000
PI: 3.1414288000e+000

Iteration: 26
Points in circle: 15707643
Total points: 20000000
PI: 3.1415286000e+000

Iteration: 27
Points in circle: 15707230
Total points: 20000000
PI: 3.1414460000e+000

Iteration: 28
Points in circle: 15708895
Total points: 20000000
PI: 3.1417790000e+000

Iteration: 29
Points in circle: 15708037
Total points: 20000000
PI: 3.1416074000e+000

Iteration: 30
Points in circle: 15706278
Total points: 20000000
PI: 3.1412556000e+000

Iteration: 31
Points in circle: 15707066
Total points: 20000000
PI: 3.1414132000e+000

Iteration: 32
Points in circle: 15705930
Total points: 20000000
PI: 3.1411860000e+000

Iteration: 33
Points in circle: 15708086
Total points: 20000000
PI: 3.1416172000e+000

Iteration: 34
Points in circle: 15710387
Total points: 20000000
PI: 3.1420774000e+000

Iteration: 35
Points in circle: 15706093
Total points: 20000000
PI: 3.1412186000e+000

Iteration: 36
Points in circle: 15707002
Total points: 20000000
PI: 3.1414004000e+000

Iteration: 37
Points in circle: 15709261
Total points: 20000000
PI: 3.1418522000e+000

Iteration: 38
Points in circle: 15709620
Total points: 20000000
PI: 3.1419240000e+000

Iteration: 39
Points in circle: 15705817
Total points: 20000000
PI: 3.1411634000e+000

Iteration: 40
Points in circle: 15708118
Total points: 20000000
PI: 3.1416236000e+000

Iteration: 41
Points in circle: 15706553
Total points: 20000000
PI: 3.1413106000e+000

Iteration: 42
Points in circle: 15710051
Total points: 20000000
PI: 3.1420102000e+000

Iteration: 43
Points in circle: 15706290
Total points: 20000000
PI: 3.1412580000e+000

Iteration: 44
Points in circle: 15706712
Total points: 20000000
PI: 3.1413424000e+000

Iteration: 45
Points in circle: 15707375
Total points: 20000000
PI: 3.1414750000e+000

Iteration: 46
Points in circle: 15707948
Total points: 20000000
PI: 3.1415896000e+000

Iteration: 47
Points in circle: 15709443
Total points: 20000000
PI: 3.1418886000e+000

Iteration: 48
Points in circle: 15707079
Total points: 20000000
PI: 3.1414158000e+000

Iteration: 49
Points in circle: 15707642
Total points: 20000000
PI: 3.1415284000e+000

Iteration: 50
Points in circle: 15706250
Total points: 20000000
PI: 3.1412500000e+000



Points in circle overall: 15706250
Total points overall: 20000000
PI: 3.1412500000e+000

    Percent Error: -1.0907002523e-002%

Answer 1

Here are a few possibilities to consider:

Avoid globals

Global variables have been known as a source of problems and errors for a long time. They're particularly pernicious when (as here) something makes relatively hidden changes to them. In this case, you haven't shown us your Point class, but it's apparent that creating a Point actually modifies your totalPoints, pointsInCircle, overallTotalPoints, and overallPointsInCircle global variables.

In most real code, I'd consider fixing this an emergency. Nearly the only thing that's more serious than this is a report from a customer that our code had lost or destroyed his data.

There are (I guess) times that globals become hard to avoid. In a few cases (especially things like configuration data that's read in while initializing, and never modified again, they aren't necessarily even a terrible thing--but that's not the kind of situation you have here.

Seeding the random number generator

It turns out to be somewhat difficult to seed a random number generator correctly. A previous question started with the right general idea, and some of the answers included minor improvements to the code in the question. Personally, I wasn't entirely excited about some the choices it embodied (though, in fairness, they weren't terrible decisions by any means).

Anyway, after some thinking on it, I'd written up a small class for generating a reasonably random seed for a C++11 random number generator:

template <class Rand> 
class Seed {
    class seeder {
        std::array < std::random_device::result_type, Rand::state_size > rand_data;
    public:
        seeder() {
            std::random_device rd;
            std::generate(rand_data.begin(), rand_data.end(), std::ref(rd));
        }

        typename std::array < std::random_device::result_type, Rand::state_size >::iterator begin() { return rand_data.begin(); }
        typename std::array < std::random_device::result_type, Rand::state_size >::iterator end() { return rand_data.end(); }
    } seed;

    std::seed_seq s;

public:
    Seed() : s(seed.begin(), seed.end()) { }

    template <class I>
    auto generate(I a, I b) { return s.generate(std::forward<I>(a), std::forward<I>(b)); }
};

As you can see, there's a bit more here than we'd really like to deal with directly most of the time. Technically speaking, this isn't absolutely required to work correctly--in particular, the seed sequence model includes a couple of things we don't provide (but virtually nobody ever uses either). We could provide them easily enough if we wanted, but they mostly don't matter so I'm omitting them at least for now.

Using that, we can create and seed our random number generator something like this:

    using Rand = std::mt19937_64;
    Seed<Rand> seed;
    Rand rng(seed);

Strictly speaking, the Rand isn't necessary, but I find it kind of handy (and it makes it easy to switch between generators, if we feel like trying that).

That brings up the next major point though:

Don't repeat yourself

You have a number things repeated throughout your code. For example, you have similar code for reporting the number of points in the circle, number of total points, computed value of Pi, etc., both inside the loop for each iteration, and after the loop to give the final value.

In such a case, you almost always want to move that code into a function that you call from two locations, rather than having separate code for each.

Avoid unnecessary code

You have a number of cases where you seem to have thrown in code that accomplishes little or nothing. For example, you have two random number distribution objects, but they're stateless, and by nature of the fact that you're dealing with a circle, they're guaranteed to be identical. You might as well just use one.

Likewise, when you create a Point object, you pass lambdas to compute random values instead of just passing the random values themselves. This doesn't appear to accomplish anything.

You also have return EXIT_SUCCESS; at the end of main. This accomplishes nothing--this will happen automatically if you don't return some other value from main.

Avoid hacks

Comments like this:

    // Divide myPI every 5 iterations to put off an overflow for as long as possible

...pretty much send chills down my spine. I'm not sure exactly what overflow we're avoiding. Given what we're doing, I hard a hard time imagining how anything should reasonably overflow. Likewise this:

    myPI += approximatePI(overallPointsInCircle, overallTotalPoints);

...looks rather problematic as well. A number of approximations of Pi, added together, shouldn't give Pi. It should give approximately Pi times the number of approximations. We need to divide by the number of approximations to get that back to an approximation. This isn't emergency-level, but it's still pretty worrisome.

Generally avoid mixing C and C++ I/O

Right now you switch off between using printf and std::cout to produce your output. You seem to gain little from doing so. I'd pick one and stick with it.

Avoid back-slashes in header paths

This code:

#include "includes\Main.hpp"

...isn't portable. It's generally considered better to use a forward slash instead:

#include "includes/Main.hpp"

This way, the code works on most current OSes (e.g., Windows, Linux, Mac OS) without modification.

Consider making generic code more explicitly generic

For example, your printPercentError uses names related to Pi for its parameters--but it applies equally well to many other quantities, not just Pi. It can also apply to various different types--might as well make it a template.

Separation of Concerns

I dislike functions like printPercentError that combine two fundamentally different types of operations--in this case, computing an error, and writing that value to standard output. I'd much rather have one piece of code that just computed an error, and separate code to print that out.

Conclusion

Putting all those together, we might end up with code something like this (note: I've reduced the number of iterations, for the sake of quicker testing):

/** Monte Carlo Simulation to estimate the value of PI
*
*/
#include <random>
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <iterator>
#include <array>

template <class T>
class Point {
    T x;
    T y;
public:
    Point(T x, T y) :x(x), y(y) {}

    T mag() { return std::hypot(x, y); }
};

template <class T>
T approximatePI(T inCircle, T total)
{
    return ((inCircle / total) * T(4));
}

template <class T>
T percentError(T computed, T actual)
{
     return T(100) * ((computed - actual) / actual);
}

template <class Rand> 
class Seed {
    class seeder {
        std::array < std::random_device::result_type, Rand::state_size > rand_data;
    public:
        seeder() {
            std::random_device rd;
            std::generate(rand_data.begin(), rand_data.end(), std::ref(rd));
        }

        typename std::array < std::random_device::result_type, Rand::state_size >::iterator begin() { return rand_data.begin(); }
        typename std::array < std::random_device::result_type, Rand::state_size >::iterator end() { return rand_data.end(); }
    } seed;

    std::seed_seq s;

public:
    Seed() : s(seed.begin(), seed.end()) { }

    template <class I>
    auto generate(I a, I b) { return s.generate(std::forward<I>(a), std::forward<I>(b)); }
};

template <class Real>
void report(unsigned long long inCircle, unsigned long long total) {
    std::cout << "Points in circle: " << inCircle << '\n';
    std::cout << "Total points: " << total << '\n';
    std::cout << "Estimated Pi: " << std::setprecision(10) << approximatePI<Real>(inCircle, total) << "\n";
}

int main()
{
    // Test Parameters
    const int M = 10;                   // Simulations per program execution
    const int N = 2'000'000;           // Points per simulation

    unsigned long long overallPointsInCircle = 0;
    unsigned long long overallPoints = 0;

    using Real = double;
    Real actualPI = 4 * std::atan(Real(1));

    for (int i = 0; i < M; ++i) {
        int pointsInCircle = 0;
        int totalPoints = 0;

        using Rand = std::mt19937_64;

        Seed<Rand> seed;

        Rand rng(seed);

        std::uniform_real_distribution<Real> uid(-1.0, 1.0);

        for (int j = 0; j < N; ++j) {
            Point<Real> p{ uid(rng), uid(rng) };

            ++totalPoints;
            ++overallPoints;

            if (p.mag() <= 1.0) {
                ++pointsInCircle;
                ++overallPointsInCircle;
            }
        }

        std::cout << "\nIteration: " << i + 1 << std::endl;
        report<Real>(pointsInCircle, totalPoints);
    }

    std::cout << "\n\nFinal values:\n";
    report<Real>(overallPointsInCircle, overallPoints);

    auto computedPi = approximatePI<Real>(overallPointsInCircle, overallPoints);
    std::cout << "Percent Error: " << std::setprecision(10) << percentError(computedPi, actualPI);
}