C----ROB SALGADO (cps615g)
C    1bd.f  10/03/92 (SET 3 [??] PROBLEM 3)
C
C    CALCULATE d-VOLUMES WITH IMPORTANCE SAMPLING


      program impdspheremonte
      integer iseed,n,r,dim
      real radius,x,l,t,logv,truevol,pi
      real hit,sum,mean,sumdevsq,monterad(1000)
      real oldhit,oldsum,oldmean,oldsumdevsq,oldmonterad(1000)
      pi=3.1415926535

      open (unit=3, file='d-sphere')
      write(3,*) 'Importance d-Sphere Monte-Carlo [ROB SALGADO]'
      print*, 'Importance d-Sphere Monte-Carlo [ROB SALGADO]'
      
      iseed=0

      print*, 'Dimension dim=?'
      read*, dim
      print*, 'Radius radius=?'
      read*, radius
C-----------------------------COMPUTE TRUE D-VOLUME
C     Use logarithms because partial factors for high dimensions
C     overflows the variables.
      if (2*int(dim/2).eq.dim) goto 8
C-------DIM IS ODD
      logv=(dim+1.0)/2.0*log(2.0)+(dim-1.0)/2.0*log(pi)
      goto 9
C-------DIM IS EVEN
 8    logv=(dim/2.0)*log(2.0*pi)
 9    logv=logv+(dim)*log(radius)
      do 5 i=dim,1,-2
           logv=logv-log(i*1.0)
 5    continue     
      truevol=exp(logv)
C-------------------------------------------------|
      print*, 'Number of samples per run?'
      read*, n
 10   print*, 'Number of runs?'
      read*, r
      if (r.gt.1000) goto 10

C------------------------------------------------------r RUNS
      sum=0.0
      oldsum=0.0
      write(3,*) 'Trial     Radius'
      do 200 j=1,r,1
         hit=0
         oldhit=0
C---------------------------------n SAMPLES
      do 100 i=1,n,1
         l=0.0
         t=0.0
         do 50 k=1,dim,1
            x=rand(iseed)
C------search for largest rectangular coord
            if (x.gt.t) t=x
            l=l+x*x
 50      continue

C------weight the randompoint in the unit d-cube
C      by (t/r)**dim, where t=max(x1,x2,x3,..,xdim) for d-dimensions
C                     where r=sqrt(l)=sqrt(x1^2+x2^2+...+xdim^2)
         hit=hit+(sqrt(t*t/l))**dim

C------use the same random numbers to Monte-Carlo the naive way
C
         if (l.le.1.0) oldhit=oldhit+1
 100  continue
C------------------------------------------|
      monterad(j)=(2*radius)**dim*hit/n
      oldmonterad(j)=(2*radius)**dim*oldhit/n
      write(3,*) j, monterad(j),oldmonterad(j)
      sum=sum+monterad(j)
      oldsum=oldsum+oldmonterad(j)
 200  continue
C-----------------------------------------------------------|

C------------------------------STATISTICS
      mean=sum/r
      oldmean=oldsum/r
      write(3,*) '----'
      write(3,*) r,' Runs, each with',n,' samples'
      write(3,*) 'Spherical volume in ',dim,'-dimensions; r=',radius
      sumdevsq=0.0
      oldsumdevsq=0.0

      do 300 j=1,r,1
         sumdevsq=sumdevsq+(mean-monterad(j))**2
         oldsumdevsq=oldsumdevsq+(oldmean-oldmonterad(j))**2
 300  continue

        
      write(3,*) oldmean,'=NAIVE      ',sqrt(oldsumdevsq/n),' std dev'
      write(3,*) mean,'=IMPORTANCE ',sqrt(sumdevsq/n),' std dev'
      write(3,*) truevol,'=TRUE volume'
      write(3,*) '% Naive Error     =', 100.0*(oldmean-truevol)/truevol
      write(3,*) '% Importance Error=', 100.0*(mean-truevol)/truevol
 
      end
