CODEWITHSUNDEEP
cmontecarlopi
Pengibaby
Pengibaby

Reputation: 373

Monte Carlo Method of finding pi using C

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;
}

Upvotes: 7

Views: 3635

Answers (5)

dreamcrash
dreamcrash

Reputation: 51583

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);
    }

Upvotes: 0

alinsoar
alinsoar

Reputation: 15813

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.

Upvotes: 5

John Bollinger
John Bollinger

Reputation: 181419

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.

Upvotes: 2

Guille
Guille

Reputation: 343

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.

Upvotes: 6

Louen
Louen

Reputation: 3677

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.

Upvotes: 3

Related Questions