2

I was asked to write a program that gets N points from a user and using a continuous distribution finds an approximation of Pi using Monte Carlo technique. This is what I wrote:

        unsigned seed = chrono::steady_clock::now().time_since_epoch().count();
        default_random_engine e (seed);
        uniform_real_distribution<> dist(0,2);
        int N = atoi(argv[1]);
        int inside = 0;
        long double appPi = 0;
        for (int i = 0; i<N; i++){
            double x = dist(e);
            double y = dist(e);
            double distance = sqrt(x*x+y*y);
            if (distance <= 1){ inside++;}
        }
        appPi = (inside/N)*4;

However after printing appPi all I get is 0. I think that algorithm by itself is ok? as it prints plausible values of x and y, but it doesn't really work for me.

Jakub Sapko
  • 304
  • 2
  • 15
  • 6
    `int` divided by `int` is an `int` – Xatyrian Mar 23 '20 at 10:06
  • Oh, that's right. Is there a way to change it to double? – Jakub Sapko Mar 23 '20 at 10:07
  • 1
    `appPi = (static_cast(inside)/N)*4;` – Yksisarvinen Mar 23 '20 at 10:09
  • Thank you, this seems to have solved the problem of division, however I'm getting kind of strange results, for example 0.7872 which is nowhere near the real PI. Is there something I did wrong? – Jakub Sapko Mar 23 '20 at 10:11
  • Why are you testing values of x and y distributed over the range (0,2)? Shouldn't that be (0,1)? Also, you can make your code run faster if you recognise that if sqrt(x\*x + y\*y) <= 1, then x\*x + y\*y <= 1. There is no need to calculate the square root. – r3mainer Mar 23 '20 at 12:13

1 Answers1

2

In addition to the integer division pointed out by Xatyrian, you lack a multiplicative factor. You are extracting random points in a square of size l = 2, and then you count how many points lie in a quarter of a circle of radius R = 1. If we define the fraction of such points f we can connect this value to the areas of the square and the quarter of a circle: pi R^2 / 4 = f l^2.

If we plug in this relation the values defined above we find pi = 16 f and not 4 f as your code seems to imply. Indeed, 0.7872 * 4 = 3.1488.

A quick and more reasonable fix than using 16 instead of 4 is to extract points in a square of size l by making the following change:

uniform_real_distribution<> dist(0,1);
lr1985
  • 1,085
  • 1
  • 9
  • 21