Monte Carlo Pi

Hi all, I am very new to this so any help would be greatly appreciated. I am currently trying to write a program that determines Pi using the Monte Carlo method and compares this to the actual value of Pi. The user inputted values a and b represent the limits of the total number of points with a step size of 100. My code works for the first value inputted but not for the rest, I feel that this is because the value of r++ does not reset after the first loop. A push in the right direction and any general comments about my style of coding would be massively helpful, thanks.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
  
#include <iostream>
#include <iomanip>
#include <cmath>
#include <ctime>
#include <cstdlib>
#define Pi 2.0*acos(0.0)

using namespace std; //avoids having to uses std:: with cout and cin

//calculates the distance of the point from the origin
double dist(double x, double y)
{
    double DistFromOrigin;
    DistFromOrigin = sqrt((x*x) + (y*y));
    //Euclidean Distance Formula from point (0,0)
    return DistFromOrigin;
}
//calculates pi using number of points inside circle versus total points
double piMC(int d, double r)
{
    double MonteCarloPi = (r / d ) * 4;
    return MonteCarloPi;
}

int main (int argc, char* argv[])
{
	system("COLOR F0"); //inverts colours
	
	
	//declare variables
	int a, b;
    double coord_x, coord_y, test;
   
    srand(time(NULL));      //creates the seed time
    test = rand() / double(RAND_MAX);
   
    double DistFromOrigin;       //distance from the origin
    double r = 0;           //number of points inside circle
	
	
	
	
	
	//Tell the user what is happening 
	cout <<"This program computes a value of Pi using the Monte Carlo mehtod";
	cout <<" and compares this value to the actual value of Pi.";
	cout <<"\nIt does this for a range of user inputted values from 80 to 850";
	cout <<" where these values represent the number of points where";
	cout << "\n0 < (x,y) < 1";
	cout << endl;
	
	//Ask user for input
	
	cout << "\n\nPlease enter a value for a, where (a >= 80): ";
	cin >> a; // enter value for a
		while (a<80 || a> 850)
	{
		cout << 
		"Value for a is not valid, please enter a value between 80 and 850: ";
		cin >> a; //enter value a
	}

	cout << "Please enter a value for b, where (a <= b <= 850): ";
	cin >> b; //enter value b
	
	
	while (b < a || b > 850)
	{
		cout << "Value for b is not valid, please enter a value higher than a";
		cout << " and less than 850: ";
		cin >> b;
	}
	
	//Print out headers 
	cout << "\n\n" << left << setfill ('d');
	cout << setw(20) << left << "Number of points ";
	cout << setw(20) << left << " $ Value of Pi ";
	cout << setw(20) << left << "    $ Difference versus Actual Pi ";
	cout << endl;
	
	for (int t = a; t <= b; t+=100)
	{
		
		
		for (int j = 1; j <= t; j++)
		{
			 coord_x = (double)rand() / (RAND_MAX);
       		 coord_y = (double)rand() / (RAND_MAX);
	  		 DistFromOrigin = dist(coord_x, coord_y);
       		 //counts how often the point is inside the circle
        	 
			 if(DistFromOrigin <= 1)
			 {
			 	r++;
			 }
		
	 	}
	 	double MonteCarloPi = piMC(t, r);
		double DifferencePi = Pi- MonteCarloPi;
		
		cout << scientific << left << setprecision(4) << setfill ('d')
		<< setw(20) << t << " $ "
		<< setw(20) << left << MonteCarloPi << " $ "
    	<< setw(20) << left << DifferencePi << " $ " 
		<< endl;	
	}
	
	//Let user see results before ending programme 
	cout << "\n\nPress enter to end the program";
	cin.get();
	cin.get();
	
}
Last edited on
any general comments about my style of coding would be massively helpful

Avoid non-descriptive variable names except when used as simple counter variable (e.g. i, j).
In the definition of your piMC function, your description seems fine, but I have no idea what t and d are supposed to mean. Give them meaningful names. Also, is r supposed to be "radius"? If so, why are you increasing the radius every time? I see from your comment that r is "number of points inside circle", but it would be better for this to be expressed in the variable name as well, ex: num_points.

cout << scientific << left << setprecision(4) << setfill ('d')
Why are you filling with the letter 'd'? It gives you output like this:

Number of points ddd $ Value of Pi ddddd    $ Difference versus Actual Pi 
90dddddddddddddddddd $ 2.8000e+00dddddddddd $ 3.4159e-01dddddddddd $ 


Also, to me at least, I don't understand your use of the 'a' and 'b' variables. What is the point of them? What should the user expect if they enter 90 and 100? They have no idea what this "step size" internal functionality is. Why not just give a single number to mean the # of iterations you want the program to run for?

Here's roughly how I would do it (clearly this isn't doing exactly what your program is doing, but does the monte carlo calcluations nonetheless). Feel free to disagree or get ideas from this:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
  
#include <iostream>
#include <iomanip>
#include <cmath>
#include <ctime>
#include <cstdlib>

using namespace std; //avoids having to uses std:: with cout and cin

const double Pi =  2.0 * acos(0.0);

