SUBROUTINE random(x)
  IMPLICIT NONE
  DOUBLE PRECISION, INTENT (out) :: x(:)

  CALL RANDOM_NUMBER(x)

  x = 2.0d0 * x - 1.0d0
END SUBROUTINE random


FUNCTION pi4(N)
  IMPLICIT NONE

  INTEGER, INTENT (in) :: N
  DOUBLE PRECISION :: pi4

  INTERFACE 
     SUBROUTINE random(x)
       IMPLICIT NONE
       DOUBLE PRECISION, INTENT (out) :: x(:)
     END SUBROUTINE random
  END INTERFACE

  DOUBLE PRECISION :: x(N), y(N), R(N)
  INTEGER :: M, i

  CALL random(x)
  CALL random(y)

  R = x*x +y*y

  M = COUNT(R < 1.0d0) 

  pi4 = DBLE(M) / DBLE(N)
END FUNCTION pi4


PROGRAM buffon
  IMPLICIT NONE

  INTERFACE
     FUNCTION pi4(N)
       IMPLICIT NONE
       INTEGER, INTENT (in) :: N
       DOUBLE PRECISION :: pi4
     END FUNCTION pi4
  END INTERFACE

  INTEGER :: N, i
  DOUBLE PRECISION, PARAMETER :: pi = 3.14159265358979323846d0

  DO i = 2, 8
     N = 10**i
     PRINT *, N, pi4(N), pi/4d0
  ENDDO

END PROGRAM buffon
