      program rantest
      parameter(NREP=50,ORDER=10,MAX3=50000,NCELLS=10)
      parameter(MAX2=100000,MAX1=100000,BMAX=10)
      integer cell(0:9,0:9,0:9)
      dimension x(0:10),g(0:10),gave(0:10),gave2(0:10)
      integer bin(0:9),seed,outbounds
      ave = 0.0
      ave2 = 0.0
      write(6,1)
    1 format(' seed ?')
      read(5,*) seed
      write(6,2) seed
    2 format('seed = ',i20)
      call ranset(seed) 
c
c/* test (1) */
c
      do 10 i=0,BMAX-1
   10 bin(i) = 0
      do 20 i=0,MAX1-1
      y = ranf(seed)
      ave = ave+y
      if(y.lt.0.or.y.gt.1) outbounds=outbounds+1
      do 30 j=0,BMAX-1
         if(y.ge.float(j)/BMAX.and.y.lt.float(j+1)/BMAX)
     c      bin(j)=bin(j)+1
   30 continue
   20 continue
      i = MAX1
      write(6,3) i
    3 format(' total number of random numbers generated = ',i10)
      write(6,4) ave/float(i),outbounds
    4 format(' average = ',f8.5,', number out of bounds = ',i10/)
      write(6,5)
    5 format(5x,'bin #',5x,'occupation'/)
      ave = 0.0
      do 40 j=0,BMAX-1
      write(6,6) j,bin(j)
    6 format(i10,5x,i10)
      ave = ave+(bin(j)-(MAX1+0.)/BMAX)**2
   40 continue
      i1 = BMAX - 1
      write(6,7) ave/((MAX1+0.)/BMAX),i1
    7 format(' chi-sqr = ',1p,e12.5,' for ',i8,' degrees of freedom'//)
c
c/* test (2) */
c
      do 50 i1=0,NCELLS-1
      do 50 i2=0,NCELLS-1
      do 50 i3=0,NCELLS-1
   50 cell(i1,i2,i3) = 0
      do 60 i=1,MAX2
      i1 = ranf(seed)*NCELLS
      i2 = ranf(seed)*NCELLS
      i3 = ranf(seed)*NCELLS
      cell(i1,i2,i3)=cell(i1,i2,i3)+1
   60 continue
      ave = MAX2/float(NCELLS*NCELLS*NCELLS)
      phi = 0.0
      do 70 i1=0,NCELLS-1
      do 70 i2=0,NCELLS-1
      do 70 i3=0,NCELLS-1
   70 phi = phi+(ave - float(cell(i1,i2,i3)))**2
      phi = phi/MAX2
      i = MAX2
      write(6,8) i,i1,i2,i3
    8 format(' #triples produced =',i10,', split into a space',
     c' of ',i4,'x',i4,'x',i4,' cells')
      write(6,9) ave,phi
    9 format(' average number per cell = ',1p,e12.5,
     c', effective chi-sqr = ',e12.5)
      write(6,11)
   11 format(' The value should be near 1'//)
c
c/* test 3 */
c
      do 80 i=0,ORDER
      g(i) = 0.0
      gave(i) = 0.0
   80 gave2(i) = 0.0
      do 90 i=1,ORDER
   90 x(i) = ranf(seed)
      i3 = ORDER
c
c
c
      do 100 j=1,NREP
c
c
c
      ave = 0.0
      ave2 = 0.0
      do 110 i=1,MAX3
	xnew = ranf(seed)
	ave = ave+xnew	
	ave2 = ave2+xnew*xnew
	i3=i3+1
        do 120 i1=1,ORDER
	   i2 = i3 - i1 
	   if(i2 .lt. 1) i2 = i2+ORDER
	   g(i1) = g(i1)+xnew*x(i2)
  120   continue
        if(i3 .gt. ORDER) i3 = 1
        x(i3) = xnew
  110 continue
      ave = ave/MAX3
      g(0) = ave2
      do 130 i=0,ORDER
	g(i) = g(i)/MAX3
	gave(i) = gave(i)+g(i)
	gave2(i) = gave2(i)+g(i)**2
	g(i) = 0.0
  130 continue
c
c
c
  100 continue
c
c
c
      ave = 0.0
      do 140 i1=0,ORDER
	gave(i1) = gave(i1)/NREP
	gave2(i1) = gave2(i1)/NREP
	gave2(i1) = gave2(i1)-gave(i1)**2
	gave2(i1) = sqrt(gave2(i1))
	ave = ave+gave(i1)
  140 continue
      ave = ave-gave(0)
      ave = ave/ORDER
      j = NREP
      i = MAX3
      write(6,12) j,i
   12 format(' for ',i6,' blocks of ',i8,' numbers ',
     c'the autocorrelation coefficients are:'/)
      write(6,13)
   13 format(5x,'order    coefficient   variance')
      do 150 i=0,ORDER
         write(6,14) i,gave(i),gave2(i)
   14    format(5x,i5,3x,1p,e12.5,3x,e12.5)
  150 continue
      write(6,15) ave
   15 format(/' average coefficient = ',1p,e12.5)
      write(6,16) 0.3/sqrt(MAX3+0.)
   16 format(' theoretically g(0) = 1/3 and the rest = 1/4',
     c'theoretical variance = ',1p,e12.5)
c
c
c
      stop
      end
c
c <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
c              Function "ranf" for random number generation
c >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c
      real function ranf(seed)
      integer seed
      seed=seed*1566083941+1
      ranf=float(seed)*2.328306e-10+0.5
      return
      end
      subroutine ranset(seed)
      integer seed
      return
      end
