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.
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).