CODEWITHSUNDEEP
c++mathrandomprobabilitymontecarlo
Zestyorangez
Zestyorangez

Reputation: 17

why does my program not approximate pi?

I'm trying to write a program that approximates pi by using the probability that two random numbers are coprime which is 6/pi^2 so i should be able to aproximate pi as sqrt(6/probability) but for some reason it doesn't tend towards pi rather it seems to tend towards 3.912. i have tried the rand() function and i have tried the mersenne twister random number generator but they both give me similar results. what's going on? here is my code in c++

#include <iostream>
#include <ctime>
#include <cstdlib>
#include <cmath>
#include <algorithm>

using namespace std;

int main()
{
    srand(time(0));
    int times;

    cout << "enter number of times: ";
    cin >> times;

    for(int i = 0; i <= times; i ++)
    {
        double pi_approx, probability;
        int num1, num2, number_of_coprime;
        num1 = rand();
        num2 = rand();
        if(__gcd(num1, num2)>1)
        {
           number_of_coprime++;
        }

        probability = double(number_of_coprime)/double(i);
        pi_approx = sqrt(6/probability);

        cout <<100*i/times<< "% number of coprimes: "<< number_of_coprime << " pi approximation: " << pi_approx<<endl;

    }

    return 0;
}

Upvotes: 0

Views: 112

Answers (2)

Arndt Jonasson
Arndt Jonasson

Reputation: 854

You are counting the wrong thing: you should increment the counter when they are coprime (gcd == 1), but you increment when the gcd is > 1.

Upvotes: 1

Yola
Yola

Reputation: 19063

You need to declare variables in a proper scope, and the right formula is that probability goes to 6/pi^2 as number of times goes to infinity. Try the following code:

int main()
{
    srand(time(0));

    //for (int times = 100; times < 1000; ++times) 
    int times = 1000000;
    {
        double pi_approx, probability;
        int number_of_coprime = 0;
        for (int i = 0; i < times; i++)
        {
            int num1, num2;
            num1 = rand() % times;
            num2 = rand() % times;
            if (__gcd(num1, num2) == 1) // increment if coprime !!
            {
                number_of_coprime++;
            }
        }
        probability = double(number_of_coprime) / double(times);
        pi_approx = sqrt(6 / probability);
        cout << pi_approx << endl;
    }

    return 0;
}

Output is:

3.14179

Upvotes: 2

Related Questions