### Author Topic: FAQ suggestion: RND  (Read 263 times)

#### mcalkins

• Hero Member
•     • Posts: 1269
• •  ##### FAQ suggestion: RND
« on: February 21, 2013, 12:27:29 PM »
Quote
FAQ: How random is RND?
FAQ: How do RND and RANDOMIZE work?
FAQ: What if I need a better algorithm?

Prerequisites.
This FAQ article assumes that you are familiar with hexadecimal and binary numbers, and with bitwise operations.

http://en.wikipedia.org/wiki/Binary_number
http://en.wikipedia.org/wiki/Bitwise_operation
http://www.qb64.net/forum/index.php?topic=5990.0

http://www.qb64.net/wiki/index.php?title=%26H

What is RND?

QBASIC and QB64 provide a pseudorandom number generator through their RANDOMIZE and RND keywords.

How random is RND? It is random enough for casual purposes where there is nothing really at stake. For example, it is sufficiently random for most games that don't involve serious competition. It is not suitable to be at the core of any cryptographic system.

At a high level, RND returns its values from a repeating sequence of 2 ^ 24 (decimal 16,777,216) numbers, with each number occurring only once. RANDOMIZE changes your position in that sequence.

RND is actually a 24 bit Linear Congruential Generator (LCG). Please read the Wikipedia article: http://en.wikipedia.org/wiki/Linear_congruential_generator

At RND's core is this algorithm:

Code: [Select]
`state = &HFD43FD * state + &HC39EC3 AND &HFFFFFF`
To use the Wikipedia terminology:
The modulus (m) is 2^24, hexadecimal: 1000000, decimal: 16777216
The multiplier (a) is hexadecimal: fd43fd, decimal: 16598013
The increment (c) is hexadecimal: c39ec3, decimal: 12820163
The initial state (X0) is hexadecimal: 50000, decimal: 327680

These numbers satisfy the requirements stated in the Wikipedia article for having a full period. So, the RND function will have a period of 2^24 for any seed.

You can think of the RND function as containing a list of 2^24 unique integers. Each integer in the range of 0 to &Hffffff occurs exactly once. Reseeding just changes your position in the sequence, but the sequence remains constant. RND, of course, does not contain such a list, but the algorithm behaves the same way.

RND's return value.

RND always returns as a floating point number the internal state divided by &h1000000 (2^24). This is why RND always returns a fractional number >= 0, but < 1.

The traditional usage of RND is:

Code: [Select]
`value = INT(RND * numberOfPossibilities) + minimumValue`
So, for example, to duplicate a die:

Code: [Select]
`die = INT(RND * 6) + 1`
This kind of usage gives more significance to the higher order bits. This is in harmony with Wikipedia's statement that the higher order bits are of better quality than the lower order bits.

RND with a parameter of 0.

When RND is called with a parameter of 0, it does not perform the LCG algorithm operation. Instead, it simply divides the state by 2^24, and returns the result.

Code: [Select]
`RND = state / &h1000000`
You can use this program to reveal QBASIC's initial state:

Code: [Select]
`CLSPRINT "&h"; HEX\$(RND(0) * &H1000000)END`
It displays &h50000.

RND with no parameter, or with a positive parameter.

The normal usage of RND is with no parameter. In this case, RND performs the LCG algorithm on the state, then returns the state divided by 2^24.

Code: [Select]
`state = &HFD43FD * state + &HC39EC3 AND &HFFFFFFRND = state / &h1000000`
Try it for yourself, starting from the initial state of &h50000:

&hfd43fd * &h50000 is &h4f253f10000
+ &hc39ec3 is &h4f254b49ec3
keep only the low 24 bits: &hb49ec3, decimal 11837123
/ &h1000000 is 0.705547511577606201171875

Let us iterate again:

&hfd43fd * &hb49ec3 is &hb2b0dec4efb7
+ &hc39ec3 is &hb2b0df888e7a
keep only the low 24 bits: &h888e7a, decimal 8949370
/ &h1000000 is 0.53342401981353759765625

This program shows that QBASIC agrees:

Code: [Select]
`CLSPRINT RND, "&h"; HEX\$(RND(0) * &H1000000)PRINT RND, "&h"; HEX\$(RND(0) * &H1000000)END`

RND with a negative parameter.

When given a negative parameter, RND uses it to reseed its state. Then it performs the LCG operation. Then it returns the state divided by 2^24.

For this reseeding, RND relies on the bit pattern of the parameter stored as a SINGLE. SINGLEs are single precision IEEE 754-1985 floating point numbers stored in little-endian byte order.

The 32 bits of the SINGLE are treated as an _UNSIGNED LONG. The new state will be the high 8 bits rightshifted by 24 bits, then added to the bit pattern, and masked to 24 bits. Then, RND proceeds as it normally would.

Since MKS\$ returns the bit pattern of a SINGLE, and CVL turns a bit pattern into a LONG, this could be represented as:

Code: [Select]
`t = CVL(MKS\$(parameter))state = (t + (t \ &H1000000)) AND &HFFFFFFstate = &HFD43FD * state + &HC39EC3 AND &HFFFFFFRND = state / &H1000000`
Try it yourself, with a parameter of -8.75:

See: http://en.wikipedia.org/wiki/IEEE_754-1985

-8.75 is binary: -1000.11   (that is, -(8 + 1/2 + 1/4))
This is: binary -1.00011 * decimal 2^3
So, to express it as a SINGLE, the integer 1 portion will be implied, so we drop it, and we add 127 to the exponent to form the biased exponent. So:

the 1 bit sign: 1
the 8 bit biased exponent: decimal 130, binary 10000010
the 23 bit fraction: binary 00011

So, -8.75 is:
1 10000010 00011000000000000000000

Numbers are stored in little-endian byte order. See: http://en.wikipedia.org/wiki/Endianness
So, the -8.75 as a SINGLE is actually stored as the series of 4 bytes:
&h00 &h00 &h0c &hc1

You can test this in QB64:

Code: [Select]
`_CONTROLCHR OFFPRINT "'"; MKS\$(-8.75); "'"END`
It will display: '  ♀┴' (blanks are nulls.)

This will be reinterpreted as an _UNSIGNED LONG. Because of the little-endian byte order, it will be:
&hc10c0000.

You can test this in QBASIC with:

Code: [Select]
`CLSPRINT "&h"; HEX\$(CVL(MKS\$(-8.75)))END`
As a reminder, this is what we are doing:

Code: [Select]
`t = CVL(MKS\$(parameter))state = (t + (t \ &H1000000)) AND &HFFFFFFstate = &HFD43FD * state + &HC39EC3 AND &HFFFFFFRND = state / &H1000000`
For a parameter of -8.75, t is: &hc10c0000
rightshifted 24 bits: &hc1
+ &hc10c0000 is &hc10c00c1
keep only the low 24 bits: &hc00c1
* &hfd43fd is: &hbdfeecc41bd
+ &hc39ec3 is: &hbdfef8fe080
keep only the low 24 bits: &h8fe080, decimal 9429120
/ &h1000000 is: 0.56201934814453125

QBASIC agrees:

Code: [Select]
`CLSPRINT RND(-8.75), "&h"; HEX\$(RND(0) * &H1000000)END`

RANDOMIZE USING.

QB64 provides RANDOMIZE USING to reseed the RND state.

RANDOMIZE USING relies on the bit pattern of its parameter stored as a DOUBLE, a double precision IEEE 754-1985 floating point number stored in little-endian byte order.

The high 32 bits of the DOUBLE are treated as an _UNSIGNED LONG. That _UNSIGNED LONG will be rightshifted 16 bits, then XORed with the low 16 bits of the original _UNSIGNED LONG, then leftshifted 8 bits, to form the new state.

It can be expressed as:

Code: [Select]
`t = CVL(RIGHT\$(MKD\$(parameter), 4))state = ((&HFFFF& AND t) XOR t \ &H10000) * &H100`
Try it for yourself with a parameter of 3.5:

3.5 is: binary 11.1   (that is, 2 + 1 + 1/2)
It is: binary 1.11 * decimal 2^1
Remember that the 1 is implied, so we drop it and keep the fraction. Also, we add 1023 to bias the exponent.

the 1 bit sign: 0
the 11 bit biased exponent: decimal 1024, binary 10000000000
the 52 bit fraction: binary 1100000000000000000000000000000000000000000000000000

So, 3.5# is:
0 10000000000 1100000000000000000000000000000000000000000000000000
&h400c000000000000

Because it is little-endian, it is stored as bytes:
&h00 &h00 &h00 &h00 &h00 &h00 &h0c &h40

The high 32 bits will be reinterpreted as an _UNSIGNED LONG. Because of little-endian byte order, it will be:
&h400c0000

(You could also think of it as reinterpreting the DOUBLE as an _UNSIGNED _INTEGER64, and rightshifting it 32 bits.)

The QB64 program:

