C  Calculating pi by Monte Carlo integration using working
C  arrays of 8K entries.
C
C  File:  pi.f

C  Tally n cycles of random samples:

      integer function sample(n)
      integer n

      integer, parameter :: V = 8192
      integer k, IR(V), INC(V)
      real, array(V) :: X, Y

      INCLUDE 'mpfrand.h'

      INC = 0      ! Counts number of points inside circle
      do, k = 1, n
        call mpf_random(IR)
        X = IR/Two31
        call mpf_random(IR)
        Y = IR/Two31
        where (X*X + Y*Y <= 1.0) INC = INC + 1
      end do
      sample = sum(INC)
      end function sample

C  Run samples and tally results

      program montepi

      INTERFACE
        integer function sample(n)
        integer n
      END INTERFACE

      INCLUDE 'mpfrand.h'
      integer, parameter :: V = 8192

      integer ng, gs, Vgs, inc, tinc, SR(V), seed, k, time, ttime
      real pi, sigma, tpi, tsigma, sqrgs

      print *, "Specify <= 0 groups to quit."
      do while (.TRUE.)
         print '(" Number of groups, group size, seed > "$)'
         read *, ng, gs
         if ( gs <= 0 ) go to 99    ! Exit
         tinc = 0     ! In circle all groups so far
         ttime = 0    ! Time for all groups so far
         Vgs = V*gs
         sqrgs = sqrt(real(Vgs))
         SR = seed + [1:V]      ! Seed vector
         print *, "  Trials    pi      sigma      tpi      tsigma    time"
         call mpf_srandom(SR)
         do, k = 1, ng          ! ng groups
            call mpTimerStart()
            inc = sample(gs)    ! Number in circle
            time = mpTimerElapsed()
            ttime = ttime + time
            pi =  real(inc)/real(Vgs)  ! Estimate for this group
            sigma = sqrt(pi - pi*pi)/sqrgs
            tinc = tinc + inc   ! Total so far in circle
            tpi = real(tinc)/real(k*Vgs)
            tsigma = sqrt((tpi - tpi*tpi)/real(k*Vgs))
            pi = 4.0*pi
            sigma = 4.0*sigma
            tpi = 4.0*tpi
            tsigma = 4.0*tsigma
            print 21, k*gs, pi, sigma, tpi, tsigma, ttime/1000.0
 21         FORMAT (1X, I7, 2(1X, F7.5, 1X, E11.5), 1X, F7.5)
        end do
      end do
 99   continue   ! For exit
      end program montepi
