Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Monte Carlo Method of finding pi using C

Tags:

c

pi

montecarlo

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;
}
like image 651
Pengibaby Avatar asked Mar 13 '19 18:03

Pengibaby


4 Answers

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.

like image 56
Guille Avatar answered Nov 10 '22 05:11

Guille


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.

like image 23
alinsoar Avatar answered Nov 10 '22 07:11

alinsoar


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.

like image 44
Louen Avatar answered Nov 10 '22 07:11

Louen


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.

like image 2
John Bollinger Avatar answered Nov 10 '22 06:11

John Bollinger