Code: [Select]
`_CONTROLCHR OFFPRINT "'"; MKD\$(3.5); "'"END`
Displays: '      ♀@' (blanks are nulls.)

The QBASIC program:

Code: [Select]
`CLSPRINT "&h"; HEX\$(CVL(LEFT\$(MKD\$(3.5), 4)))PRINT "&h"; HEX\$(CVL(RIGHT\$(MKD\$(3.5), 4)))END`
confirms our results so far.

As a reminder, we are doing:

Code: [Select]
`t = CVL(RIGHT\$(MKD\$(parameter), 4))state = ((&HFFFF& AND t) XOR t \ &H10000) * &H100`
So, for a parameter of 3.5, t is &h400c0000
&h400c0000 rightshifted 16 bits is: &h400c
The low 16 bits of &h400c0000 are: 0
0 XOR &h400c is: &h400c
leftshifted 8 bits is: &h400c00

We can test in QBASIC. (QBASIC doesn't have RANDOMIZE USING, but if we start with RANDOMIZE, the result is the same. You'll see why in the next section.):

Code: [Select]
`CLSRANDOMIZE 3.5PRINT "&h"; HEX\$(RND(0) * &H1000000)END`

RANDOMIZE.

RANDOMIZE, without USING, is almost the same as RANDOMIZE USING. The only difference is that RANDOMIZE ORs the result with the low 8 bits of the old state.

Code: [Select]
`t = CVL(RIGHT\$(MKD\$(parameter), 4))state = ((&HFFFF& AND t) XOR t \ &H10000) * &H100 OR &HFF AND state`
Therefore, as long as the low 8 bits of the existing state are 0, RANDOMIZE and RANDOMIZE USING produce the same result. Recall that the initial state is &h50000. Therefore, RANDOMIZE and RANDOMIZE USING have identical results if the state is at its initial value. This is in harmony with the purpose of RANDOMIZE USING, to conveniently reseed as if using RANDOMIZE at the start of the program.

For our experiment, from the initial state of &h50000, we will iterate the LCG once. This produces a state of &hb49ec3, as shown earlier. From here, we will RANDOMIZE with a parameter of 5.25:

5.25 is: binary 101.01   (that is, 4 + 1 + 1/4)
It is: binary 1.0101 * decimal 2^2

Code: [Select]
`sign          fraction0 10000000001 0101000000000000000000000000000000000000000000000000  biased exponent`
&h4015000000000000

MKD\$(5.25) would be: "      §@" (blanks are nulls.)

As a reminder, we are doing:

Code: [Select]
`t = CVL(RIGHT\$(MKD\$(parameter), 4))state = ((&HFFFF& AND t) XOR t \ &H10000) * &H100 OR &HFF AND state`
the existing state is: &hb49ec3
t is: &h40150000
&h40150000 rightshifted 16 bits is: &h4015
the low 16 bits of &h40150000 are: 0
0 XOR &h4015 is: &h4015
leftshifted 8 bits is: &h401500
the low 8 bits of &hb49ec3 are: &hc3
&h401500 OR &hc3 is: &h4015c3

This is confirmed in QBASIC:

Code: [Select]
`CLSPRINT "&h"; HEX\$(RND * &H1000000)RANDOMIZE 5.25PRINT "&h"; HEX\$(RND(0) * &H1000000)END`
So far, the low 16 bits of t have been 0. Let's use a more complicated parameter for RANDOMIZE:

4*ATN(1)
MKD\$(4 * ATN(1)) is "↑-DT√!○@"
&h400921fb54442d18

From the state of &h4015c3 in the last example, let us iterate the LCG once, producing a state of &h82297a, and then try the RANDOMIZE 4 * ATN(1):

Code: [Select]
`t = CVL(RIGHT\$(MKD\$(parameter), 4))state = ((&HFFFF& AND t) XOR t \ &H10000) * &H100 OR &HFF AND state`
the existing state is: &h82297a
t is: &h400921fb
&h400921fb rightshifted 16 bits is: &h4009
the low 16 bits of &h400921fb are: &h21fb
&h21fb XOR &h4009 is: &h61f2
leftshifted 8 bits is: &h61f200
the low 8 bits of &h82297a are: &h7a
&h61f200 OR &h7a is: &h61f27a

QBASIC confirms:

Code: [Select]
`CLSPRINT "&h"; HEX\$(RND * &H1000000)RANDOMIZE 5.25PRINT "&h"; HEX\$(RND * &H1000000)RANDOMIZE 4 * ATN(1)PRINT "&h"; HEX\$(RND(0) * &H1000000)END`
The traditional usage of RANDOMIZE is to have:

Code: [Select]
`RANDOMIZE TIMER`
once before using RND. Additional occasional uses of RANDOMIZE TIMER might help improve the results slightly, by relocating yourself in the sequence of RND states. However, I doubt that the benefit is enough to justify it, in most cases.

Sample implementation.

Code: [Select]
`' derived from: QB64's libqbx.cpp (versions 0.942 and 0.954) and http://www.network54.com/Forum/178387/message/1046747461/' I do not guarantee correctness.DIM SHARED qState AS _UNSIGNED LONG: qState = &H50000test: ENDFUNCTION qRnd!qState = &HFD43FD * qState + &HC39EC3 AND &HFFFFFFqRnd = qState / &H1000000END FUNCTIONSUB qRandomizeUsing (s AS DOUBLE)qState = CVL(RIGHT\$(MKD\$(s), 4))qState = ((&HFFFF& AND qState) XOR qState \ &H10000) * &H100END SUBSUB qRandomize (s AS DOUBLE)DIM l AS _UNSIGNED LONGl = CVL(RIGHT\$(MKD\$(s), 4))qState = ((&HFFFF& AND l) XOR l \ &H10000) * &H100 OR &HFF AND qStateEND SUBFUNCTION qRndPrevious!qRndPrevious = qState / &H1000000END FUNCTIONFUNCTION qRndSeed! (s AS SINGLE)DIM l AS _UNSIGNED LONGl = CVL(MKS\$(s))qState = (l + (l \ &H1000000)) AND &HFFFFFF ' I apply the bitmask after the addition. libqbx applies it before the addition. I don't know if that could cause a carry into the 25th bit in libqbx.qState = &HFD43FD * qState + &HC39EC3 AND &HFFFFFFqRndSeed = qState / &H1000000END FUNCTIONSUB testDIM i AS LONGPRINT "initial default seed": GOSUB pRANDOMIZE 5.25: qRandomize 5.25PRINT "randomize": GOSUB pRANDOMIZE USING -3.5: qRandomizeUsing -3.5PRINT "randomize using": GOSUB pPRINT "rnd(0)"PRINT RND(0), qRndPreviousPRINT RND(0), qRndPreviousPRINT "rnd(negative)"PRINT RND(-8.75), qRndSeed(-8.75)GOSUB pEXIT SUBp:FOR i = 0 TO 2: PRINT RND, qRnd: NEXTRETURNEND SUB`

What if I need a better algorithm?

http://en.wikipedia.org/wiki/Pseudorandom_number_generator
http://en.wikipedia.org/wiki/List_of_pseudorandom_number_generators
http://en.wikipedia.org/wiki/Cryptographically_secure_pseudorandom_number_generator

I'm not very familiar with the other PRNGs, so I can't offer very much specific advice.

The algorithms listed on those pages vary significantly in simplicity/complexity, execution speed, memory usage, quality, and cryptographic security. Note that some of them are not any more secure than LCGs.

You might find some of them already implemented on the forums.

I have sample code to use the Windows CryptGenRandom function here:
http://www.qb64.net/forum/index.php?topic=6891.msg70629#msg70629
However, I offer it with no guarantees. If you use it, examine it yourself, and use it at your own risk.

http://en.wikipedia.org/wiki/Linear_congruential_generator
http://www.network54.com/Forum/178565/message/1046721807/
http://www.network54.com/Forum/178387/message/1046747461/

http://www.qb64.net/wiki/index.php?title=RND
http://www.qb64.net/wiki/index.php?title=RANDOMIZE
http://www.qb64.net/wiki/index.php?title=INT

http://www.qb64.net/wiki/index.php?title=CVL
http://www.qb64.net/wiki/index.php?title=MKS\$
http://www.qb64.net/wiki/index.php?title=MKD\$
http://en.wikipedia.org/wiki/IEEE_754-1985
http://en.wikipedia.org/wiki/Endianness

http://en.wikipedia.org/wiki/Pseudorandom_number_generator
http://en.wikipedia.org/wiki/List_of_pseudorandom_number_generators
http://en.wikipedia.org/wiki/Cryptographically_secure_pseudorandom_number_generator

Clippy: Please wait at least a few days before putting this in the Wiki. It might need to be revised, and other people may have suggestions.

If people have implementations of other PRNGs, it might be appropriate to link to them in the "What if I need a better algorithm?" section.

I'm too tired to double check the technical stuff. I'll just give it a read through to make sure it sounds right...

