!     Last change:  DOS  31 Jul 2000    8:23 pm
!                                       *** copyright 2000 ***
! *** filename ranf1.f95                *** John F. Monahan **
!                                       **********************
      PROGRAM PRANF1
!  BASIC CHI-SQUARE TEST OF UNIFORM PSEUDORANDOM GENERATOR
      IMPLICIT NONE
      REAL CHISQ,U,EXP,FKCELL,RANF1
      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 RANF1 UNIFORM GENERATOR'/ &
     &    ' WITH',I8,' OBS IN',I4,' CELLS GIVES CHISQ =',F12.4)
!
      OPEN( UNIT=6, FILE='ranf1.out' )
!
!                   INITIALIZE CELL COUNTERS TO ZERO
      CELL = 0
!                   INITIALIZE GENERATOR
      EXP = RANF1(5151917)
!                   EXPECTED NUMBER IN EACH CELL
      FKCELL = REAL(KCELLS)
      EXP = REAL(N)/FKCELL
!                   FILL UP THE CELLS                                                    
      DO I = 1,N
      U = RANF1(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 PRANF1
      REAL FUNCTION RANF1(IDUM)
!  UNIFORM PSEUDORANDOM NUMBER GENERATOR
!   A LINEAR CONGRUENTIAL GENERATOR X(N+1) = MOD( A*X(N), 2**31 - 1 )
!  CODING FOLLOWS ALGORITHM 1 OF HORMANN AND DERFLINGER, WITH
!   PERSONAL MODIFICATIONS TO AVOID OVERFLOWS
!
!  W. HORMANN & G. DERFLINGER (1993) "A PORTABLE RANDOM NUMBER GENERATOR
!   WELL SUITED FOR THE REJECTION METHOD," ACM TOMS VOL 19, PP.489-495.
!
!  MULTIPLIER FROM
!  G. FISHMAN & L. MOORE (1982) "A STATISTICAL EVALUATION OF MULITPLICA-
!   TIVE CONGRUENTIAL RANDOM NUMBER GENERATORS WITH MODULUS 2**32-1,"
!   J AMERICAN STATISTICAL ASSN, VOL 77, PP. 129-136.
!
!  ARGUMENT
!   IDUM    FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS
!
!  NUMBERS WRITTEN HERE IN TWO PIECES, X = XHI*(2**16) + XLOW
!     AND A = AHI*(2**15) + ALOW          ***** A MUST BE < 2**30 *****
!     2**15 = 32768 AND 2**16 = 65536
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: IDUM
      INTEGER, PARAMETER :: AHI = 19237          ! PART OF MULTIPLIER
      INTEGER, PARAMETER :: ALOW = 2000          ! PART OF MULTIPLIER
      INTEGER, PARAMETER :: ALOW2 = 4000         ! PART OF MULTIPLIER
      INTEGER, PARAMETER :: B15 = 32768           ! 2**15
      INTEGER, PARAMETER :: B16 = 65536           ! 2**16
      INTEGER, PARAMETER :: P = 2147483647        ! MODULUS 2**31 - 1
      INTEGER, SAVE :: XHI = 0
      INTEGER, SAVE :: XLOW = 0
      INTEGER MID1,MID2,MID,X,K
!
!                   IF NOT FIRST CALL, THEN SKIP SETTING SEED
      IF( (XHI .EQ. 0) .AND. (XLOW .EQ. 0) ) THEN
      XHI = IDUM / 65536
      XLOW = IDUM - XHI*65536
                                             END IF ! ( FIRST CALL )
!                   MULTIPLIER IS A = 630360016 = AHI*(2**15) + ALOW
      MID1 = AHI*XLOW
      MID2 = ALOW2*XHI                                                     
!                   TEST FOR OVERFLOW
      IF( MID1-1 .GT. P-MID2 ) THEN
!                   HERE MID IS > 2**31, SO WRITE AS MID - 2**31
      MID = (MID1-1) - (P-MID2)
      K = 1
                               ELSE
!                   HERE MID IS < 2**31, SO WRITE AS POSITIVE
      MID = MID1 + MID2
      K = 0
                               END IF ! ( MID1-1 .GT. P-MID2 )
!                   NOW GET TO MAIN STEP
!                    SUBTRACT P = 2**31 - 1 IN ADVANCE
      MID2 = MID / B16
      X = (ALOW*XLOW - P) + AHI*XHI + MID2 + K*B15
      IF( X .LT. 0 ) X = X + P
      X = X + ( ( MID - MID2*B16 )*B15 - P )
      IF( X .LT. 0 ) X = X + P
!                   WILL NEED TWO PARTS OF X FOR NEXT CALL
      XHI = X / B16
      XLOW = X - XHI*B16
!                                          2**16      2**15
      RANF1 = ( REAL(XHI) + ( REAL(XLOW) / 65536. ) )/32768.
      RETURN
      END FUNCTION RANF1
! *** end of filename ranf1.f95         *********************