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