!     Last change:  DOS  31 Jul 2000    8:30 pm
!                                       *** copyright 2000 ***
! *** filename ranft.f95                *** John F. Monahan **
!                                       **********************
!
!  SOME COMPILERS REQUIRE MODULE TO PREVIOUSLY EXIST OR BE
!   AT THE BEGINNING
      MODULE GFSRFT
      INTEGER TABPTJ
      INTEGER, DIMENSION(521) :: TABLE
      END MODULE GFSRFT
!
      PROGRAM PRANFT
!  BASIC CHI-SQUARE TEST OF UNIFORM PSEUDORANDOM GENERATOR
      USE GFSRFT
      IMPLICIT NONE
      REAL CHISQ,U,EXP,FKCELL,RANFT
      INTEGER, DIMENSION(64) :: CELL
!                   DO 2**10 = 1024 UNIFORMS IN 64 CELLS
      INTEGER, PARAMETER :: N = 1024
      INTEGER, PARAMETER :: KCELLS = 64
      INTEGER I,J
!
  20  FORMAT(6I12)
  21  FORMAT(' RESULT OF CHI-SQUARE TEST ON RANFT UNIFORM GENERATOR'/ &
     &    ' WITH',I8,' OBS IN',I4,' CELLS GIVES CHISQ =',F12.4)
!
      OPEN( UNIT=6, FILE='ranft.out' )
      OPEN( UNIT=5, FILE='gfsrft.tb0' )            ! original table
      OPEN( UNIT=7, FILE='gfsrft.tb1' )            ! store table at end
!
!                   INITIALIZE CELL COUNTERS TO ZERO
      CELL = 0
!                   INITIALIZE GENERATOR WITH TABLE
      READ(5,20) TABPTJ,TABLE
!                   EXPECTED NUMBER IN EACH CELL
      FKCELL = REAL(KCELLS)
      EXP = REAL(N)/FKCELL
!                   FILL UP THE CELLS                                                    
      DO I = 1,N
      U = RANFT(I)                                                      
!                   J POINTS TO CELL
      J = 1 + INT( U*FKCELL )
!                   INCREMENT COUNT IN CELL
      CELL(J) = CELL(J) + 1
      END DO  ! LOOP ON I
!                   SAVE TABLE
      WRITE(7,20) TABPTJ,TABLE
!                   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 PRANFT
      REAL FUNCTION RANFT(IDUM)
!   FUSHIMI-TEZUKA GFSR UNIFORM RANDOM NUMBER GENERATOR
!
!  USES PRIMITIVE TRINOMIAL WITH P=521 AND Q=32 AS USED BY BRIGHT &
!   ENISON AND ARVILLIAS & MARITSAS BUT WITH RANDOM SEED MATRIX
!  FUSHIMI & TEZUKA GIVE RULES FOR TESTING K-DISTRIBUTION OF
!   SEQUENCE -- THE ORIGINAL SEED TABLE HAS BEEN CHECKED AND
!   31 BIT NUMBERS ARE 16-DISTRIBUTED (BEST POSSIBLE)
!
!  M. FUSHIMI & S. TEZUKA (1983) THE K-DISTRIBUTION OF GENERALIZED
!   FEEDBACK SHIFT REGISTER PSEUDORANDOM NUMBERS, COMMUNICATIONS OF
!   THE ACM, VOLUME 26, NUMBER 7, PP. 516-523
!
!  ARGUMENT IS A DUMMY AND NEVER USED
!
!  INPUT POINTER AND TABLE PASSED IN MODULE GFSRFT
!  ORIGINAL TABLE IN THE FILE 'gfsrft.tb0' IN 6I12 FORMAT
!  CAN WRITE POINTER AND TABLE IN FILE AS SEED FOR NEXT EXPERIMENT
!
      USE GFSRFT              ! MODULE HOLDS TABLE AND POINTER
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: IDUM
!                   PARAMETERS OF TAUSWORTHE SEQUENCE
      INTEGER, PARAMETER :: P = 521
      INTEGER, PARAMETER :: Q = 32
      INTEGER K
      REAL, PARAMETER :: FN = 2147483648.         ! 2**31
!
!                   UPDATE POINTER
      TABPTJ = TABPTJ + 1                                                   
      IF ( TABPTJ .GT. P ) TABPTJ = 1
!                   UPDATE DELAY POINTER
      K = TABPTJ + Q                                                        
      IF ( K .GT. P ) K = K - P                                                 
!                   COMPUTE EXCLUSIVE OR OF TWO TABLE ENTRIES
!                    AND REPLACE WITH NEW ONE
      TABLE(TABPTJ) = IEOR( TABLE(K), TABLE(TABPTJ) )                                                   
!                   CONVERT BIG INTEGER TO FLOATING POINT NUMBER
      RANFT = REAL( TABLE(TABPTJ) ) / FN
      RETURN
      END FUNCTION RANFT
! *** end of filename ranft.f95         *********************