Monte Carlo Method of finding pi using C

Question

I have written a function that takes in a long long value n and uses that as the number of iterations to go through. The function should give a good estimate of pi, however, all the values for large n tends towards 3.000, and not 3.1415,so I am not sure what is going on?

Is there anything that I did wrong?

This is my code:

double estimate_pi(long long n){
    double randomx, randomy, equation, pi;
    long long i, incircle = 0;

    for(i = 0; i < n; i++){
        randomx = (double)(rand() % (1+1-0) + 0);
        randomy = (double)(rand() % (1+1-0) + 0);

        equation = randomx * randomx + randomy * randomy;

        if(equation <= 1){
            incircle++;
        }
    }

    pi = (long double)4 * (long double)incircle / (long double)n;

    return pi;
}

in the main function, to print 10 values of pi:

int main(void){

    long long N;
    double pi_approx;
    int i;

    printf("Input a value of N: ");
    if(scanf("%ld", &N) != 1){
        printf("Error, input must be an integer!\n");
        exit(EXIT_SUCCESS);
    } 
    if(N < 1){
        printf("Error, the integer must be positive!\n");
        exit(EXIT_SUCCESS); 
    }

    srand(time(NULL));
    for(i = 0; i < 10; i++){
        pi_approx = estimate_pi(N);
        printf("%.10f\n", pi_approx);
    }
    return 0;
}

score 6 · Answer 1 · answered Mar 13 '19 at 18:29

It works as it should. The problem is the implementation.

The C rand() function returns an integer in the range 0 to RAND_MAX. The keyword there is integer.

You then calculate the result of that integer modulo 2, which can be 0 or 1. That leaves you with 4 possible points: (0,0), (0,1), (1,0), (1,1).

Of those 4 points, only 1 lies outside of the circle of radius 1: (1,1). That is, of out of 4 possible points, 3 lie in the circle.

You should replace that code to use floating point values, not integers, so you calculate the proportion of points inside and outside of the circle.

Thank you so much for your explanation, this has helped me understand this a bit more! ^^ — Pengibaby, Mar 13 '19 at 19:17

alinsoar · Accepted Answer · 2019-03-13T19:02:08.187

5

You need to use floating-point randomization, or otherwise to use a circle with a very big radius.

So instead of

    randomx = (double)(rand() % (1+1-0) + 0);
    randomy = (double)(rand() % (1+1-0) + 0);

you use

    randomx = rand();
    randomy = rand();

and you consider if it falls inside the circle of radius RAND_MAX

   #define RMAX ((double)RAND_MAX*(double)RAND_MAX)
   equation <= RMAX;

You do the details. Read man 3 rand to see that rand() returns integer.

edited Mar 13 '19 at 19:02

answered Mar 13 '19 at 18:29

alinsoar

15,386
4
57
74

Thank you so much! This is my first time writing C program, so still unfamiliar with it. But this really helped! ^^ – Pengibaby Mar 13 '19 at 19:16

score 3 · Answer 3 · answered Mar 13 '19 at 18:29

Your randomx and randomy variables are constrained to an integer value, because the rand() functions returns an integer.

See it live here.

As a consequence, each your two variables will be either 1 or 0, so your point will be randomly one of (0,0), (1,0), (0,1), (1,1), which has a 3:4 chance of being in the circle. Hence your result of 3.

You can look up How to generate random float number in C if you want a random number between 0 and 1.

Thank you for the link, and your explanation! ^^ – Pengibaby Mar 13 '19 at 19:18 — Pengibaby, Mar 13 '19 at 19:18

score 2 · Answer 4 · answered Mar 13 '19 at 18:31

2

For your approach to work, you need to generate double values drawn from a uniform distribution on the interval [0,1] (or approximately so). You are instead generating random integers drawn from the two-element set {0, 1}, and converting them to type double. This does not yield anything remotely like the distribution you need for your purposes.

answered Mar 13 '19 at 18:31

John Bollinger

160,171
8
81
157

Thank you for your explanation! ^^ – Pengibaby Mar 13 '19 at 19:20

score 0 · Answer 5 · answered Jan 27 '21 at 20:05

The others already pointed out your mistake. I am providing a optimized version based on another SO Thead

   int main(void)
    {
        long points = 1000000000; // Some input
        long m = 0;
        unsigned long HAUSNUMERO = 1;
        double DIV1byMAXbyMAX = 1. / RAND_MAX / RAND_MAX;
    
        unsigned int aThreadSpecificSEED_x = HAUSNUMERO + 1
        unsigned int aThreadSpecificSEED_y = HAUSNUMERO - 1
        for(long i = 0; i < points; i++)
        {
            double x = rand_r( &aThreadSpecificSEED_x );
            double y = rand_r( &aThreadSpecificSEED_y );
            m += (1  >= ( x * x + y * y ) * DIV1byMAXbyMAX);
        }
        printf("Pi is roughly %lf\n", (double) 4*m / (double) points);
    }

Monte Carlo Method of finding pi using C

5 Answers5