i have to modify this program to calculate the PI. I understand basicaly that the monte carlo methods are used to calculate numbers in a statistical simulation, but i don't know how to apply the methods to the problem. here is the program:
// MONCARLO.CPP
//
// This program calculates the perimeter of an ellipse.
// One quarter of the perimeter of an ellipse with the
// larger radius 1 is equal to the so-called
// "elliptic integral" E -- the area under the curve
// y = sqrt(1 - (k*sin(x))^2) on the interval 0 <= x <= pi/2.
// The parameter k describes the shape of the ellipse and is
// usually expressed in terms of the so-called "modular angle":
// k = sin(alpha), where alpha is the angle between the minor
// axis and the line that connects its end to a focus of
// the ellipse. For a circle, alpha = 0 and k = 0.
// The integral is calculated by using the Monte Carlo method.
//
// Author: B. Random
// Date: 2/29/1999
//
#include <iostream.h>
#include <iomanip.h>
#include <stdlib.h> // Declares rand() and RAND_MAX
#include <math.h> // Declares sin(x)
const double PI = 3.14159265358979323846;
// Modular angle (in degrees and radians):
const double alpha = 45;
const double alpha_radians = alpha * PI / 180.;
// k defines elongation of the ellipse:
const double k = sin(alpha_radians);
//************************************************** **************
double MonteCarlo (long nPoints)
// This function calculates the elliptic integral E.
// It generates "nPoints" (pseudo)random points in the
// rectangle { 0 <= x <= pi/2; 0 <= y <= 1 }, finds
// the fraction of all points that fall under the curve,
// and calculates the corresponding fraction of the area of
// the rectangle.
{
double x, y, f;
long n, count = 0;
for (n = nPoints; n > 0; n--) {
// Generate a random point:
// 0 <= x <= pi/2
// 0 <= y <= 1
x = double(rand()) * (.5 * PI / RAND_MAX);
y = double(rand()) / RAND_MAX;
// Calculate 1 - (k*sin(x))^2
f = k * sin(x);
f *= f; // Square f
f = 1 - f;
// Increment the counter if the point (x,y) is on or under
// the curve:
if (y*y <= f) count++;
// (more efficient than: if (y <= sqrt(f)) ...)
}
// count E
// ------- = -----------------------------------------
// nPoints area of rect {0 <= x <= pi/2; 0 <= y <= 1}
return (double(count) / double(nPoints)) * (.5 * PI);
}
//************************************************** **************
int main()
{
long n;
double ellipsePerim;
cout << "Enter the number of random points in Monte Carlo ==> ";
cin >> n;
ellipsePerim = 4. * MonteCarlo(n);
cout << "The perimeter of the ellipse with a modular angle of "
<< alpha << " degrees\n"
<< " and a major axis of 2 units is approximately = "
<< setprecision(3) << ellipsePerim << endl;
return 0;
}