! This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/byterand.f90
!
!=====Sample driver for portable pseudo-random number generator, RAND.
!     Purpose: Portable quasi-random number generator package with
!              good statistical properties that utilizes only 16-bit
!              integers (good for mainframe and microcomputers).
!     Some Properties: - Portable Fortran
!                      - uses only 16 bit integer constants
!                      - very long sequence; good statistical properties
!                      - random repeats (a rare property)
!                      - easily adaptable to other languages
!                      - well-documented in qualified literature.
!     H. D. Knoble, August 1987, hdkLESS at SPAM psu dot edu
!     Penn State University Center, for Academic Computing
!
!  Author's statement of intent: Program units following this notice
!  may be freely copied by anyone under the condition that they be used
!  only for peaceful purposes that do not threaten human dignity and
!  life or the environment.
!
!---This driver program is used to illustrate how to use the subprograms
!   RAND, RPERM, and IRAND by interactively generating quasi-random
!   permutations of the first N integers.
!
!   To quit this sample driver program simply press Enter
!   in response to any prompt.
!
      INTEGER X(999),N
!---Maximum value of N is dimension of X.
      INTEGER SEED(3)
      CHARACTER*80 NC
      COMMON/RANDOM/SEED
!
      WRITE(*,*) '** Random Permutation of the first N Integers **'
      WRITE(*,*) 'Please enter the value of N (or press Enter to quit):'
      NC=' '
      READ(*,1) NC
 1    FORMAT(A)
      IF (NC.EQ.' ') GOTO 5
      READ(NC,*,END=5) N
      IF (N.EQ.0.OR.N.GT.999) GOTO 5
      DO NSETS=1,9999
!---Generate a random permutation of the 1ST N integers.
!!--Generate a "random" seed using Time-Of-Day clock (VS Fortran only)
!*      DOUBLE PRECISION TOD
!*      INTEGER ISEED
!*!---TOD is returned as the Time Of Day in 1/1000000ths of a second.
!*    VS Fortran only.
!*      CALL CLOCKX(TOD)
!*      ISEED=MAX(1+MOD(INT(TOD),1000)/10, 1+MOD(INT(TOD),100))
        WRITE(*,*) 'Please enter 3 integer seeds, each >0 & <30000:'
        NC=' '
        READ(*,1) NC
        IF(NC.EQ.' ') GOTO 5
        READ(NC,*,END=5) SEED
        IF(SEED(1).EQ.0) GOTO 5
!---Generate a random permutation, X
        CALL RPERM(X,N)
!---Output the random permutation.
        WRITE(*,2) N
 2      FORMAT(' Random permutation of the 1st ',I4,' integers:')
        WRITE(*,4) (X(I),I=1,N)
 4      FORMAT(20I4)
      END DO
 5    STOP
      END
      SUBROUTINE RPERM(P,N)
!===Generate a random permutation, P, of the first N integers.
!   (Shuffling, equivalent to sampling WITHOUT REPLACEMENT).
!   Adaptation of Knuth Volume 2, Algorithm 3.4.2P.
      INTEGER N,P(N), K,J,I,IPJ,ITEMP,M
      REAL U(100)
      DO I=1,N
        P(I)=I
      END DO
!---Generate up to 100 U(0,1) numbers at a time.
      DO I=1,N,100
        M=MIN(N-I+1,100)
        CALL RAND(U,M)
        DO J=1,M
          IPJ=I+J-1
          K=INT(U(J)*(N-IPJ+1))+IPJ
          ITEMP=P(IPJ)
          P(IPJ)=P(K)
          P(K)=ITEMP
        END DO
      END DO
      RETURN
      END
      SUBROUTINE IRAND(S,N,LOW,HI)
!===Generate a random integer sequence: S(1),S(2), ... ,S(N)
!   such that each element is in the closed interval <LOW,HI> and
!   sampled WITH REPLACEMENT.                            HDK, JUNE 1971.
      INTEGER N,S(N),LOW,HI,IX,I
      REAL U(1)
      DOUBLE PRECISION X
      DO I=1,N
        CALL RAND(U,1)
!---Use DP arithmetic to effect a more precise transformation.
        X=DBLE((HI+1)-LOW)*U(1) + DBLE(LOW)
        IX=X
        IF(X.LT.0 .AND. IX.NE.X) IX=X-1.D0
        S(I)=IX
      END DO
      RETURN
      END
      BLOCK DATA
!=======================================================================
!  Portable pseudo-random integer generator, especially for
!  microcomputers with a limitation of 16 bit integers. Translated from
!  original Pascal version(1) to Fortran 77 by H. D. Knoble, PSU.
!
!  Corrected typo (change 30308.0 to 30307.0) courtesy of Ed Wynn
!  and Dave Dodson, 3 May 2002.
!
!   The supporting paper is:
!   (1) B. Wichmann & D. Hill, "Building a Random-Number Generator",
!             BYTE, March, 1987, 12(3):127-128.
!
!   Also see the following related works:
!   (2) Wichmann, B.A. and Hill, I.D., "An Efficient and Portable",
!             Pseudo-random Number Generator", Applied Statistics,
!             Volume 31, 1982, pages 188-190.
!   (3) Haas, Alexander, "The Multiple Prime Random Number Generator",
!             ACM Transactions on Mathematical Software; December,
!             1987, 13(4):368-381.
!   (4) L'Ecuyer, Pierre, "Efficient and Portable Combined Random Number
!             Generators", Communications of the ACM; June, 1988,
!             31(6):742-749,774.
!
! Use...
!      CALL RAND(U,N)
!          To generate a sequence, U, of N Uniform(0,1) numbers.
!          Cycle length is ((30269-1)*(30307-1)*(30323-1))/4  or
!          6953607871644  > 6.95E+12.
!
!     To access the SEED vector in the calling program use statements:
!     INTEGER SEED(3)
!     COMMON/RANDOM/SEED
!
!  The common variable SEED is the array of three current seeds.
      INTEGER SEED(3)
      COMMON/RANDOM/SEED
      DATA SEED(1),SEED(2),SEED(3)/1,10000,3000/
      END
!=======================================================================
      SUBROUTINE RAND(U,N)
      INTEGER N,X,Y,Z
      REAL U(N),V
      COMMON/RANDOM/X,Y,Z
      IF(N.LE.0) RETURN
      DO  I=1,N
        X=171*MOD(X,177)-2*(X/177)
        IF(X.LT.0) X=X+30269
        Y=172*MOD(Y,176)-35*(Y/176)
        IF(Y.LT.0) Y=Y+30307
        Z=170*MOD(Z,178)-63*(Z/178)
        IF(Z.LT.0) Z=Z+30323
        V=X/30269.0 + Y/30307.0 + Z/30323.0
        U(I)=V-INT(V)
      END DO
      RETURN
      END