C This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/knuths.for
C Purpose: Select K records from sequential Input File at random
C          from a set of N records,  where 0 < K <= N. Write
C          the random sample of size K to Output File.
C Reference: Knuth, "Art of Compouter Programming", Vol 2,
C            Seminumerical Algorithms, Algorithm 3.4.2-S.
C Results: Exactly K records will be chosen;
C          Sample will be unbiased; that is the probability of
C          any record  being selected is K/N.
C HDK - 06/20/1997.
C
C---Arbitrary 80 character records.
      CHARACTER*80 RECORD
C---FIN is Input Filespec; FOUT is Output Filespec.
      CHARACTER*60 FIN,FOUT
      INTEGER I,K,M,N,SEED(3)
      REAL U(1)
      COMMON/RANDOM/SEED
C---Define N and K.
      WRITE(*,*) 'Please enter two integers, N and K:'
      READ(*,*) N, K
      IF(N.LT.K .OR. K.LT.0) THEN
        WRITE(*,*) 'Values of N and K must be such that:'
        WRITE(*,*) 'N >= K, K>0. Input Values are N=',N,'  K=',K
        STOP
      ENDIF
C---Open INPUT and OUTPUT files.
      FIN='KNUTHS.IN'
      FOUT='KNUTHS.OUT'
      OPEN(UNIT=50,STATUS='OLD',FILE=FIN)
      OPEN(UNIT=60,STATUS='UNKNOWN',FILE=FOUT)
      WRITE(*,*) 'INPUT FILE OF ',N,' RECORDS IS:', FIN
      WRITE(*,*) 'Please enter 3 random Integer Seeds;'
C---Initialize Random Number Generator, RAND, starting place.
      WRITE(*,*) '(Integers between 1 and 32760):'
      READ(*,*) (SEED(I),I=1,3)
C===============================================================
C---BEGIN Algorithm 3.4.2-S
C---Initialize
1     T=0
      M=0
C---Generate a U(0,1) number.
2     CALL RAND(U,1)
C---Test.
      IF ((N-T)*U(1) .GE. K-M) THEN
C---Skip next record.
        READ(50,'(A)')
        T=T+1
        GO TO 2
      ELSE
C---Select next Record and Copy it.
        READ(50,'(A80)') RECORD
        WRITE(60,'(A80)') RECORD
        M=M+1
        T=T+1
        IF (M .LT. K) GO TO 2
C---Terminate, K records were selected at random.
      ENDIF
      WRITE(*,*) 'OUTPUT FILE OF ',M,' RECORDS IS:', FOUT
      STOP
      END
      SUBROUTINE RAND(U,N)
C=======================================================================
C  Portable pseudo-random integer generator, especially for
C  microcomputers with a limitation of 16 bit integers. Translated from
C  original Pascal version(1) to Fortran 77 by H. D. Knoble, PSU.
C
C   The supporting paper is:
C   (1) B. Wichmann & D. Hill, "Building a Random-Number Generator",
C             BYTE, March, 1987, 12(3):127-128.
C
C   Also see the following related works:
C   (2) Wichmann, B.A. and Hill, I.D., "An Efficient and Portable",
C             Pseudo-random Number Generator", Applied Statistics,
C             Volume 31, 1982, pages 188-190.
C   (3) Haas, Alexander, "The Multiple Prime Random Number Generator",
C             ACM Transactions on Mathematical Software; December,
C             1987, 13(4):368-381.
C   (4) L'Ecuyer, Pierre, "Efficient and Portable Combined Random Number
C             Generators", Communications of the ACM; June, 1988,
C             31(6):742-749,774.
C
C     Purpose: Portable quasi-random number generator package with
C              good statistical properties that utilizes only 16-bit
C              integers (good for mainframe and microcomputers).
C     Some Properties: - Portable Fortran
C                      - uses only 16 bit integer constants
C                      - very long sequence; good statistical properties
C                      - random repeats (a rare property)
C                      - easily adaptable to other languages
C                      - well-documented in qualified literature.
C     H. D. Knoble, August 1987, HDK at PSUVM
C     Penn State University Center, for Academic Computing
C
C
C Use...
C      1) Define SEED(1), SEED(2), SEED(2) as three integers between
C         1 and 2**15-1.
C      2) CALL RAND(U,N)
C         To generate a sequence, U, of N Uniform(0,1) numbers.
C
C      Cycle length is ((30269-1)*(30307-1)*(30323-1))/4  or
C      6953607871644  > 6.95E+12.
C
C     To access the SEED vector in the calling program use statements:
C     INTEGER SEED(3)
C     COMMON/RANDOM/SEED
C
C  The common variable SEED is the array of three current seeds.
      INTEGER N,X,Y,Z
      REAL U(N),V
      COMMON/RANDOM/X,Y,Z
      IF(N.LE.0) RETURN
      DO 1 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/30308.0 + Z/30323.0
1       U(I)=V-INT(V)
      RETURN
      END