CFrom HDK@psuvm.psu.edu Thu Aug 25 09:11:00 MDT 1994
C
C A PORTABLE QUASI-RANDOM NUMBER GENERATOR
C
CThere have been many requests for a portable Fortran program that
Cgenerates quasi-random (quasi- meaning "apparent" "near" versus pseudo-
Cmeaning "false") sequences of numbers with good statistical properties.
CBy "portable" here is meant limited to use of 16-bit integers as well as
Cportable Fortran 77. The advantages of this generator over others is its
Cstatistical properties. In addition it has a very long cycle length
C(6.9E+12), and produces quasi-randomly REPEATED numbers in its output
Csequences. The algorithm has been covered rather extensively in refereed
Cliterature which is cited in the comments of the BLOCK DATA subprogram.
C
CThe subroutine RAND and companion BLOCK DATA subprogram comprise the
Calgorithm itself. With them are also a sample driver (main) program, and
Ctwo subroutines for random sampling - one with replacement (IRAND) and
Cthe other (more difficult case, RPERM) without replacement. The Fortran
Ccode follows:
C=====Sample driver for portable pseudo-random number generator, RAND.
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 PSU.EDU
C Penn State University Center, for Academic Computing
C
C Author's statement of intent: Program units following this notice
C may be freely copied by anyone under the condition that they be used
C only for peaceful purposes that do not threaten human dignity and
C life or the environment.
C
C---This driver program is used to illustrate how to use the subprograms
C RAND, RPERM, and IRAND by interactively generating quasi-random
C permutations of the first N integers. To run it on the Penn State
C VM/CMS system issue:
C 1) For WATFOR77 issue: WATFOR77 BYTERAND
C or 2) For VS Fortran issue: FORTVS2 BYTERAND
C LOAD BYTERAND (START
C or 3) For IBM PC WATFOR77, download and copy this program to
C d:\path, and issue: WATFOR77 d:\path\BYTERAND.FOR
C
C To quit this sample driver program simply press Enter
C in response to any prompt.
INTEGER X(999),N
C---Maximum value of N is dimension of X.
INTEGER SEED(3)
COMMON/RANDOM/SEED
C
WRITE(6,*) '** Random Permutation of the first N Integers **'
WRITE(6,*) 'Please enter the value of N (or press Enter to quit):'
READ(5,*,END=5) N
IF (N.EQ.0.OR.N.GT.999) STOP
DO 9999 NSETS=1,9999
C---Generate a random permutation of the 1ST N integers.
*C---Generate a "random" seed using Time-Of-Day clock (VS Fortran only)
* DOUBLE PRECISION TOD
* INTEGER ISEED
*C---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(6,*) 'Please enter 3 integer seeds, each >0 & <30000:'
READ(5,*,END=5) SEED
IF(SEED(1).EQ.0) STOP
C---Generate a random permutation, X
CALL RPERM(X,N)
C---Output the random permutation.
WRITE(6,2) N
2 FORMAT(' Random permutation of the 1st ',I4,' integers:')
WRITE(6,4) (X(I),I=1,N)
4 FORMAT(20I4)
9999 CONTINUE
5 STOP
END
SUBROUTINE RPERM(P,N)
C===Generate a random permutation, P, of the first N integers.
C (equivalent to sampling WITHOUT REPLACEMENT).
C Adaptation of Knuth Volume 2, Algorithm 3.4.2P.
INTEGER N,P(N), K,J,I,IPJ,ITEMP,M
REAL U(100)
DO 1 I=1,N
1 P(I)=I
C---Generate up to 100 U(0,1) numbers at a time.
DO 3 I=1,N,100
M=MIN(N-I+1,100)
CALL RAND(U,M)
DO 2 J=1,M
IPJ=I+J-1
K=INT(U(J)*(N-IPJ+1))+IPJ
ITEMP=P(IPJ)
P(IPJ)=P(K)
2 P(K)=ITEMP
3 CONTINUE
RETURN
END
SUBROUTINE IRAND(S,N,LOW,HI)
C===Generate a random integer sequence: S(1),S(2), ... ,S(N)
C such that each element is in the closed interval and
C sampled WITH REPLACEMENT. HDK, JUNE 1971.
INTEGER N,S(N),LOW,HI,IX,I
REAL U(1)
DOUBLE PRECISION X
DO 1 I=1,N
CALL RAND(U,1)
C---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
1 CONTINUE
RETURN
END
BLOCK DATA
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 Use...
C CALL RAND(U,N)
C To generate a sequence, U, of N Uniform(0,1) numbers.
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 SEED(3)
COMMON/RANDOM/SEED
DATA SEED(1),SEED(2),SEED(3)/1,10000,3000/
END
C=======================================================================
SUBROUTINE RAND(U,N)
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