! Last change: DOS 31 Jul 2000 8:18 pm ! *** copyright 2000 *** ! *** filename ranbf.f95 *** John F. Monahan ** ! ********************** PROGRAM PRANBF ! BASIC CHI-SQUARE TEST OF UNIFORM PSEUDORANDOM GENERATOR IMPLICIT NONE REAL CHISQ,U,EXP,FKCELL,RANBF 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 RANBF UNIFORM GENERATOR'/ & & ' WITH',I8,' OBS IN',I4,' CELLS GIVES CHISQ =',F12.4) ! OPEN( UNIT=6, FILE='ranbf.out' ) ! ! INITIALIZE CELL COUNTERS TO ZERO CELL = 0 ! INITIALIZE GENERATOR EXP = RANBF(5151917) ! EXPECTED NUMBER IN EACH CELL FKCELL = REAL(KCELLS) EXP = REAL(N)/FKCELL ! FILL UP THE CELLS DO I = 1,N U = RANBF(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 PRANBF REAL FUNCTION RANBF(IDUM) ! PORTABLE IMPLEMENTATION OF UNIFORM PSEUDORANDOM NUMBER GENERATOR ! LEWIS, GOODMAN, & MILLER MULTIPLICATIVE GENERATOR ! X(N+1) = MOD( 16807*X(N), 2**31-1 ) ! ! P. BRANTLEY, B.L. FOX, L. SCHRAGE (1983) A GUIDE TO SIMULATION ! SPRINGER-VERLAG, NEW YORK. PP. 200-202. ! UPDATED VERSION OF ! LINUS SCHRAGE (1979)'A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR' ! ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 5, PP. 132-138. ! ! ARGUMENT ! IDUM INTEGER FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS ! IMPLICIT NONE INTEGER, INTENT(IN) :: IDUM REAL, PARAMETER :: RP = 4.656612875E-10 ! 1/P INTEGER, PARAMETER :: A = 16807 ! MULTIPLIER INTEGER, PARAMETER :: B = 127773 ! B = P / A INTEGER, PARAMETER :: C = 2836 ! C = P MOD A INTEGER, PARAMETER :: P = 2147483647 ! MODULUS 2**31 - 1 INTEGER, SAVE :: IX = 0 INTEGER K1 ! ! IF NOT FIRST CALL, THEN SKIP SETTING SEED IF( IX .EQ. 0) IX = IDUM ! WRITE NUMBER AS ALPHA*2**16 + BETA K1 = IX / B IX = A*( IX - K1*B) - K1*C ! ABOVE DOES A*IX MOD B -K1*C IF( IX .LT. 0 ) IX = IX + P ! RP BELOW IS RECIPROCAL OF P RANBF = REAL(IX)*RP RETURN END FUNCTION RANBF ! *** end of filename ranbf.f95 *********************