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