<<Back

Lab 1: Estimation of Constant PI (== 3.14...) using Monte Carlo Simulation Techniques

Requirements

1. Understand how to generate pseudorandom numbers under Linux using C/C++.

2. Implement two algorithms to estimate the constant PI (== 3.14 ...) using Monte Carlo simulation techniques.

3. Compare the two algorithms and tell which one is better.

Random Number Generation under Linux 

See "man random" for information about how to use the standard random number generator library function under Linux. We will discuss the algorithms used by such functions in detail, later in the quarter, but for now we need to understand the followings:

1. Repeated calls to the function return "pseudorandom" (which means algorithmically generated, but random "appearing") integer values. For example, the following program:

    #include <iomanip.h>
    #include<stdlib.h>
    int main(){
        cout << random() << endl << random() << endl;
    }

produces this output:

    846930886
    1804289383

2. Because these values are algorithmically generated, you get the same sequence of values every time you run the program. If you want to change this, you must call the function srandom() to change the "seed" of the random number generator.

3. These integer values are uniformly distributed over the range from 0 to RAND_MAX.  Thus the expression random()/RAND_MAX will generate pseudorandom numbers uniformly distributed between 0 and 1.

Algorithm 1

Repeat the following calculation N times:

a) generate x and y as two uniform random real values between 0 and 1

b) use these two values as the coordinates of a point randomly distributed over the unit square, which goes from (0,0) on the bottom left corner to (1,1) on the top right corner.

c) test whether this random point is also inside a unit circle (i.e., radius 1, center at 0,0) by finding its distance from the origin:

   distance = x^2 + y^2

d) if distance <= 1 then increment a success counter, S, otherwise do nothing.

After completing the N iterations, form the ratio (S / N).  Since the area of the square is exactly 1, while the area of the quarter of the unit circle that is contained within that square is PI/4, the value of ratio should be approximately (PI/4)/(1), or equivalently

    PI~= 4 * S / N

After implementing this algorithm, they should run it several times for n=10, 100, 1000 etc up to 1 million to see how well it works.

Algorithm 2

Instead of looking at the overlap between a 1/4 circle and the unit square, let's split the square diagonally from upper left [i.e., point (0,1)] to lower right [i.e., point (1,0)]. Notice that the boundary of the 1/4 circle is entirely contained in the upper triangle [i.e., away from the origin point (0,0)] so we can restrict our experiment to generating points in that triangle.

Notice that the equation of the diagonal straight line forming the lower boundary of the triangle is

     x+y=1

and whenever x+y<1 then we have that:

    (1-x) + (1-y)
= 2 - (x+y)
> 2 - 1
= 1

Furthermore, the area of the upper triangle is clearly 1/2, whereas the area of the "lens" shaped intersection of the 1/4 circle and the upper triangle is

    PI/4 - 1/2  == (PI - 2) / 4

Thus, if we modify step (b) of the algorithm 1 like this:

b) use these two values as the coordinates of a point randomly distributed over the unit square, which goes from (0,0) on the bottom left corner to (1,1) on the top right corner.  However, if x+y < 1, then modify them as follows:
    x = 1-x
    y = 1-y

Then the ratio (S/N) becomes an estimate of [(PI - 2) / 4] / [(1 / 2)], or equivalently:

    PI ~= 2 * (S/N) + 2

Notice that with this second version of the algorithm, the estimate of PI must be between 2 and 4 (rather than between 0 and 4).

<<Back

Last Update: 04/22/2004