Calculate Pi using Monte Carlo

Question

#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

Answer 1

One could consider at least the following points:

• Instead of including <stdlib.h>, I'd include <cstdlib>.

• In getRandDart(), it might in this case be more readable to do static_cast<double>(rand()) / RAND_MAX; instead of multiplying by 1.0.

• In the for loop, all of x, y and d can be const, so I'd make them const. This has the potential to protect the programmer from unintended mistakes, and can sometimes allow the compiler to optimize better.

• When you increment by one (in countInSquare += 1;), it makes more sense to use the ++ operator, i.e., to just write ++countInSquare. This is more idiomatic and protects us from unintended mistakes: ++ conveys the meaning of increment (by one), whereas with += we might accidentally write += 2; and that would be perfectly valid (but not what we wanted).

• Regardless of the above point, notice that during the for-loop, it holds that iteration == countInSquare. So strictly speaking, the variable countInSquare is unnecessary and could be replaced by just iteration when needed.

• You could consider making the number of iterations and the second operand of the % operand constants to allow for easier modification and perhaps to slightly improve readability.

• Instead of typing (4.0 * countInCircle / countInSquare) twice, we could make a function that takes the two variables as parameters. This would allow us to save some typing, and again to protect us from unintended mistakes.

Answer 2

#include <random>
std::default_random_engine              generator;

If you require portable and reproducible results from your random engine, consider specifying the exact random engine you want instead of relying on an implementation-defined alias.

// GCC & Clang
using minstd_rand0 = linear_congruential_engine<uint_fast32_t, 16807, 0, 2147483647>;
using default_random_engine = minstd_rand0;

// MSVC
using default_random_engine = mt19937;

---

std::uniform_real_distribution<double>  distribution(0,1);

For distributions, the standard doesn't require the result to be reproducible and identical across all implementations. If you care for portability, you'll need to write/copy the specific distribution behavior you want.

---

inline double getRandDart() {return distribution(generator);}

// 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)

This is an example of where the comment doesn't match the function. The function getRandDart() doesn't return a dart. It returns a single magnitude instead of a Dart/Vec2d.

---

    for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {

Did you intend for ten billion and one iterations (<=)?

---

    long     countInSquare   = 0;
    for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {
        countInSquare   += 1;
        if (iteration % 10'000'000 == 0) {
             std::cout << iteration << " " << (4.0 * countInCircle / countInSquare) << "\n";

Do you need countInSquare? Both iteration and countInSquare are maintaining the same count.

When you don't care for the size in which an integer can represent because the value is countable, simply use int. In this case, you went with long as you needed a larger integral type to hold a value that couldn't be represented by an int. Compiling in an environment where long is 32-bit (windows) would obviously be bad. In these cases, use a specific fixed width integer type from <cstdint>. auto also works in determining the correct type (64-bit long on gcc/clang, 64-bit long long on msvc).

---

        countInCircle   += (d >= 1.0) ? 0 : 1;

Compilers are nice nowadays and optimized your ternary add.

        countInCircle   += (d < 1.0);

---

You can tighten the inner loop by reorganizing the operations. Rather than checking every iteration to report, tile the iterations into groups that calculate between reports.

    constexpr auto totalIterations = 10'000'000'001;
    constexpr auto reportInterval = 10'000'000;
    constexpr auto reports = std::div(totalIterations, reportInterval);

    while (reports.quot--) {
        for (auto iteration = 0; iteration < reportInterval; ++iteration) {
            const auto x = getRandDart();
            const auto y = getRandDart();
            const auto d = x * x + y * y;
            countInCircle += (d < 1.0);
        }
        countInSquare += reportInterval;
        std::cout << countInSquare << " " << (4.0 * countInCircle / countInSquare) << "\n";
    }

    while (reports.rem--) {
        const auto x = getRandDart();
        const auto y = getRandDart();
        const auto d = x * x + y * y;
        countInCircle += (d < 1.0);
    }

    std::cout << "\n\n" << std::setprecision(9) << (4.0 * countInCircle / totalIterations) << "\n";