Tweet
'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