C  Calculating volume of the unit ball in d dimensions by Monte
C  Carlo with importance sampling, using  arrays of 8K entries and
C  a Moebius transformation to modify distribution.  Integration is
C  carried out over the unit cube [0..1]^d .
C
C  File:  pimprt.f

C  Tally n cycles of random samples for dimension d, using importance
C  sampling according to
C
C                 x
C     q(x) = ----------   where 0 < a < 1
C            ax + 1 - a
C
C  Results are:  hits = no. samples in ball
C                sf   = SUM F
C                sf2  = SUM F^2


      subroutine sample(n, d, a, hits, sf, sf2)
      integer n, d, hits
      real a, sf, sf2

      integer, parameter :: V = 8192
      integer k, r, IR(V), INC(V)
      real, array(V) :: X, Y, R2, qp, dn, SUMF, SUMF2
      real a1

      INCLUDE 'mpfrand.h'

      a1 = 1.0 - a ! For calculating q^-1(y)
      INC = 0      ! Counts number of points inside circle
      SUMF = 0.0   ! Sum for <F>
      SUMF2 = 0.0  ! Sum for <F^2>
      do, k = 1, n
        R2 = 0.0      ! Accumulates SUM X^2
        dn = 1.0      ! Forms PROD q'(x)
        do, r = 1, d
          call mpf_random(IR)
          Y  = IR/Two31           ! Y uniform in [0..1]
          X  = (a1*Y)/(1.0 - a*Y) ! X in [0..1]
          qp = a*X + a1
          dn = dn*a1/(qp*qp)      ! Next q'(x)
          R2 = R2 + X*X
        end do ! r
        where (R2 <= 1.0)
          INC = INC + 1           ! Tally hits
          qp  = 1.0/dn            ! F(y) = f(x1,...,xd)/(q'(x1)...q'(xd))
          SUMF = SUMF + qp        ! For <F>
          SUMF2 = SUMF2 + qp*qp   ! For <F^2>
        end where
      end do  ! k
      hits = SUM(INC)
      sf   = SUM(SUMF)            
      sf2  = SUM(SUMF2)
      end subroutine sample

C  Run samples and tally results

      program montepimprt

      INTERFACE
        subroutine sample(n, d, a, hits, sf, sf2)
        integer n, d, hits
        real a, sf, sf2
      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, a, sf, sf2, tsf, tsf2
      real kVgs

      print *, " Parameter a should satisfy 0 <= a < 1.  a = 0 gives"
      print *, " uniform sampling, a --> 1 is strongly biased."
      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, a > "$)'
         read *, ng, gs, seed, a
         if ( gs <= 0 ) go to 99    ! Exit
         tinc  = 0     ! In circle all groups so far
         ttime = 0     ! Time for all groups so far
         tsf   = 0.0   ! Sum of F so far
         tsf2  = 0.0   ! Sum of F^2 so far
         Vgs = V*gs   ! Number of samples/group
         twod = 2.0**d          ! 2^d for scaling results
         sqrgs = sqrt(real(Vgs))
         SR = seed + [1:V]      ! Seed vector
         print *,
     *   "  Trials    hits    v(d)     sigma     tv(d)     tsigma     time"
C          iiiiiii iiiiiiii fffffff eeeeeeeeeee fffffff eeeeeeeeeee fffffffff
         call mpf_srandom(SR)
         do, k = 1, ng          ! ng groups
            call mpTimerStart()
            call sample(gs, d, a, inc, sf, sf2)
            time = mpTimerElapsed()
            ttime = ttime + time
            vd =  sf/real(Vgs)     ! Estimate for this group
            sigma = sqrt((sf2/Vgs - vd*vd)/Vgs)
            tinc = tinc + inc   ! Total so far in circle
            tsf  = tsf + sf     ! SUM F so far
            tsf2 = tsf2 + sf2   ! SUM F^2 so far
            tvd = tsf/real(k*Vgs)  ! Estimate so far
            kVgs = k*Vgs
            tsigma = sqrt((tsf2/kVgs - tvd*tvd)/kVgs)
            vd = twod*vd
            sigma = twod*sigma
            tvd = twod*tvd
            tsigma = twod*tsigma
            print 21, k*gs, inc, vd, sigma, tvd, tsigma, ttime/1000.0
 21         FORMAT (1X, I7, 1X, I8, 2(1X, F7.5, 1X, E11.5), 1X, F9.5)
        end do
        print 22, tinc, 100.0*tinc/(ng*Vgs)
 22     FORMAT (" Total hits = ", I7, " = ", F5.1, "%")
      end do
 99   continue   ! For exit
      end program montepimprt