//calculates the Euclidean distance of the point from the origin
double dist(double x, double y)
{
	return sqrt(x*x + y*y);
}

int main (int argc, char* argv[])
{
	srand(time(NULL));

	cout << "Enter number of points to generate: ";
	int num_iterations;
	cin >> num_iterations;
	
	int num_points_in_circle = 0;
	
	for (int i = 0; i < num_iterations; i++)
	{
	    double x_coord = static_cast<double>(rand()) / RAND_MAX;
	    double y_coord = static_cast<double>(rand()) / RAND_MAX;
	    
	    if (dist(x_coord, y_coord) < 1)
	    {
	        num_points_in_circle++;  	        
	    }
	    
	    // The ratio of the # of points inside the 1st quadrant of the circle to the # of iterations.
	    // converges to the ratio of the area of the 1st quadrant of the circle to the tatal square area.
	    // (or, Pi / 4)
	    double ratio = static_cast<double>(num_points_in_circle) / i;
	    
	    cout << "New estimate for Pi: " << 4 * ratio << "\n";
	}
	
	//Let user see results before ending programme 
	cout << "\n\nPress enter to end the program";
	cin.get();
	
}


Edit, note that the Monte Carlo method seems to converge extremely slowly... there are faster algorithms for calculating Pi that aren't probabilistic

New estimate for Pi: 3.14998
New estimate for Pi: 3.15003
New estimate for Pi: 3.15008
New estimate for Pi: 3.15013
New estimate for Pi: 3.15017
New estimate for Pi: 3.15022
New estimate for Pi: 3.15004
New estimate for Pi: 3.14986
New estimate for Pi: 3.14968
New estimate for Pi: 3.14973
New estimate for Pi: 3.14978
New estimate for Pi: 3.14983
New estimate for Pi: 3.14988
New estimate for Pi: 3.14993
New estimate for Pi: 3.14997
New estimate for Pi: 3.14979
New estimate for Pi: 3.14962
New estimate for Pi: 3.14966
New estimate for Pi: 3.14971
New estimate for Pi: 3.14976
New estimate for Pi: 3.14981
New estimate for Pi: 3.14986
New estimate for Pi: 3.14991
New estimate for Pi: 3.14995
New estimate for Pi: 3.15
New estimate for Pi: 3.15005
New estimate for Pi: 3.1501
New estimate for Pi: 3.14992
New estimate for Pi: 3.14997
New estimate for Pi: 3.15002
New estimate for Pi: 3.14984
New estimate for Pi: 3.14989


https://en.wikipedia.org/wiki/Category:Pi_algorithms
Last edited on
I'm somewhat confused about what you've done. To begin with, what is the purpose of deltaC?

Maybe reading this might help.
http://www.cplusplus.com/forum/unices/75447/#msg405443
Last edited on
Hi guys and thank you for the replies. I'll try to explain in more detail about what the program is supposed to do. The program requires two user inputted integers a and b. These integers represent the lower and upper limit for the total number of points used, increasing in steps of 100 each time. ie if I pick a = 100 and b = 500, the first loop estimates Pi for 100 points, then the second loop for 200 and so on until it has reached 500. I have not decided to do this, it is simply a requirement of the program I am asked to make. Again the fill character is simply a requirement of the program, likewise with the character $. In my head this is what I'm trying to get the program to do:

1) User inputs upper and lower limit (a and b)
2) A number of random points (a) are generated between 0 and 1 for x and y
3) The program then determines which points are inside the circle of radius 1
4) Calculates Monte Carlo Pi by dividing points in circle by total points and multiplies this by 4
5) This process is repeated for a value a+100 until this number reaches the upper limit b, with a different value for Monte Carlo Pi obtained for each total points value

I understand that this may not be the most useful method but again it is simply what I have been asked to do

Also thank you kbw for pointing out my deltaC, I have changed it above
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#include <iostream>
#include <iomanip>
#include <ctime>
#include <cstdlib>
using namespace std;

int main()
{
   const double RAD = RAND_MAX;
   const double RADSQ = RAD * RAD;
   double x, y;

   srand( time( 0 ) );
   cout << fixed << setprecision( 6 );
   cout << setw( 12 ) << "N" << "  " << setw( 12 ) << "pi" << '\n';

   for ( int N = 1000; N <= 1000000000; N *= 10 )
   {
      int success = 0;
      for ( int i = 0; i < N; i++ )
      {
         x = rand();
         y = rand();
         if ( x * x + y * y < RADSQ ) success++;
      }
      double pi = 4.0 * success / N;
      cout << setw( 12 ) << N << "  " << setw( 12 ) << pi << '\n';
   }
}


           N            pi
        1000      3.148000
       10000      3.152800
      100000      3.140160
     1000000      3.143264
    10000000      3.142640
   100000000      3.141521
  1000000000      3.141539


Even if rand() were unbiased, the standard deviation would be about 1.6/sqrt(N): you need a lot of points to get a good estimate.
Thank you lastchance, this is much simpler than mine and makes a lot more sense
Topic archived. No new replies allowed.