Regards,
Michael
The QBASIC Forum Community: http://www.network54.com/index/10167 Includes off-topic subforums.
QB64 Off-topic subforum: http://qb64offtopic.freeforums.org/

#### SMcNeill

• Hero Member
•     • • Posts: 2414
•  ##### Re: FAQ suggestion: RND
« Reply #1 on: February 21, 2013, 01:43:43 PM »
Quote
If people have implementations of other PRNGs, it might be appropriate to link to them in the "What if I need a better algorithm?" section.

Here's a really simple little random routine which I've used in the past.  It doesn't tend to have the same repeatable pattern RND does, as we shuffle it each time when we get to the end of the list again.  Make your limit as large as you want, as larger numbers tend to give a smoother spread of distributions (due to rounding issues).  With 10,000 we generate a nice spread of the numbers 0 -9 with 1000 of 8 of them, 1001 of 1 of them, and 999 of the last.  Not quite a perfect distribution of the values, but mighty dang close!

The only issue with larger Limits? More memory required to store the array, and more time to shuffle it once we run though the sequence once.   Of course, there'd be fewer calls to resequence the array, but when they did occur they'd take longer to do.

Code: [Select]
`CONST Limit = 10000#DIM SHARED PRNG(Limit) AS DOUBLEFOR i = 1 TO Limit    x = INT(Rndm * 10)    spread(x) = spread(x) + 1    IF i < 31 THEN PRINT x,NEXTPRINT: PRINTFOR i = 0 TO 9    PRINT spread(i),NEXTFUNCTION RndmSTATIC init, countIF NOT init THEN    init = -1    RANDOMIZE TIMER    FOR i = 0 TO Limit        PRNG(i) = i / (Limit + 1)    NEXT    FOR i = 0 TO Limit        x = INT(RND * Limit)        SWAP PRNG(i), PRNG(x)    NEXTEND IFcount = count + 1IF count > Limit THEN    RANDOMIZE TIMER    count = 0    FOR i = 0 TO Limit        x = INT(RND * Limit)        SWAP PRNG(i), PRNG(x)    NEXTEND IFRndm = PRNG(count)END FUNCTION`
http://bit.ly/TextImage -- Library of QB64 code to manipulate text and images, as a BM library.
http://bit.ly/Color32 -- A set of color CONST for use in 32 bit mode, as a BI library.

http://bit.ly/DataToDrive - A set of routines to quickly and easily get data to and from the disk.  BI and BM files

#### Clippy

• Hero Member
•     • • Posts: 16439
• I LOVE π = 4 * ATN(1)    Use the QB64 WIKI >>>
• •  ##### Re: FAQ suggestion: RND
« Reply #2 on: February 21, 2013, 02:25:07 PM »
There is no such thing as more random or less random in CHAOS. Any pattern could be replicated even in pure randomness, whatever that is. For its purposes, RND does the job 99% of the time!

We don't care about the 1%, they can take care of themselves .... QB64 WIKI: Main Page

#### Mrwhy

• Hero Member
•     • Posts: 2893
• My Dad called me Mr Why when I was 5.
•  ##### Re: FAQ suggestion: RND
« Reply #3 on: February 21, 2013, 04:38:56 PM »
How random is RND. Fourteen tests it shouuld  (but does not) pass:-

