!     Last change:  DOS  31 Jul 2000    8:13 pm
!                                       *** copyright 2000 ***
! *** filename ran50.f95                *** John F. Monahan **
!                                       **********************
      PROGRAM PRAN50
!  BASIC CHI-SQUARE TEST OF UNIFORM PSEUDORANDOM GENERATOR
      IMPLICIT NONE
      REAL CHISQ,U,EXP,FKCELL,RAN50
      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 RAN50 UNIFORM GENERATOR'/ &
     &    ' WITH',I8,' OBS IN',I4,' CELLS GIVES CHISQ =',F12.4)
!
      OPEN( UNIT=6, FILE='ran50.out' )
!
!                   INITIALIZE CELL COUNTERS TO ZERO
      CELL = 0
!                   INITIALIZE GENERATOR
      EXP = RAN50(5151917)
!                   EXPECTED NUMBER IN EACH CELL
      FKCELL = REAL(KCELLS)
      EXP = REAL(N)/FKCELL
!                   FILL UP THE CELLS                                                    
      DO I = 1,N
      U = RAN50(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 PRAN50
      REAL FUNCTION RAN50(IDUM)
!  UNIFORM PSEUDORANDOM NUMBER GENERATOR
!   A LINEAR CONGRUENTIAL GENERATOR X(N+1) = MOD( A*X(N)+C, 2**48 )
!
!  ARGUMENT
!   IDUM    FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS
!
!  NUMBERS WRITTEN HERE IN FOUR PIECES,
!     X = X0 + (2**12)*X1 + (2**24)*X2 + (2**36)*X3
!     A = A0 + (2**12)*A1 + (2**24)*A2 + (2**36)*A3
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: IDUM
      INTEGER, PARAMETER :: A0 = 3781             ! MULTIPLIER
      INTEGER, PARAMETER :: A1 = 3010             ! MULTIPLIER
      INTEGER, PARAMETER :: A2 = 418              ! MULTIPLIER
      INTEGER, PARAMETER :: A3 = 11               ! MULTIPLIER
      INTEGER, PARAMETER :: C = 1                 ! CONSTANT
      INTEGER, PARAMETER :: TWO12 = 4096          ! 2**12
!                   INITIAL VALUE OF ZERO TO ALLOW SETTING SEED
!                    WITH FIRST CALL
      INTEGER, SAVE :: X0 = 0
      INTEGER, SAVE :: X1 = 0
      INTEGER, SAVE :: X2 = 0
      INTEGER, SAVE :: X3 = 0
      INTEGER CARRY,P0,P1,P2,P3
      REAL, PARAMETER :: FTWO12 = 4096.
!
!                   MULTIPLIER IS A = B1A2BC2EC5(HEX) = 5**17
!                     A3 =   B(HEX) =   11 (TEN)
!                     A2 = 1A2(HEX) =  418 (TEN)
!                     A1 = BC2(HEX) = 3010 (TEN)
!                     A0 = EC5(HEX) = 3781 (TEN)
!                   CONSTANT IS C = 1 (HEX) = 1 (TEN)
!
!
!                   IF NOT FIRST CALL, THEN SKIP SETTING SEED
      IF( X0 + X1 + X2 + X3 .EQ. 0 ) THEN
      X1 = IDUM / TWO12
      X0 = IDUM - X1*TWO12
                                     END IF ! ( FIRST CALL )
!                   FIRST GET ALL THE PRODUCTS
      P0 = X0*A0 + C
      P1 = X0*A1 + X1*A0
      P2 = X0*A2 + X1*A1 + X2*A0
      P3 = X0*A3 + X1*A2 + X2*A1 + X3*A0
!                   NOW DO MOD 2**12 AND CARRY TO NEXT DIGIT
      CARRY = P0 / TWO12
      X0 = P0 - CARRY*TWO12
      P1 = P1 + CARRY 
      CARRY = P1 / TWO12
      X1 = P1 - CARRY*TWO12
      P2 = P2 + CARRY 
      CARRY = P2 / TWO12
      X2 = P2 - CARRY*TWO12
      P3 = P3 + CARRY 
      CARRY = P3 / TWO12
      X3 = P3 - CARRY*TWO12
!                   NOW FORM UNIFORM FROM THE FOUR PIECES
      RAN50 = ( REAL(X3) + ( REAL(X2) + ( REAL(X1) +     &
     &           REAL(X0) / FTWO12 ) / FTWO12 ) / FTWO12 ) / FTWO12
      RETURN
      END FUNCTION RAN50
! *** end of filename ran50.f95         *********************