!     Last change:  DOS  31 Jul 2000    8:15 pm
!                                       *** copyright 2000 ***
! *** filename ran69.f95                *** John F. Monahan **
!                                       **********************
      PROGRAM PRAN69
!  BASIC CHI-SQUARE TEST OF UNIFORM PSEUDORANDOM GENERATOR
      IMPLICIT NONE
      REAL CHISQ,U,EXP,FKCELL,RAN69
      INTEGER, DIMENSION(64) :: CELL
!                   DO 2**10 = 1024 UNIFORMS IN 64 CELLS
      INTEGER, PARAMETER :: N = 1024
      INTEGER, PARAMETER :: KCELLS = 64
      INTEGER I,J
!
  21  FORMAT(' RESULT OF CHI-SQUARE TEST ON RAN69 UNIFORM GENERATOR'/ &
     &    ' WITH',I8,' OBS IN',I4,' CELLS GIVES CHISQ =',F12.4)
!
      OPEN( UNIT=6, FILE='ran69.out' )
!
!                   INITIALIZE CELL COUNTERS TO ZERO
      CELL = 0
!                   INITIALIZE GENERATOR
      EXP = RAN69(5151917)
!                   EXPECTED NUMBER IN EACH CELL
      FKCELL = REAL(KCELLS)
      EXP = REAL(N)/FKCELL
!                   FILL UP THE CELLS                                                    
      DO I = 1,N
      U = RAN69(I)                                                      
!                   J POINTS TO CELL
      J = 1 + INT( U*FKCELL )
!                   INCREMENT COUNT IN CELL
      CELL(J) = CELL(J) + 1
      END DO  ! LOOP ON I
!                   TEST FOR UNIFORMITY
      CHISQ = 0.
      DO J = 1,KCELLS
      CHISQ = CHISQ + (CELL(J) - EXP)*(CELL(J) - EXP)/EXP
      END DO  ! LOOP ON J
!                   WRITE OUT RESULTS
      WRITE(6,21) N,KCELLS,CHISQ
      STOP                                                              
      END PROGRAM PRAN69
      REAL FUNCTION RAN69(IDUM)
!  UNIFORM PSEUDORANDOM NUMBER GENERATOR
!   A LINEAR CONGRUENTIAL GENERATOR X(N+1) = MOD( 69069*X(N), 2**31)
!
!  NUMBERS WRITTEN HERE IN TWO PIECES, H*(2**15) + L
!  ARGUMENT 
!   IDUM    INTEGER  FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: IDUM
      INTEGER, PARAMETER :: B15 = 32768           ! 2**15
      INTEGER, PARAMETER :: B16 = 65536           ! 2**16
      INTEGER, SAVE :: H = 0
      INTEGER, SAVE :: L = 0
!
!                   IF NOT FIRST CALL, THEN SKIP SETTING SEED
      IF( (H .EQ. 0) .AND. (L .EQ. 0) ) THEN
      H = IDUM / B15
      L = IDUM - H*B15
                                        END IF ! ( FIRST CALL )
!                   MULTIPLIER IS 69069 = 2**16 + 3533
      H = 3533*H + 2*L                                                     
      L = 3533*L                                                           
!                   MODULUS IS 2**31
      H = MOD( H + (L/B15) , B16)
      L = MOD( L , B15 )
!                                     2**15      2**16
      RAN69 = ( REAL(H) + ( REAL(L) / 32768. ) )/65536.
      RETURN
      END FUNCTION RAN69
! *** end of filename ran69.f95         *********************