C  Calculating volume of the unit ball in d dimensions by Monte
C  Carlo methods, using  arrays of 8K entries.  Integration is
C  carried out over the unit cube [0..1]^d .
C
C  File:  pid.f

C  Tally n cycles of random samples for dimension d:

      integer function sample(n, d)
      integer n, d

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

      INCLUDE 'mpfrand.h'

      INC = 0      ! Counts number of points inside circle
      do, k = 1, n
        R2 = 0.0      ! Accumulates SUM X^2
        do, r = 1, d
          call mpf_random(IR)
          X  = IR/Two31
          R2 = R2 + X*X
        end do ! r
        where (R2 <= 1.0) INC = INC + 1
      end do
      sample = sum(INC)
      end function sample

C  Run samples and tally results

      program montepid

      INTERFACE
        integer function sample(n, d)
        integer n, d
      END INTERFACE

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

      integer d, ng, gs, Vgs, inc, tinc, SR(V), seed, k, time, ttime
      real vd, sigma, tvd, tsigma, sqrgs, twod

      print *, "Specify dimensions <= 0 to quit."
      do while (.TRUE.)
         print '(" Dimensions > "$)'
         read *, d
         if ( d <= 0 ) go to 99
         print '(" Number of groups, group size, seed > "$)'
         read *, ng, gs, seed
         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
         twod = 2.0**d          ! 2^d for scaling results
         sqrgs = sqrt(real(Vgs))
         SR = seed + [1:V]      ! Seed vector
         print *, "  Trials   v(d)     sigma     tv(d)     tsigma     time"
         call mpf_srandom(SR)
         do, k = 1, ng          ! ng groups
            call mpTimerStart()
            inc = sample(gs, d)    ! Number in circle
            time = mpTimerElapsed()
            ttime = ttime + time
            vd =  real(inc)/real(Vgs)  ! Estimate for this group
            sigma = sqrt(vd - vd*vd)/sqrgs
            tinc = tinc + inc   ! Total so far in circle
            tvd = real(tinc)/real(k*Vgs)  ! Estimate so far
            tsigma = sqrt((tvd - tvd*tvd)/real(k*Vgs))
            vd = twod*vd
            sigma = twod*sigma
            tvd = twod*tvd
            tsigma = twod*tsigma
            print 21, k*gs, vd, sigma, tvd, tsigma, ttime/1000.0
 21         FORMAT (1X, I7, 2(1X, F7.5, 1X, E11.5), 1X, F9.5)
        end do
      end do
 99   continue   ! For exit
      end program montepid

