C mpfortran -pesize=8 -d_lines mpfrand.o randtest.f
C
C----ROB SALGADO (cps615g)
C    randtest.f 09/27/92 (09/28/92 PROBLEM 2.2)
C
C    TEST THE RANDOM NUMER GENERATOR
C
      program randtest
      integer iseed,n
      real psum,pavg
      real x(2000000),sum,mean,sumdevsq
      integer p(20,20)

      iseed=0

      open (unit=3,file='rand-test')
      write(3,*) 'Random-Number-Generator Test [ROB SALGADO]'
      print*, 'Random-Number-Generator Test [ROB SALGADO]'

      print *, 'Number of samples (even number up to 2,000,000)'
      read*, n
      n=int(n/2)*2
      print*, n
C------------------------------------PAIRWISE COLLECT RANDOM NUMBERS
      sum=0.0
      do 100 j=1,n,2
         x(j)=rand(iseed)
         x(j+1)=rand(iseed)
C----SAVE SUM FOR STATISTICS
         sum=sum+x(j)+x(j+1)
C----BIN PAIRWISE IN 20x20 ARRAY
         k1=int(20*x(j)+1)
         k2=int(20*x(j+1)+1)
         p(k1,k2)=p(k1,k2)+1
 100  continue
C------------------------------------------------------------------|

C-------------------------------------------CRUDE 20x20 PLOT
      write(3,*) ' '
      write(3,*) 'ABSOLUTE 20x20 PLOT'
      write(3,*) ' '
      do 150 j=1,20,1
         write(3,999) (p(i,21-j), i=1,20)
 150  continue

      write(3,*) ' '
      write(3,*) '^ x(j+1)'
      write(3,*) '|'
      write(3,*) '|'
      write(3,*) ' ----> x(j)'
      write(3,*) ' '
 
C-------------------------------------------CRUDE PLOT WRT ARRAY AVERAGE
      write(3,*) '20x20 PLOT RELATIVE TO ARRAY AVERAGE'
      write(3,*) ' '
      psum=0
      do 160 j=1,20
         do 161 i=1,20
            psum=psum+p(i,j)
 161     continue
 160  continue
      pavg=psum/(20*20)
      do 170 j=1,20
         write(3,998) (p(i,21-j)/pavg, i=1,20)
 170  continue
C----------------------------------------------------------|

C--------------------------STATISTICS
      mean=sum/n
      write(3,*) ' '
      write(3,*) n, ' Total Random Numbers'
      write(3,*) 'Mean=   ',mean

      sumdevsq=0.0
      do 200 j=1,n,1
         sumdevsq=sumdevsq+(mean-x(j))*(mean-x(j))
 200  continue
      write(3,*) 'Std Dev=',sqrt(sumdevsq/n)
 998  format(20F4.1)
 999  format(20I4)
      end