Here are the results obtained by my friend Delmer in USA
Code: [Select]
`                        Is Your RND Really Random?                       ÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ                       Testing the BASIC RND Function                        by   Delmer D. Hinrichs                           2116 S.E. 377th Avenue                           Washougal, Wash.  98671-9725         If you want to be absolutely correct, the answer to the title questionmust be `No'. The Basic RND function generates each new `random' numberfrom the previous number, starting from an initial `seed' number, by asimple calculation.   Thus these apparently random numbers are morecorrectly called `pseudo-random' numbers.   For simplicity, from now on wewill just call them `random' numbers.   Which still leaves us with thequestion:   Are these apparently random numbers close enough to being trulyrandom to be useful?      But useful for what?   If you use the RND function merely to get adifferent starting point for a Lunar Lander game, it doesn't make muchdifference how random it is.   If you use RND to simulate dealing cards fora poker game, you may lose more games than you should if the RND is notrandom, but it still isn't really serious.   However, if you use RND in aMonte Carlo simulation to determine the pricing policy for your business,or for some other serious use, you must make sure that your RND isadequately random before you can trust the results.   For these criticaluses, or even to check on the honesty of your poker dealer, there are teststhat can be run on your RND function.       How can we check how random our RND function is?   Strangely enough,as the value of each individual random number becomes more uncertain (moretruly random), the statistics of a large number of these random numbersbecome more certain.   Therefore we can let our RND function produce alarge number of random numbers and check their statistics.      Which statistics, and how do we check them?   There are many differentstatistics that could be used, as random numbers have many different randomcharacteristics.   The `RNG' test program has fourteen different routinesthat test for sixteen different random number characteristics.   A trulyrandom number generator (RNG) would pass all of these tests.   Any one RNDfunction may pass some tests easily, but fail others.   The uses to whichyou put your RND determine which tests are most important to you.      How can we tell if an RNG fails to pass a test?   We cannot reallytrust our unaided judgement on this, as a series of numbers that appearrandom may not be random, and vice versa.   For example, if we throw asingle die, each of the six faces should come up about one-sixth of thetime.   If we throw the die 600 times, would we expect each face to come upexactly 100 times?   I think that we can agree that this is highlyunlikely.   But just how much should the number of times that each facecomes up differ from 100?   One of the most useful statistical tests forthis type of problem is the chi square test.      To apply the chi square test to our die-throwing experiment, wecompute:                         (n(i) - 100)ý                  X = ä ÄÄÄÄÄÄÄÄÄÄÄ                           100 where:     X = the chi square statistic  (Greek letter chi)           ä = summation symbol, indicating that all six computed values               are to be added together  (Greek letter sigma)          n(i) = number of times that face `i' comes up (i = integers 1-6)         100 = `expected' number of times that each face should come up      Now that we have the chi square, how do we use it?   To interpret thechi square, we need a chi square distribution table such as Table 1.There we find another new term, `degrees of freedom.' This is simply thenumber of independent choices possible, or one less than the number ofpossible categories.   For throwing a die, there are five degrees offreedom (if the face is not numbers 1 through 5, then it must be six, sothere are only 5 independent choices).      In Table 1, the 50% probability column is the `Ideal' chi squarelevel, with the results neither more uniform (probability > 50%, too smalla chi square) nor more scattered (probability < 50%, too large a chisquare) than would be obtained from truly random events.   If each of thesix faces of the die came up exactly 100 times, the chi square would bezero, or much more uniform than would be expected from random events.   Ifwe had three faces coming up 110 times each, and three faces coming up 90times each, the chi square would be 6.000, or just slightly more scatteredthan the ideal chi square value of 4.351.   Actually, any chi square valuefrom 1.15 to 11.1 (with probabilities of 95% to 5%) would be quitereasonable for the die-throwing experiment.   If the probability is greaterthan 99% or less than 1%, it is normally grounds for the rejection of theresult.   Of course, if you run hundreds of trials, you can expect aboutone of each 100 chi squares to have a probability greater than 99%, and thesame for a probability less than 1%.      Random events that occur as a continuous distribution instead of as asmall number of discrete categories are not suitable for the use of the chisquare statistic.   In some of these cases, the `Z score' statistic may beused to check how well the observed events fit the theoretical probabilities.The `Ideal' Z score is zero, with a 5% probability of its exceeding ñ1.96,and a 1% probability of its exceeding ñ2.58.      Some other types of random events are more difficult to treat probab-ilistically, but we may still gain some added information on randomness byjust looking at the results.      Now how do we apply these methods to testing a random number generator?Let us go through the 14 routines of the `RNG' program one by one.   I haveincluded the degrees of freedom (DF) for the chi square tests.   Note thatall test routines pause after display of the heading until you press a key,and may be stopped at any point during the run by pressing a key.   Allinterpreted Basic tests were run for 1,000 batches, and all compiled Basictests were run for 10,000 batches. 1.  Distribution Test   (DF = 9)      The distribution test (sometimes called the frequency test) checks ifthe RNG gives a uniform distribution of numbers.   That is, if all numberswithin its range are equally likely to occur.   This is one of the simplesttests, and most RNGs pass it.      Random integers (1 through 10) should each occur about 10% of thetime.   The program display of Table 5 shows the batch number, the countin each of the ten bins for 1000 random integers, the chi square of thatbatch, then the cumulative chi square for all batches taken together.   Theheading shows the Ideal (expected) value of each number.      The batch chi square measures the short-term randomness while thecumulative chi square measures the long-term randomness.   The RNG may passone chi square test but not the other.   For example, if the RNG produced aslight systematic excess of high integers, it could still pass the chisquare test for each batch, but the cumulative chi square would graduallybecome larger and larger, showing a greater and greater variation from truerandomness. 2.  Mean and Standard Deviation      Random numbers (0 to 1) should average about 0.5000, with a standarddeviation of about 0.2887.   For each batch of 1000 random numbers theprogram displays the batch number, the batch mean, the batch standarddeviation, the cumulative mean of all batches, and the cumulative standarddeviation of all batches. 3.  Correlation Test      Random numbers (0 to 1) should have a serial correlation coefficientof approximately zero.   That is, each new number should not depend uponthe previous number in any regular way.   A positive correlation shows atendency to run, while a negative correlation shows a tendency toalternate.      For each batch of 1000 sequential pairs of random numbers, the programshows the batch number, the correlation coefficient, the batch Z score, andthe cumulative Z score for all batches taken together. 4.  Serial Test   (DF = 8)      Random integers (0, 1, or 2) should occur in pairs in equal numbers.That is, the pair 0,0 should occur about as often as 0,1 or 0,2 or 1,0 etc.      For each batch of 1000 pairs, the program shows the batch number, thecount of each of the nine possible pairs, the batch chi square, and thecumulative chi square of all batches taken together. 5.  Run Test   (DF = 4)      Random numbers (0 to 1) should have `runs up,' where each succeedingnumber is larger than the last, in a predictable manner.   For this test,the count of the number of random numbers needed to give 1000 `runs up'will also vary a little.      For each batch of 1000 `runs up,' the program shows the batch number,the count for the batch, the number of `runs up' of length zero (the nextrandom number was smaller), one, two, three, and four or more, then thebatch chi square and the cumulative chi square for all batches. 6.  Gap Test   (DF = 5)      Random numbers (0 to 1) should have `gaps' where no random number isless than 0.5 in a predictable manner.   For this test, the number ofrandom numbers needed to give 1000 `gaps' will also vary a little.      For each batch of 1000 `gaps,' the program shows the batch number, thecount for the batch, the number of `gaps' of length zero (succeeding randomnumbers were less than 0.5), one through four, and five or more, then thebatch chi square and the cumulative chi square for all batches. 7.  Permutation Test   (DF = 5)      If three different random numbers (0 to 1) are created, there are sixdifferent possible permutations.   That is, LMH, LHM, MLH, MHL, HLM, andHML are equally likely, where L = low, M = medium, and H = high.      For each batch of 1000 sets of permutations, the program shows thebatch number, the bin count for each permutation, the batch chi square, andthe cumulative chi square for all batches. 8.  Mean Square Successive Difference Test      Successive random numbers (0 to 1) should have a predictable averagesquared difference.   This routine calculates both eta and the Z score.Eta is the sum of the squared differences of the successive iterationsdivided by the sum of the square of the differences between each iterationand the mean.   The Z score is then calculated from eta.   A positive Zscore indicates long trends, while a negative Z score indicates shortoscillations.      For each batch of 1000 differences, the program shows the batchnumber, eta, the Z score, and the cumulative Z score. 9.  Poker Test   (DF = 3)      A `hand' of four random integers (1 through 10) has a predictablenumber of duplications of integers.   The matches may be none (garbage),one (1 pair), two (2 pairs), three (3 of a kind), or four (4 of a kind).Since four matches occur so seldom, they are combined with three matchesfor the chi square calculation.      For each batch of 1000 `hands,' the program shows the batch number,the number of matches for each category, the batch chi square, and thecumulative chi square. 10.  Coupon Collector's Test      For random integers (1 through 10), to get at least one of eachinteger takes an average of 29.29 samples.   For each batch of 100 completesets, the program shows the batch number, the average number of samples forthe batch, and the cumulative average. 11.  Pascal's Triangle Test   (DF = 8)      For a ten-level Pascal's triangle, if a `dropping ball' has a 50%chance of being deflected right or left at each level (like a pin-ballmachine), the distribution of balls in the eleven bins at the bottom ispredictable (see Table 4).   Since there will be so few balls in the endbins, these are combined with the next-to-end bins for the chi squarecalculation.      For each batch of 1024 `balls', the program shows the batch number,the number of balls in each bin, the batch chi square, and the cumulativechi square. 12.  Maximum of `T' Test   (DF = 3)      For three successive random integers (1 through 10), each should bethe largest 28.5% of the time, with no maximum (duplication of the highinteger) 14.5% of the time.      For each batch of 1000 sets of three random integers, the programshows the batch number, the number of `no maximum,' the number of timesthat each of the three sampled integers was maximum, the batch chi square,and the cumulative chi square. 13.  Rectilinear Random Walk (Drunkard's Walk)      If walls are set up at X equals ñ20 and at Y equals ñ20, and a randomwalk started at X = 0 and Y = 0, on the average it should take 472 randomsteps (right or left, up or down) to hit a wall, with a standard deviationof 260 steps.   After 25 walks, the average number of steps should bebetween 420 and 524.   Also, X and Y, plus and minus walls should be hitequally often.      For each random walk, the program shows the walk number, the X and Yvalues after a wall is hit, the number of steps required for this walk, theaverage number of steps for all walks, and (after the first walk) thestandard deviation of the number of steps.   Note that this test is quitevariable. 14.  Looping Test      A random number generator, producing floating point random numbersbetween zero and one, should not reproduce its original random number, norany random number produced later.   If it does, it will loop (produce thesame sequence of numbers again).   This test checks the random numbers forduplication at user-specified intervals.   Up to 100 random numbers may bestored and checked against all new random numbers for duplication.   If aduplication occurs, the sequence numbers of the duplicate random numbersare displayed.   Otherwise, only the number of random numbers that havebeen checked is displayed.      Note that if one duplication is found, many others will follow, as thesame loop of numbers is generated again and again.   To test one millionrandom numbers, you must set the interval at least as large as 10000. Results      So, after running all of these RNG tests, how good are the MS-DOSBasic and compiled Basic random number generators?   First, a digression:      While Table 2 shows only six columns for two interpreted Basics andfour compiled Basics, more were tested.   BASICA version 1.19 and GWBasicsversions 2.0, 2.02, 3.22, and 3.23 gave identical numeric results, so the`GWBasic' column gives the results for all five.   With the compiledBasics, BASCOM and QuickBasic 1 gave identical numeric results, as shown inthe `BASCOM' column.   In addition, the QuickBasics 3, 4, 4.5, plus BasicCompiler 7 (sometimes called QuickBasic Extended) and Visual Basic for DOSall gave identical numeric results, all listed in the `QB 4.00' column.  Also,QB64 and the QBasic interpreter gave identical results, as listed in the`QBasic' column.      Back to how good are the random number generators.   As you can see inTable 2, that depends:   If you are using QB64 or the QBasic interpretedBasic, it passed all tests, at least up to one million numbers (10 millionfor QB64), though there are 4 tests for which the results were `questionable'.      However, random number generators in all other Basics failed at leastone test.   The worst was the old BASCOM compiler, which failed threetests, Run, Gap, and Poker.   The best was the Turbo Basic compiler, whichfailed only the Looping test.   If you need `only' 2.5 million or so randomnumbers, Turbo Basic would be best.   QuickBasic 4 and its clones failedtwo tests and was questionable on two others.   The Microway Basic compilerfailed only one test, and was questionable on one other test.      The two tests that were most difficult for the random numbergenerators to pass were the Poker test, failed by three RNGs, and thePascal's Triangle test, failed by two RNGs.      One other interesting finding was the wide difference in how fast theRNGs in the various Basic compilers worked.   Table 3 shows that QuickBasic 4 was far faster than seven of the other compilers, over ten times asfast as QuickBasic 4.5 or Turbo Basic.   The Microway compiler was nearlyas fast as QuickBasic 4.      Since most of the MS-DOS Basics were from Microsoft, the variety ofresults obtained is quite surprising.   I would have expected them to finda good random number generator and to then stick with it.   But no, aftergetting a pretty good RNG for QBasic, they had a couple of very poor onesfor their BASCOM and QuickBasic 1 compilers.   Then for QuickBasic 3 andlater, they apparently use the same, not too good RNG routine.   Odd, tosay the least.      I have heard that more computer time has been spent in testing RNGsthan in actually using them.   After evaluating these RNGs, I can believeit! Other Systems:      The `RNG' program is correct for running on most Basic interpreters orfor compiling on the most compilers.   One minor change was needed for itrun on the old BASCOM compiler:   Delete all CLS statements.   The BASCOMcompiler does not support CLS.   Also another minor change was used forrunning the program with the Microway compiler:   Add a statement to retaincertain variables in the numeric coprocessor's registers.   This speeds upoperation by about 12.5%. If Your RND Fails the Test:      If the RND function in your BASIC (or compiled BASIC) is not goodenough for your requirements, must you give up?   No, you can stillsimulate the RND function with a few BASIC statements.   The best formof RNG is called a `mixed congruential' RNG.   To implement this form ofRNG, you simply select a seed number `R' (which may be zero), multiplyit by a selected number, add a different selected number, then save onlythe decimal portion of the resulting number for both the random numberand for the seed number for the next random number.   That sounds rathercomplicated, but it's really quite simple:   Just use something likeP = 21 and Q = 0.3271, and the formula:                         R(n+1) = FRAC(R(n)*P+Q) A subroutine such as:                     9000 R=R*P+Q :R=R-INT(R) :RETURN might be used after setting the values of P, Q & R in the main program.Note that the values selected for P and Q are critical, and the greaterthe precision of calculation, the better.   You might try P = 9821,Q = 0.211327 and R = 0;  I tried these in the subroutine above in acompiled BASIC version of the RNG test program, and after one millionnumbers, it passed all but the Maximum of `T' test.   Of course, I useddouble precision for P, Q, and R to give enough significant digits inthe subroutine calculations.   Experiment, or check the references.   Ican't guarantee that these values will work for you, as each BASIC doesits calculations a little differently.   Good luck! Conclusions:      The subject of random numbers may seem simple, but it actually getsheavily into statistics.   By the use of these RNG tests, you may checkwhether your BASIC's RND function is good enough for you to base seriousprograms on it, or whether it can be used only for selecting a startingpoint for a game.   Even if your RND fails, you may be able to replaceit with a subroutine. References: Heyman, Victor K.  65 Notes   Vol. 4   No. 8   pp. 1 - 6 Knuth, Donald E.   The Art of Computer Programming,   Vol. 2,                   Seminumerical Algorithms,   pp. 1 - 160 Note: 65 Notes is now called:   The PPC Journal It is the publication of:   The PPC Club                            2541 W. Camden Place                            Santa Ana, Calif.  92704`

Code: [Select]
`10 CLS20 PRINT"-  -  -  -  -  -  Random Number Generator (RNG) Tests  -  -  -  -  -  -30 PRINT40  '              (c)  by Delmer D. Hinrichs  199750  '                      2116 S. E. 377th Ave.60  '                      Washougal, Wash. 9867170  '80  '   Written in MS-DOS GWBASIC for IBM PC/XT/AT/PS2 and compatibles90  '100 DEFINT I-N :S=.25 :H=.5 :T=.75 :B=0110 IF II=0 THEN DIM A(100), C(10), J(10), P(10) :II=1120 FOR I=0 TO 10 :C(I)=0 :J(I)=0 :P(I)=0 :NEXT I130 FOR I=0 TO 100 :A(I)=0 :NEXT I140 A=RND(-1)150 PRINT"Select the test that you want from the following list" :PRINT160 PRINT"     1.  Distribution Test"170 PRINT"     2.  Mean and Standard Deviation"180 PRINT"     3.  Correlation Test"190 PRINT"     4.  Serial Test"200 PRINT"     5.  Run Test"210 PRINT"     6.  Gap Test"220 PRINT"     7.  Permutation Test"230 PRINT"     8.  Mean Square Successive Difference Test"240 PRINT"     9.  Poker Test"250 PRINT"    10.  Coupon Collector's Test"260 PRINT"    11.  Pascal's Triangle Test"270 PRINT"    12.  Maximum of `T' Test"280 PRINT"    13.  Rectilinear Random Walk Test"290 PRINT"    14.  Loop Test"300 PRINT"     0.  Quit Tests" :PRINT310 PRINT TAB(15); :INPUT"Which Test No."; A\$ :A=VAL(A\$) :IF A\$="0" THEN END320 ON A GOTO 350,590,800,1040,1290,1570,1840,2170,2410,2750,2970,3290,3580,3810330 PRINT"Your Entry `"; A\$; "' is Illegal.  Try Again"; :GOTO 310340  '350 CLS360 PRINT"(1) *   *   *   *   *   Distribution Test for RNGs   *   *   *   *   *"370 PRINT"Bin-->  1    2    3    4    5    6    7    8    9   10      Chi  Cum.Chi"380 PRINT"Ideal  100  100  100  100  100  100  100  100  100  100     8.34   8.34"390 WHILE INKEY\$="" :WEND400 FOR J=1 TO 1000410   IR=INT(RND*10)+1 :J(IR)=J(IR)+1420   NEXT J430 X=0 :C=0 :B=B+1440 FOR I=1 TO 10450   X=X+(J(I)-100)^2/100460   C(I)=C(I)+J(I)470   C=C+(C(I)-100*B)^2/(100*B)480   NEXT I490 PRINT USING"#####  ###  ###  ###  ###  ###  ###  ###  ###  ###  ###   ###.## ###.##"; B; J(1); J(2); J(3); J(4); J(5); J(6); J(7); J(8); J(9); J(10); X; C500 FOR I=1 TO 10 :J(I)=0 :NEXT I510 IF INKEY\$<>"" THEN PRINT :GOTO 540520 IF B<1000 GOTO 570530   PRINT"1000 Trials"540   PRINT"Bin-->  1    2    3    4    5    6    7    8    9   10      Chi  Cum.Chi"550   PRINT"Ideal  100  100  100  100  100  100  100  100  100  100     8.34   8.34"560   WHILE INKEY\$="" :WEND :GOTO 10570 GOTO 400580  '590 CLS600 PRINT"(2)  *   *   *   *   Mean and Standard Deviation of RNGs   *   *   *   *"610 PRINT"                  Mean     Std.D     Cum.Mean  Cum.SD"620 PRINT"Ideal            0.5000   0.2887      0.5000   0.2887"630 WHILE INKEY\$="" :WEND :C=0 :D=0640 R1=0 :R2=0650 FOR I=1 TO 1000660   R=RND :R1=R1+R :R2=R2+R*R670   NEXT I680 B=B+1 :T=1000*B :C=C+R1 :D=D+R2690 SD=SQR((R2-R1*R1/1000)/999)700 CD=SQR((D-C*C/T)/(T-1))710 PRINT USING"#####            #.####   #.####      #.####   #.####"; B, R1/1000, SD, C/T, CD720 IF INKEY\$<>"" THEN PRINT :GOTO 750730 IF B<1000 GOTO 780740   PRINT"1000 Trials"750   PRINT"                  Mean     Std.D     Cum.Mean  Cum.SD"760   PRINT"Ideal            0.5000   0.2887      0.5000   0.2887"770   WHILE INKEY\$="" :WEND :GOTO 10780 GOTO 640790  '800 CLS810 PRINT"(3)  *   *   *   *   *   Correlation Test for RNGs   *   *   *   *   *"820 PRINT, " Cor.Coef", " Z Score", "Cum. Z Score"830 PRINT"Ideal", " 0.000000", " 0.000000", " 0.000000"840 WHILE INKEY\$="" :WEND :C1=0 :C2=0 :D1=0 :D2=0 :CD=0850 X1=0 :X2=0 :Y1=0 :Y2=0 :XY=0 :X=RND860 FOR I=1 TO 1000870   Y=X :X=RND880   X1=X1+X :X2=X2+X*X :Y1=Y1+Y :Y2=Y2+Y*Y :XY=XY+X*Y890   NEXT I900 B=B+1 :C1=C1+X1 :C2=C2+X2 :D1=D1+Y1 :D2=D2+Y2 :CD=CD+XY910 C=(1000*XY-X1*Y1)/SQR((1000*X2-X1*X1)*(1000*Y2-Y1*Y1))920 Z=(C+1/999)/((1/999)*SQR(997000!/1001))930 W=1000*B :CC=(W*CD-C1*D1)/SQR((W*C2-C1*C1)*(W*D2-D1*D1))940 CZ=(CC+1/(W-1))/((1/(W-1))*SQR(W*(W-3)/(W+1)))950 PRINT USING"#####         ##.######     ##.######     ##.######";B, C, Z, CZ960 IF INKEY\$<>"" THEN PRINT :GOTO 990970 IF B<1000 GOTO 1020980   PRINT"1000 Trials"990    PRINT, " Cor.Coef", " Z Score", "Cum. Z Score"1000   PRINT"Ideal", " 0.000000", " 0.000000", " 0.000000"1010   WHILE INKEY\$="" :WEND :GOTO 101020 GOTO 8501030  '1040 CLS1050 PRINT"(4)   *   *   *   *   *   Serial Test for RNGs   *   *   *   *   *   *"1060 PRINT"Pair->    1    2    3    4    5    6    7    8    9      Chi  Cum.Chi"1070 PRINT"Ideal    111  111  111  111  111  111  111  111  111     7.34   7.34"1080 WHILE INKEY\$="" :WEND1090 FOR I=1 TO 10001100   K=3*INT(3*RND) + INT(3*RND)1110   J(K)=J(K)+11120   NEXT I1130 B=B+1 :X=0 :C=01140 FOR I=0 TO 81150   X=X+(J(I)-111)^2/1111160   C(I)=C(I)+J(I)1170   C=C+(C(I)-111*B)^2/(111*B)1180   NEXT I1190 PRINT USING"#####    ###  ###  ###  ###  ###  ###  ###  ###  ###   ###.## ###.##"; B, J(0), J(1), J(2), J(3), J(4), J(5), J(6), J(7), J(8), X, C1200 FOR I=0 TO 8 :J(I)=0 :NEXT I1210 IF INKEY\$<>"" THEN PRINT :GOTO 12401220 IF B<1000 GOTO 12701230   PRINT"1000 Trials"1240   PRINT"Pair->    1    2    3    4    5    6    7    8    9      Chi  Cum.Chi"1250   PRINT"Ideal    111  111  111  111  111  111  111  111  111     7.34   7.34"1260   WHILE INKEY\$="" :WEND :GOTO 101270 GOTO 10901280  '1290 CLS1300 PRINT"(5)  *   *   *   *   *   *   Run Test for RNGs   *   *   *   *   *   *"1310 PRINT"         Count        0    1    2   3    4+      Chi    Cum.Chi"1320 PRINT"Ideal    2718        500  333  125  33   8       3.36     3.36"1330 P(0)=500 :P(1)=333.3 :P(2)=125 :P(3)=33.33 :P(4)=8.3331340 WHILE INKEY\$="" :WEND :M=01350 FOR I=1 TO 10001360   D=RND :K=01370   R=RND :IF R>D THEN D=R :K=K+1 :GOTO 13701380   IF K>4 THEN K=41390   J(K)=J(K)+1 :M=M+K+11400   NEXT I1410 B=B+1 :X=0 :C=01420 FOR I=0 TO 41430   X=X+(J(I)-P(I))^2/P(I)1440   C(I)=C(I)+J(I)1450   C=C+(C(I)-P(I)*B)^2/(P(I)*B)1460   NEXT I1470 PRINT USING"#####    ####        ###  ###  ###  ##  ##      ##.##    ##.##"; B, M+1000, J(0), J(1), J(2), J(3), J(4), X, C1480 FOR I=0 TO 4 :J(I)=0 :NEXT I1490 IF INKEY\$<>"" THEN PRINT :GOTO 15201500 IF B<1000 GOTO 15501510   PRINT"1000 Trials"1520   PRINT"         Count        0    1    2   3    4+      Chi    Cum.Chi"1530   PRINT"Ideal    2718        500  333  125  33   8       3.36     3.36"1540   WHILE INKEY\$="" :WEND :GOTO 101550 M=0 :GOTO 13501560  '1570 CLS1580 PRINT"(6)  *   *   *   *   *   *   Gap Test for RNGs   *   *   *   *   *   *"1590 PRINT"        Count      0    1    2   3   4   5+      Chi    Cum.Chi"1600 PRINT"Ideal   2000      500  250  125  63  31  31      4.35     4.35"1610 P(0)=500 :FOR I=1 TO 4 :P(I)=P(I-1)/2 :NEXT I :P(5)=P(4)1620 WHILE INKEY\$="" :WEND :K=0 :L=01630 FOR I=1 TO 10001640   K=K+1 :IF RND>H THEN L=L+1 :GOTO 16401650   IF L>5 THEN L=51660   J(L)=J(L)+1 :L=01670   NEXT I1680 B=B+1 :X=0 :C=01690 FOR I=0 TO 51700   X=X+(J(I)-P(I))^2/P(I)1710   C(I)=C(I)+J(I)1720   C=C+(C(I)-P(I)*B)^2/(P(I)*B)1730   NEXT I1740 PRINT USING"#####   ####      ###  ###  ### ### ### ###     ##.##    ##.##"; B, K, J(0), J(1), J(2), J(3), J(4), J(5), X, C1750 FOR I=0 TO 5 :J(I)=0 :NEXT I1760 IF INKEY\$<>"" THEN PRINT :GOTO 17901770 IF B<1000 GOTO 18201780   PRINT"1000 Trials"1790   PRINT"        Count      0    1    2   3   4   5+      Chi    Cum.Chi"1800   PRINT"Ideal   2000      500  250  125  63  31  31      4.35     4.35"1810   WHILE INKEY\$="" :WEND :GOTO 101820 K=0 :GOTO 16301830  '1840 CLS1850 PRINT"(7)  *   *   *   *   *   Permutation Test for RNGs   *   *   *   *   *"1860 PRINT"       Bin-->     1    2    3    4    5    6       Chi    Cum.Chi"1870 PRINT"Ideal            167  167  167  167  167  167      4.35     4.35"1880 P=166.67 :K=01890 WHILE INKEY\$="" :WEND1900 FOR I=1 TO 10001910   D=RND :E=RND :F=RND1920   IF D>E GOTO 19501930   IF F>D GOTO 19801940   K=2 :GOTO 19901950   IF F>D THEN K=3 :GOTO 19901960   IF E>F THEN K=0 ELSE K=11970   GOTO 19901980   IF E>F THEN K=4 :ELSE K=51990   J(K)=J(K)+12000   NEXT I2010 B=B+1 :X=0 :C=02020 FOR I=0 TO 52030   X=X+(J(I)-P)^2/P2040   C(I)=C(I)+J(I)2050   C=C+(C(I)-B*P)^2/(B*P)2060   NEXT I2070 PRINT USING"#####            ###  ###  ###  ###  ###  ###     ##.##    ##.##"; B, J(0), J(1), J(2), J(3), J(4), J(5), X, C2080 FOR I=0 TO 5 :J(I)=0 :NEXT I :K=02090 IF INKEY\$<>"" THEN PRINT :GOTO 21202100 IF B<1000 GOTO 21502110   PRINT"1000 Trials"2120   PRINT"       Bin-->     1    2    3    4    5    6       Chi    Cum.Chi"2130   PRINT"Ideal            167  167  167  167  167  167      4.35     4.35"2140   WHILE INKEY\$="" :WEND :GOTO 102150 GOTO 19002160  '2170 CLS2180 PRINT"(8)  *   *   *   Mean Square Successive Difference RNG Test   *   *   *"2190 PRINT, "   Eta", " Z Score", "Cum. Z Score"2200 PRINT"Ideal", " 2.000000", " 0.000000", " 0.000000"2210 WHILE INKEY\$="" :WEND :CD=0 :C1=0 :C2=02220 X=RND :D2=0 :X1=0 :X2=02230 FOR I=1 TO 10002240   Y=X :X=RND2250   D2=D2+(Y-X)*(Y-X) :X1=X1+X :X2=X2+X*X2260   NEXT I2270 B=B+1 :W=1000*B :CD=CD+D2 :C1=C1+X1 :C2=C2+X22280 ETA=D2/(X2-X1*X1/1000)2290 Z=(1-ETA/2)/SQR(998/999999!)2300 CE=CD/(C2-C1*C1/W)2310 CZ=(1-CE/2)/SQR((W-2)/(W*W-1))2320 PRINT USING"#####         ##.######     ##.######     ##.######"; B, ETA, Z, CZ2330 IF INKEY\$<>"" THEN PRINT :GOTO 23602340 IF B<1000 GOTO 23902350   PRINT"1000 Trials"2360   PRINT, "   Eta", " Z Score", "Cum. Z Score"2370   PRINT"Ideal", " 2.000000", " 0.000000", " 0.000000"2380   WHILE INKEY\$="" :WEND :GOTO 102390 GOTO 22202400  '2410 CLS2420 PRINT"(9)    *   *   *   *   *   Poker Test for RNGs   *   *   *   *   *   *"2430 PRINT"     Matchs-->    0    1   2   3   4              Chi   Cum.Chi"2440 PRINT"Ideal            504  432  27  36  1              2.37    2.37"2450 P(0)=504 :P(1)=432 :P(2)=27 :P(3)=372460 WHILE INKEY\$="" :WEND2470 FOR I=1 TO 1000 :K=02480   J1=INT(RND*10)+1 :J2=INT(RND*10)+1 :J3=INT(RND*10)+1 :J4=INT(RND*10)+12490   IF J1=J2 THEN K=12500   IF J1=J3 THEN K=K+12510   IF J1=J4 THEN K=K+12520   IF J2=J3 THEN K=K+12530   IF J2=J4 THEN K=K+12540   IF J3=J4 THEN K=K+12550   IF K>4 THEN K=42560   J(K)=J(K)+12570   NEXT I2580 B=B+12590 J(3)=J(3)+J(4) :X=0 :C=02600 FOR I=0 TO 32610   X=X+(J(I)-P(I))^2/P(I)2620   C(I)=C(I)+J(I)2630   C=C+(C(I)-P(I)*B)^2/(P(I)*B)2640   NEXT I2650 PRINT USING"#####            ###  ### ### ### ##             ##.##   ##.##"; B, J(0), J(1), J(2), J(3), J(4), X, C2660 FOR I=0 TO 4 :J(I)=0 :NEXT I2670 IF INKEY\$<>"" THEN PRINT :GOTO 27002680 IF B<1000 GOTO 27302690   PRINT"1000 Trials"2700   PRINT"     Matchs-->    0    1   2   3   4              Chi   Cum.Chi"2710   PRINT"Ideal            504  432  27  36  1              2.37    2.37"2720   WHILE INKEY\$="" :WEND :GOTO 102730 GOTO 24702740  '2750 CLS2760 PRINT"(10)  *   *   *   *   Coupon Collector's Test for RNGs   *   *   *   *"2770 PRINT, "Average", "Cum.Avg"2780 PRINT"Ideal", 29.29, 29.292790 WHILE INKEY\$="" :WEND :C=0 :M=02800 FOR I=1 TO 1002810   FOR L=1 TO 102820     K=INT(RND*10)+1 :M=M+1 :IF J(K)=1 GOTO 28202830     J(K)=12840     NEXT L2850   FOR K=1 TO 10 :J(K)=0 :NEXT K2860   NEXT I2870 B=B+1 :C=C+M2880 PRINT USING"#####         ###.##        ###.##"; B, M/100, C/(100*B)2890 IF INKEY\$<>"" THEN PRINT :GOTO 29202900 IF B<1000 GOTO 29502910   PRINT"1000 Trials"2920   PRINT, "Average", "Cum.Avg"2930   PRINT"Ideal", 29.29, 29.292940   WHILE INKEY\$="" :WEND :GOTO 102950 M=0 :GOTO 28002960  '2970 CLS2980 PRINT"(11)   *   *   *   Pascal's Triangle Test for RNGs   *   *   *   *   *"2990 PRINT"Bin->    1   2   3   4    5    6    7    8    9  10  11    Chi   Cum.Chi"3000 PRINT"Ideal    1  10  45  120  210  252  210  120  45  10  1     7.34    7.34"3010 P(1)=11 :P(2)=45 :P(3)=120 :P(4)=210 :P(5)=2523020 P(6)=210 :P(7)=120 :P(8)=45 :P(9)=113030 WHILE INKEY\$="" :WEND3040 FOR I=1 TO 1024 :K=103050   FOR L=1 TO 103060     R=RND3070     IF R>H THEN K=K+1 ELSE K=K-13080     NEXT L :J(K/2)=J(K/2)+13090   NEXT I3100 B=B+13110 J0=J(0) :J1=J(1) :J9=J(9) :JT=J(10)3120 J(1)=J(1)+J(0) :J(9)=J(9)+J(10)3130 X=0 :C=03140 FOR I=1 TO 93150   X=X+(J(I)-P(I))^2/P(I)3160   C(I)=C(I)+J(I)3170   C=C+(C(I)-P(I)*B)^2/(P(I)*B)3180   NEXT I3190 PRINT USING"#####   ##  ##  ##  ###  ###  ###  ###  ###  ##  ## ##    ##.##   ##.##"; B, J0, J1, J(2), J(3), J(4), J(5), J(6), J(7), J(8), J9, JT, X, C3200 FOR I=0 TO 10 :J(I)=0 :NEXT I3210 IF INKEY\$<>"" THEN PRINT :GOTO 32403220 IF B<1000 GOTO 32703230   PRINT"1000 Trials"3240   PRINT"Bin->    1   2   3   4    5    6    7    8    9  10  11    Chi   Cum.Chi"3250   PRINT"Ideal    1  10  45  120  210  252  210  120  45  10  1     7.34    7.34"3260   WHILE INKEY\$="" :WEND :GOTO 103270 GOTO 30403280  '3290 CLS3300 PRINT"(12)  *   *   *   *   Maximum of `T' Test for RNGs   *   *   *   *   *"3310 PRINT"     Max.->     None   1    2    3                Chi   Cum.Chi"3320 PRINT"Ideal            145  285  285  285               2.37    2.37"3330 P(0)=145 :P(1)=285 :P(2)=285 :P(3)=2853340 WHILE INKEY\$="" :WEND3350 FOR I=1 TO 10003360   J=INT(RND*10)+1 :K=INT(RND*10)+1 :L=INT(RND*10)+13370   IF J>K THEN IF J>L THEN J(1)=J(1)+13380   IF K>J THEN IF K>L THEN J(2)=J(2)+13390   IF L>J THEN IF L>K THEN J(3)=J(3)+13400   NEXT I3410 J(0)=1000-J(1)-J(2)-J(3)3420 B=B+1 :X=0 :C=03430 FOR I=0 TO 33440   X=X+(J(I)-P(I))^2/P(I)3450   C(I)=C(I)+J(I)3460   C=C+(C(I)-P(I)*B)^2/(P(I)*B)3470   NEXT I3480 PRINT USING"#####            ###  ###  ###  ###              ##.##   ##.##"; B, J(0), J(1), J(2), J(3), X, C3490 FOR I=0 TO 3 :J(I)=0 :NEXT I3500 IF INKEY\$<>"" THEN PRINT :GOTO 35303510 IF B<1000 GOTO 35603520   PRINT"1000 Trials"3530   PRINT"     Max.->     None   1    2    3                Chi   Cum.Chi"3540   PRINT"Ideal            145  285  285  285               2.37    2.37"3550   WHILE INKEY\$="" :WEND :GOTO 103560 GOTO 33503570  '3580 CLS3590 :PRINT"(13)   *   *   *   Rectilinear Random Walk RNG Test   *   *   *   *   *"3600 PRINT"                 X and Y        Steps    Avg.      Std.Dev"3610 PRINT"Ideal            Values          472    472.00     260.00"3620 WHILE INKEY\$="" :WEND :C=0 :C2=03630 D=0 :J=0 :K=03640 R=RND :D=D+1 :IF R>H GOTO 36803650 IF R>S THEN J=J+1 ELSE J=J-13660 IF ABS(J)<20 GOTO 36403670 GOTO 37003680 IF R>T THEN K=K+1 ELSE K=K-13690 IF ABS(K)<20 GOTO 36403700 B=B+1 :C=C+D :C2=C2+D*D3710 IF B>1 THEN SD=SQR((C2-C*C/B)/(B-1))3720 PRINT USING"#####          ###   ###        ####   ####.##     ###.##"; B, J, K, D, C/B, SD3730 IF INKEY\$<>"" THEN PRINT :GOTO 37603740 IF B<1000 GOTO 37903750   PRINT"1000 Trials"3760   PRINT"                 X and Y        Steps    Avg.      Std.Dev"3770   PRINT"Ideal            Values          472    472.00     260.00"3780   WHILE INKEY\$="" :WEND :GOTO 103790 GOTO 36303800  '3810 CLS3820 PRINT"(14)   *   *   *   *   *   *   RNG Looping Test   *   *   *   *   *   *"3830 R=RND :A(0)=R3840 PRINT :INPUT"Test `RND' at Intervals of (10000 for 1 Million RNDs) "; W3850 IF W>0 THEN IF W=INT(W) GOTO 38703860 PRINT"Your Entry "; W;" is Illegal.  Try Again" :GOTO 38403870 PRINT :PRINT"Number of RNDs Tested = "3880 FOR I=1 TO 100 :J=I-1 :P=W*J :PRINT P;3890   FOR U=1 TO W3900     R=RND3910     FOR K=0 TO J3920       IF R=A(K) THEN PRINT :PRINT :PRINT"Looping Test FAILED.  RND #"; P+U;           "Duplicates"; K*W :END3930       NEXT K3940     NEXT U3950   A(I)=R3960   NEXT I3970 PRINT W*100 :PRINT :PRINT"Looping Test Passed."3980 WHILE INKEY\$="" :WEND :GOTO 103990 END `
Code: [Select]
` Table 2.   Results of RNG Tests with various Basics and Compiled Basics.          Basic Interpreters                Basic Compilers         __________________     ___________________________________________         GWBasic     QBasic      BASCOM     Turbo      QB 4.00     Microway         _______    _______     ________   ________    ________    ________ 1.  Distribution     DF = 9    Ideal Chi = 8.34   (1000)          7.21       1 _7. _03 _       4.54       13.98       4.40        10.49 2.  Mean    Ideal Mean = 0.5000    (1000)         0.5002     0.5001      0.5000     0.4999      0.4999      0.5001     Standard Deviation     Ideal Std.Dev. = 0.2887   (1000)         0.2886     0.2888      0.2887     0.2887      0.2887      0.2887 3.  Correlation   Ideal Z-Score = 0.0000   (999)        -0.4326    -0.8499      0.4036     0.3462     -0.7113     -0.2695 4.  Serial   DF = 8   Ideal Chi = 7.34   (500)          10.55      5.11        10.83      15.09       10.76       14.90 5.  Run    DF = 4    Ideal Chi = 3.36    (368)          3.54        _0. _45 _        0 _. _1 _9 _        3.66        0.91        2.48 6.  Gap    DF = 5    Ideal Chi = 4.35    (500)          1.35        _0. _94 _        0 _. _3 _4 _        5.15         _0. _81 _        4.90 7.  Permutation    DF = 5    Ideal Chi = 4.35    (333)          1.26        _1. _10 _        2.87        3.76         _0. _68 _        1 _3. _12 _ 8.  Mean Square Successive Difference    Ideal Z-Score = 0.0000    (1000)        -0.4311    -0.8334      0.3736      0.3909      -0.6978    -0.3600 9.  Poker    DF = 3    Ideal Chi = 2.37    (250)          1 _1 _. _4 _4 _      4.31        9 _8 _. _0 _7 _       4.43        9 _5 _. _5 _0 _       1.33 10. Coupon Collector    Ideal Average = 29.29    (1000)          29.31      29.32       29.28       29.28       29.29       29.30 11. Pascal's Triangle    DF = 8    Ideal Chi = 7.34    (98)          4.57        6.09        4.96        7.85        1 _. _0 _5 _        0 _. _7 _3 _ 12. Maximum of `T'    DF = 3    Ideal Chi = 2.37    (333)          1.12        1.28        0.80        0.70        0.46        6.08 13. Random Walk, Average    Ideal Average = 472    (2119)          471         465         469          475         470        474     Standard Deviation    Ideal Std.Dev. = 260    (2119)          327         328         329          336         333        333 14. Looping    Ideal Looping = Never    (1,000,000)       1 Million   1 Million  10 Million    F _a _i _l _e _d _!  10 Million  10 Million                                         [2,653,368=                                           400,000] Values in parentheses are number of batches to test one million random Nos. The underlined values are unsatisfactory;  the RNG failed that test. The half-underlined values are questionable;  the RNG gave unlikely results`

Code: [Select]
`             Table 1.   ChiﾄSquare Distribution Table            ﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍﾍ         P=99%   P=95%   P=75%   P=50%   P=25%    P=5%    P=1%      ﾄﾄﾄﾄﾄﾄﾄ ﾄﾄﾄﾄﾄﾄﾄ ﾄﾄﾄﾄﾄﾄﾄ ﾄﾄﾄﾄﾄﾄﾄ ﾄﾄﾄﾄﾄﾄﾄ ﾄﾄﾄﾄﾄﾄﾄ ﾄﾄﾄﾄﾄﾄﾄ DF=1  0.00016 0.00393  0.1015  0.4549  1.323   3.841   6.635 DF=2  0.02010 0.1026   0.5753  1.386   2.773   5.991   9.210 DF=3  0.1148  0.3518   1.213   2.366   4.108   7.815  11.34 DF=4  0.2971  0.7107   1.923   3.357   5.385   9.488  13.28 DF=5  0.5543  1.1455   2.675   4.351   6.626  11.07   15.09 DF=6  0.8720  1.635    3.455   5.348   7.841  12.59   16.81 DF=7  1.239   2.167    4.255   6.346   9.037  14.07   18.48 DF=8  1.646   2.733    5.071   7.344  10.22   15.51   20.09 DF=9  2.088   3.325    5.899   8.343  11.39   16.92   21.67 DF=10 2.558   3.940    6.737   9.342  12.55   18.31   23.21 DF=11 3.053   4.575    7.584  10.34   13.70   19.68   24.73 DF=12 3.571   5.226    8.438  11.34   14.84   21.03   26.22         "DF" is the Degrees of Freedom.        "P" is the Probability.   50% is the "ideal" value,       showing that the scatter from the most likely result       is just what would be expected.   Greater than 50%       (low chi square values) shows scatter from the most       likely result to be less than expected, while less       than 50% (high chi square values) shows scatter to be       greater than expected.   Values outside the 1% to 99%       range should be rejected, while values outside the 5%       to 95% range are suspect.`

Code: [Select]
`            Table 5.   Start of Distribution Test for GWBasic  (1) *   *   *   *   *   Distribution Test for RNGs   *   *   *   *   *Bin-->  1    2    3    4    5    6    7    8    9   10      Chi  Cum.ChiIdeal  100  100  100  100  100  100  100  100  100  100     8.34   8.34    1  105  103  107   93  100   90   96   98  112   96     4.12   4.12    2   98   82   98   97  118   95  102  105  105  100     7.44   6.13    3   87   91  108  114   85   94  103  102  110  106     9.20   6.90    4   97  109  104  122   81   87  108   98  100   94    12.24   9.01    5  102  106  102   99   98  100   93  102  104   94     1.54   7.52    6  103  106   92   91  100  109   98  103   97  101     2.94   3.88    7   81  101   94  118  101  107  104   98  105   91     8.98   5.70    8  112  101   77  104  113  124   95   89   90   95    17.06   3.91    9  116   95   97   88   80   88   94  115  110  117    16.28   3.38   10   83  112  102   79  106  116   87  101  102  112    14.88   2.88   11  111  100  104   99  108   89  108  100   82   99     7.12   0.89   12   95   94  104  105  110   97  101  104  112   78     8.56   1.49   13  107   95   95  113  100   96   99  101   90  104     4.02   1.36   14   87  100   86   93  113  101   91  105  116  108    10.10   2.68   15  105   85   91  108   90  113  110  101   88  109     9.90   2.51`

Code: [Select]
`                                             1                               1     1                            1     2     1                         1     3     3     1                      1     4     6     4     1                   1     5    10    10     5     1                1     6    15    20    15     6     1             1     7    21    35    35    21     7     1          1     8    28    56    70    56    28     8     1       1     9    36    84   126   126    84    36     9     1    1    10    45   120   210   252   210   120    45    10     1     Table 4.   Pascal's Triangle  --  Ten-Level.   Each number is the sum         of the two numbers directly above it.   Also, the sum of the         numbers at each level is twice the sum of the numbers of the         preceding level.   If the triangle is considered as a pinball         machine, the numbers represent the relative probabilities of         a ball ending up in that bin (at any level).`
[code]

Table 3.   RNG Program Speed of Operation with Various Basic Compilers

Time, sec., to run 10,000 Batches
Basic Compiler, version and date     of #8, Mean Square Success. Difference
_________________________________    ______________________________________

BASCOM, ver. 5.35, 17 Mar 83                        244

Quick Basic 1, ver. 1.02, 26 Jan 86                 230.3

Quick Basic 3, ver. 3.00, 4 Jul 87                  108.7

Quick Basic 4, ver. 4.00, 8 Oct 87                   37.6

Quick Basic 4.5, ver. 4.50, 28 Sept 88              416.4

Basic Compiler 7, ver. 7.00, 9 Nov 89               298.7

Visual Basic for DOS, ver. 1.00, 19 Aug 92          374.4

Turbo Basic, ver. 1.0, 4 Feb 87

#### mcalkins

• Hero Member
•     • Posts: 1269
• •  ##### Re: FAQ suggestion: RND
« Reply #4 on: February 21, 2013, 11:14:44 PM »
Clippy has suggested that I put this in the Tutorials section, but I didn't write it as a tutorial. It's not really in the style of a tutorial. I'll have to consider it...

SMcNeill: Thanks for the contribution. However:

I would prefer to link to implementations of well known algorithms, such as the ones listed on the Wikipedia pages. These should have been reviewed by many people by now, and their properties should be understood.

Yours uses RND and TIMER to shuffle and occasionally reshuffle an array. I don't know how much this improves the randomness quality, or if the improvement is worth the overhead. It is not cryptographically secure.

MrWhy: QBASIC does not predate QuickBASIC 1. QBASIC was derived from QuickBASIC 4.x. Wikipedia says from 4.5, but it might have been from 4.0.

The Wikipedia article has, for Visual BASIC <=6:
m = 2^24
a = &h43FD43FD
c = &hC39EC3

Which is effectively identical to QBASIC. Wikipedia gives as a reference:
http://support.microsoft.com/kb/231847
which indicates that &h50000 is also the initial state in Visual BASIC <=6.

http://support.microsoft.com/kb/28150/EN-US
Which gives the following values:
m = 2^24
a = &h343fd
c = &h269ec3

And seems to imply that those numbers are used for more versions of BASIC than they actually are. However, I can confirm that those are the numbers in GW-BASIC 3.23, and that the initial state is &h4fc752.

Regards,
Michael
The QBASIC Forum Community: http://www.network54.com/index/10167 Includes off-topic subforums.
QB64 Off-topic subforum: http://qb64offtopic.freeforums.org/

#### Clippy

• Hero Member
•     • • Posts: 16439
• I LOVE π = 4 * ATN(1)    Use the QB64 WIKI >>>
• •  ##### Re: FAQ suggestion: RND
« Reply #5 on: February 22, 2013, 05:09:59 AM »
Put it in there anyway. Then Steve can only vote on it...

That's a good place to put stuff like that because I can link it in the WIKI without other commentary. We could also put special Linus or Mac install procedures in there.
QB64 WIKI: Main Page

#### mcalkins

• Hero Member
•     • Posts: 1269
• •  ##### Re: FAQ suggestion: RND
« Reply #6 on: February 24, 2013, 07:37:42 AM »
I want to research and probably implement some of the other algorithms...

(Wikipedia is reminding me how much I suck at math...)

Regards,
Michael
The QBASIC Forum Community: http://www.network54.com/index/10167 Includes off-topic subforums.
QB64 Off-topic subforum: http://qb64offtopic.freeforums.org/