'RandINT.BAS -- A Swiss Army linear congruential pseudo-random number generator

' by E.W. Menzel Jr. March 2005.

' THIS IS FREEWARE

' Language: PowerBASIC PBCC

' Purpose: Provide a choice of random number generators -- all with known, open-source,

' well-tested algorithms, & some with much longer periods than PB's built-in RND. You can

' even switch your algorithms on the fly, should you care to do so.

' NOTES:

' 1. To "RANDOMIZE" simply set Rand.Seed = any positive integer between 1 and Rand.M -2.

' PB's RANDOMIZE and RND are not, however, used here.

' 2. For the theory behind these and other generators see:

' S.K. Park & K.W. Miller, "Random number generators: Good ones are hard to find",

' Communications of the ACM, Oct 1988, vol 31, 1192-1201.

' Donald Knuth, "The art of computer programming", 1998, vol 2, Chapter 3, 3rd edition.

' 3. Knuth lists many good (and a few "infamous") pairs of constants, together with their

' performances on tests of randomness. Here I pick four good pairs. If you need more than

' that, it will be easy to add them.

------------------

[This message has been edited by Emil Menzel (edited March 23, 2005).]

' by E.W. Menzel Jr. March 2005.

' THIS IS FREEWARE

' Language: PowerBASIC PBCC

' Purpose: Provide a choice of random number generators -- all with known, open-source,

' well-tested algorithms, & some with much longer periods than PB's built-in RND. You can

' even switch your algorithms on the fly, should you care to do so.

' NOTES:

' 1. To "RANDOMIZE" simply set Rand.Seed = any positive integer between 1 and Rand.M -2.

' PB's RANDOMIZE and RND are not, however, used here.

' 2. For the theory behind these and other generators see:

' S.K. Park & K.W. Miller, "Random number generators: Good ones are hard to find",

' Communications of the ACM, Oct 1988, vol 31, 1192-1201.

' Donald Knuth, "The art of computer programming", 1998, vol 2, Chapter 3, 3rd edition.

' 3. Knuth lists many good (and a few "infamous") pairs of constants, together with their

' performances on tests of randomness. Here I pick four good pairs. If you need more than

' that, it will be easy to add them.

Code:

