Monte Carlo integration to find pi with a certain precision in FORTRAN

Question

I'm taking a course in Numerical Methods and I've been requested to implement the famous Monte Carlo algorithm to find pi that you can find here.

I had no difficulties in writing the code with an arbitrary number of trials:

REAL(8) FUNCTION distance(xvalue, yvalue) RESULT(dist)
IMPLICIT NONE
REAL(8), INTENT(in) :: xvalue, yvalue
dist = SQRT(xvalue**2 + yvalue**2)
END FUNCTION distance 

PROGRAM ass2
  IMPLICIT NONE

  INTEGER, DIMENSION(1) :: SEED
  REAL(8) :: p, x, y

  REAL(8), EXTERNAL :: distance

  REAL(8) :: pi_last, pi
  INTEGER :: npc, npt, i

  npc = 0
  npt = 0
  pi = 1.0

  SEED(1) = 12345
  CALL RANDOM_SEED

  DO i=1, 1000000000

     CALL RANDOM_NUMBER(p)
     x = p
     CALL RANDOM_NUMBER(p)
     y = p

     npt = npt + 1

     IF (distance(x, y) < 1.0)  THEN
          npc = npc + 1
     END IF

     pi_last = pi
     pi = 4.0*(npc*1.0)/(npt*1.0)

  END DO

  PRINT*, 'Pi:', pi

END PROGRAM ass2

I noticed that it converges approximately as sqrt(N of steps). Now I have to stop the algorithm at a certain precision, so I created an endless DO loop with an EXIT inside an IF statement:

    REAL(8) FUNCTION distance(xvalue, yvalue) RESULT(dist)
IMPLICIT NONE
REAL(8), INTENT(in) :: xvalue, yvalue
dist = SQRT(xvalue**2 + yvalue**2)
END FUNCTION distance 

PROGRAM ass2
  IMPLICIT NONE

  INTEGER, DIMENSION(1) :: SEED
  REAL(8) :: p, x, y

  REAL(8), EXTERNAL :: distance

  REAL(8) :: pi_last, pi
  INTEGER :: npc, npt, i

  npc = 0
  npt = 0
  pi = 1.0

  SEED(1) = 12345
  CALL RANDOM_SEED

  DO

     CALL RANDOM_NUMBER(p)
     x = p
     CALL RANDOM_NUMBER(p)
     y = p

     npt = npt + 1

     IF (distance(x, y) < 1.0)  THEN
          npc = npc + 1
     END IF

     pi_last = pi
     pi = 4.0*(npc*1.0)/(npt*1.0)

  IF ( ABS(pi - pi_last) < 0.000001 .AND. pi - pi_last /= 0)    THEN
    EXIT
  END IF

  END DO

  PRINT*, 'Pi:', pi

END PROGRAM ass2

The problem is that this returns a value of pi that doesn't have the precision I asked for. I get the logic behind it: if I get two consecutive values far from pi but close to each other, the condition will be satisfied and the program will exit the DO statement. The problem is that I don't know how to modify it in order to get a precision decided by me. So the question is:

How do I implement this algorithm in a manner such that I can decide the precision of pi in output?

EDIT: Ok, I implemented both of your solutions and they work, but only for 10^(-1), 10^(-3) and 10^(-5). I think it's a problem of the pseudorandom sequence if for 10^(-2) and 10^(-4) it returns an incorrect value of pi.

Answer 1

It's not enough to specify a desired precision -- you also need to allow some chance that the precision target is not met. Then you can solve for the number of trials in (e.g.) Hoeffding's inequality to meet your desired precision with the desired probability (as you've observed, n needs to be about the square root of one over the precision to succeed with constant probability).

Answer 2

In an ideal setting (perfect random number generator generating real numbers in the mathematical sense) your variable npc is a random variable with binomial distribution B(n,π/4) where n is npt from your program. Its expected value is n*π/4, so you correctly compute the approximation of π as pi=4*npc/npt. Now this approximation can take all values from 0 to 4 no matter how many loop iterations you calculate, because npc can take all values from 0 to npt. For a range around π you can only give a probability (using c as shorthand for npc; P denotes the probability of an event):

P(|pi - π| < d) = P(-d < pi - π < d) = P(-d < 4*c/n - π < d) = P(n*(π-d)/4 < c < n*(π+d)/4) =

= P(c < n*(π+d)/4) - P(c < n*(π-d)/4) ~=

~= F_N(n*(π+d)/4) - F_N(n*(π-d)/4) = 2F(d*√(n/(π(4-π))))-1

where F_N is the probability function of the normal distribution N(nπ/4;nπ/4(1-π/4)) which approximates the above binomial distribution and F is the probability function of the standard normal distribution. Now given a deviation d and a probability p you can compute n s.t. the last term does not go below p:

n = ceil(π(4-π)(F^-1((p+1)/2)/d)^2))

Then with n loop iterations you can compute the approximation pi of π to a desired precision with given probability. If we want to achieve a probability of p=99%, then the above formula simplifies to

n ~= 17.89/d² ,

so for precision d=0.0001 roughly n=1.789E9 iterations are necessary!

Note: since a computer cannot meet the above ideal setting there is also a (theoretical) limit on the precision you can reach with this algorithm. There are only finitely many floating point numbers representable in a computer, so your points (x,y) lie on a kind of grid. The best approximation of π that can be computed with this algorithm boils down to performing your loop over all grid points in [0,1]x[0,1]. The good old C-function rand() has a resolution of 31 bits (at least in the VS stdlib). So it does not make sense to compute more than n=31² points, which gives a maximal precision of √(17.89/n) = 1.97E-9 when asking for 99% correctness.