#COMPILE EXE #DIM ALL DECLARE FUNCTION RandInt() AS QUAD DECLARE FUNCTION RNG() AS EXT DECLARE FUNCTION RandPeriod() AS QUAD DECLARE FUNCTION RandPick(whichOne AS LONG, seed AS QUAD) AS QUAD 'seed & constants for linear congruential random number generators TYPE RandType seed AS QUAD A AS QUAD 'constant a M AS QUAD 'constant m END TYPE GLOBAL Rand AS RandType 'Function RandInt (plus the Type definition above) might be all you 'really need; the rest is window-dressing. It returns a 10-digit integer. 'Note that if the generator has not been defined, default values are used. FUNCTION RandINT AS QUAD DIM Rand AS GLOBAL RandType IF Rand.M < 1&& THEN Rand.A=16807&&: Rand.M=(2^31)-1&& IF Rand.seed < 1&& THEN Rand.seed=TIMER Rand.seed = (Rand.seed * Rand.A) MOD Rand.M FUNCTION = Rand.seed END FUNCTION 'RNG is an equivalent of PB's RND, 'in yielding uniformly distributed numbers ranging from 0 to 1, 'with theoretical exact mean=0.5 and var=(1/12) or 0.08333... FUNCTION RNG AS EXT FUNCTION = RandINT / Rand.M END FUNCTION 'which generator & starting seed number do you want? FUNCTION RandPick(whichOne AS LONG, seed AS QUAD) AS QUAD DIM A AS STRING, N AS LONG N =4 'however many generators are listed below IF seed <1 OR whichOne<1 OR whichOne>N THEN CLS LOCATE 18,1 END IF IF WhichOne >N OR WhichOne <1 THEN 'Rand. is a GLOBAL TYPE variable. Function RandInt DIM's it and will screen it. PRINT "Which random number generator do you want?" PRINT "(1)=Lavaux&Jannsens (2)=Cray X-MP (3) Marsaglia (4)Park&Miller " PRINT "[default = ";N; LINE INPUT "] ";A$ WhichOne=VAL(A$) END IF SELECT CASE WhichOne CASE 1: Rand.A=31167285&&: Rand.M=2^48 'see Knuth, 1998, p.106, line 23 CASE 2: Rand.A=44485709377909: Rand.M=2^46 'see Knuth, 1998, p.106, line 22 CASE 3: Rand.A=69069: Rand.M=2^32 'see Knuth, 1998, p.106, line 22; G.Marsaglia CASE ELSE: Rand.A=16807&&: Rand.M=(2^31)-1 'see Knuth, 1998, p. 106, line 19; Park-Miller END SELECT IF seed THEN Rand.seed=seed ELSE LINE INPUT "Enter a seed number (0 or neg. = TIMER) ";A$ Rand.seed = VAL(a$) END IF IF Rand.seed < 1 THEN Rand.seed=INT(TIMER) FUNCTION=Rand.seed END FUNCTION 'Note: the period for PB's RND is 2^32 or over 4 billion, and 'the period for the generator using Rand.A=31167285&& & Rand.M=2^48 'is 2^48. Don't expect Function RandPeriod to reach the latter number 'any time soon. (I quit after a day or so, and count of a mere 350 billion.) FUNCTION RandPeriod AS QUAD DIM X(1 TO 5) AS QUAD, N AS QUAD, Hit AS LONG, UB AS LONG DIM Y AS LONG, I AS LONG, J AS LONG UB=UBOUND(X) Y=24 LOCATE Y-1,10: PRINT "Testing for RandINT's period or length of cycle" FOR N=1 TO UB X(N)=RandInt NEXT N DO FOR I=1 TO 1000000 INCR N IF RandInt=X(1) THEN Hit=1 FOR J=2 TO UB IF RandInt=X(J) THEN INCR Hit NEXT J IF Hit >= UB THEN EXIT DO END IF NEXT I LOCATE Y,10: PRINT N; LOOP UNTIL INSTAT IF INKEY$<>"" THEN N = -N LOCATE Y,10: PRINT N PRINT "press any key to continue" WAITKEY$ FUNCTION = N END FUNCTION 'demo 'Note, from RandINT above, that Rand. is a GLOBAL TYPE variable & that Function RandInt DIM's it. 'You can change Rand.seed, Rand.a, & Rand.m values as you wish -- but don't say I didn't warn you. FUNCTION PBMAIN () AS LONG DIM I AS LONG, A AS STRING, Period AS QUAD, RandConstM AS QUAD DIM freq(9) AS LONG, Y AS LONG DIM X AS EXT, sumX AS EXT, sumX2 AS EXT, N AS EXT, mean AS EXT, var AS EXT Y=RandPick(4,666) 'pick generator# & seed.. If you don't, RandINT will use its defaults Period=RandPeriod 'Start FUNCTION to determine how many randoms we can draw before the same 'sequence repeats. (Note: This could take many days for some generators) Rand.seed = 1 're-seed & start over DO PRINT FOR I=1 TO 10000 X= RNG 'same as RandINT / Rand.M sumX = sumX+X sumX2=sumX2+X*X: N=N+1 A$=FORMAT$(X,".000000000000000000 ") IF I>9900 THEN PRINT A$; A$=MID$(A$,2,1) ' #2 is first digit within the fractional number Y=VAL(A$) INCR freq(Y) NEXT PRINT PRINT "frequencies of the digits 0-9" FOR I=0 TO UBOUND (freq) PRINT freq(I); NEXT PRINT PRINT "current seed num ";Rand.seed 'If initial seed was 1 and seed #10,001 is not 1043618065, exactly, then ' the Park & Miller minimal standard RNG was not implemented IF N=10000 AND Rand.seed=1043618065 THEN PRINT "Alright already... You used seed=1, a=16807, m=(2^31)-1. " END IF mean=sumX/N: var=(sumX2-(sumX^2/N)) /(N-1) A$="mean "+ STR$(mean##,18)+" variance "+STR$(var##,18)+ " N "+STR$(N) PRINT A$ PRINT "Enter q to quit, any number 1-"; Rand.M-1; LINE INPUT " for new seed ";A$ IF VAL(A$) THEN Rand.seed = VAL(a$) LOOP UNTIL UCASE$(A$)="Q" END FUNCTION

[This message has been edited by Emil Menzel (edited March 23, 2005).]

## Comment