• Print

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

mcalkins

  • Hero Member
  • *****
  • Posts: 1269
    • qbasicmichael.com
    • Email
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/Hexadecimal
http://en.wikipedia.org/wiki/Bitwise_operation
http://www.qb64.net/forum/index.php?topic=5990.0

This article follows the QBASIC tradition of using the &h prefix for hexadecimal numbers.
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]
CLS
PRINT "&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 &HFFFFFF
RND = 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]
CLS
PRINT 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 &HFFFFFF
state = &HFD43FD * state + &HC39EC3 AND &HFFFFFF
RND = 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
hexadecimal: c10c0000

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 OFF
PRINT "'"; 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]
CLS
PRINT "&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 &HFFFFFF
state = &HFD43FD * state + &HC39EC3 AND &HFFFFFF
RND = 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]
CLS
PRINT 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 OFF
PRINT "'"; MKD$(3.5); "'"
END

Displays: '      ♀@' (blanks are nulls.)

The QBASIC program:

Code: [Select]
CLS
PRINT "&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]
CLS
RANDOMIZE 3.5
PRINT "&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          fraction
0 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]
CLS
PRINT "&h"; HEX$(RND * &H1000000)
RANDOMIZE 5.25
PRINT "&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]
CLS
PRINT "&h"; HEX$(RND * &H1000000)
RANDOMIZE 5.25
PRINT "&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 = &H50000

test: END

FUNCTION qRnd!
qState = &HFD43FD * qState + &HC39EC3 AND &HFFFFFF
qRnd = qState / &H1000000
END FUNCTION

SUB qRandomizeUsing (s AS DOUBLE)
qState = CVL(RIGHT$(MKD$(s), 4))
qState = ((&HFFFF& AND qState) XOR qState \ &H10000) * &H100
END SUB

SUB qRandomize (s AS DOUBLE)
DIM l AS _UNSIGNED LONG
l = CVL(RIGHT$(MKD$(s), 4))
qState = ((&HFFFF& AND l) XOR l \ &H10000) * &H100 OR &HFF AND qState
END SUB

FUNCTION qRndPrevious!
qRndPrevious = qState / &H1000000
END FUNCTION

FUNCTION qRndSeed! (s AS SINGLE)
DIM l AS _UNSIGNED LONG
l = 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 &HFFFFFF
qRndSeed = qState / &H1000000
END FUNCTION

SUB test
DIM i AS LONG
PRINT "initial default seed": GOSUB p
RANDOMIZE 5.25: qRandomize 5.25
PRINT "randomize": GOSUB p
RANDOMIZE USING -3.5: qRandomizeUsing -3.5
PRINT "randomize using": GOSUB p
PRINT "rnd(0)"
PRINT RND(0), qRndPrevious
PRINT RND(0), qRndPrevious
PRINT "rnd(negative)"
PRINT RND(-8.75), qRndSeed(-8.75)
GOSUB p
EXIT SUB
p:
FOR i = 0 TO 2: PRINT RND, qRnd: NEXT
RETURN
END SUB


What if I need a better algorithm?

Read the following articles:

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.


Links.

http://en.wikipedia.org/wiki/Linear_congruential_generator
http://www.network54.com/Forum/178565/message/1046721807/
http://www.network54.com/Forum/178387/message/1046747461/
See also libqbx.cpp: rnd_seed, sub_randomize, func_rnd

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
    • Email
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 DOUBLE

FOR i = 1 TO Limit
    x = INT(Rndm * 10)
    spread(x) = spread(x) + 1
    IF i < 31 THEN PRINT x,
NEXT

PRINT: PRINT
FOR i = 0 TO 9
    PRINT spread(i),
NEXT


FUNCTION Rndm
STATIC init, count


IF 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)
    NEXT
END IF

count = count + 1
IF count > Limit THEN
    RANDOMIZE TIMER
    count = 0
    FOR i = 0 TO Limit
        x = INT(RND * Limit)
        SWAP PRNG(i), PRNG(x)
    NEXT
END IF
Rndm = 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 >>>
    • Pete's Qbasic Site
    • Email
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 ....  ;D
QB64 WIKI: Main Page
Download Q-Basics Code Demo: Q-Basics.zip
Download QB64 BAT, IconAdder and VBS shortcuts: QB64BAT.zip
Download QB64 DLL files in a ZIP: Program64.zip

Mrwhy

  • Hero Member
  • *****
  • Posts: 2893
  • My Dad called me Mr Why when I was 5.
    • Email
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 question
must be `No'. The Basic RND function generates each new `random' number
from the previous number, starting from an initial `seed' number, by a
simple calculation.   Thus these apparently random numbers are more
correctly called `pseudo-random' numbers.   For simplicity, from now on we
will just call them `random' numbers.   Which still leaves us with the
question:   Are these apparently random numbers close enough to being truly
random to be useful?
 
     But useful for what?   If you use the RND function merely to get a
different starting point for a Lunar Lander game, it doesn't make much
difference how random it is.   If you use RND to simulate dealing cards for
a poker game, you may lose more games than you should if the RND is not
random, but it still isn't really serious.   However, if you use RND in a
Monte Carlo simulation to determine the pricing policy for your business,
or for some other serious use, you must make sure that your RND is
adequately random before you can trust the results.   For these critical
uses, or even to check on the honesty of your poker dealer, there are tests
that 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 (more
truly random), the statistics of a large number of these random numbers
become more certain.   Therefore we can let our RND function produce a
large number of random numbers and check their statistics.
 
     Which statistics, and how do we check them?   There are many different
statistics that could be used, as random numbers have many different random
characteristics.   The `RNG' test program has fourteen different routines
that test for sixteen different random number characteristics.   A truly
random number generator (RNG) would pass all of these tests.   Any one RND
function may pass some tests easily, but fail others.   The uses to which
you 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 really
trust our unaided judgement on this, as a series of numbers that appear
random may not be random, and vice versa.   For example, if we throw a
single die, each of the six faces should come up about one-sixth of the
time.   If we throw the die 600 times, would we expect each face to come up
exactly 100 times?   I think that we can agree that this is highly
unlikely.   But just how much should the number of times that each face
comes up differ from 100?   One of the most useful statistical tests for
this type of problem is the chi square test.
 
     To apply the chi square test to our die-throwing experiment, we
compute:
 
                        (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 the
chi square, we need a chi square distribution table such as Table 1.
There we find another new term, `degrees of freedom.' This is simply the
number of independent choices possible, or one less than the number of
possible categories.   For throwing a die, there are five degrees of
freedom (if the face is not numbers 1 through 5, then it must be six, so
there are only 5 independent choices).
 
     In Table 1, the 50% probability column is the `Ideal' chi square
level, with the results neither more uniform (probability > 50%, too small
a chi square) nor more scattered (probability < 50%, too large a chi
square) than would be obtained from truly random events.   If each of the
six faces of the die came up exactly 100 times, the chi square would be
zero, or much more uniform than would be expected from random events.   If
we had three faces coming up 110 times each, and three faces coming up 90
times each, the chi square would be 6.000, or just slightly more scattered
than the ideal chi square value of 4.351.   Actually, any chi square value
from 1.15 to 11.1 (with probabilities of 95% to 5%) would be quite
reasonable for the die-throwing experiment.   If the probability is greater
than 99% or less than 1%, it is normally grounds for the rejection of the
result.   Of course, if you run hundreds of trials, you can expect about
one of each 100 chi squares to have a probability greater than 99%, and the
same for a probability less than 1%.
 
     Random events that occur as a continuous distribution instead of as a
small number of discrete categories are not suitable for the use of the chi
square statistic.   In some of these cases, the `Z score' statistic may be
used 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 by
just 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 have
included the degrees of freedom (DF) for the chi square tests.   Note that
all 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.   All
interpreted Basic tests were run for 1,000 batches, and all compiled Basic
tests were run for 10,000 batches.
 
1.  Distribution Test   (DF = 9)
 
     The distribution test (sometimes called the frequency test) checks if
the RNG gives a uniform distribution of numbers.   That is, if all numbers
within its range are equally likely to occur.   This is one of the simplest
tests, and most RNGs pass it.
 
     Random integers (1 through 10) should each occur about 10% of the
time.   The program display of Table 5 shows the batch number, the count
in each of the ten bins for 1000 random integers, the chi square of that
batch, then the cumulative chi square for all batches taken together.   The
heading shows the Ideal (expected) value of each number.
 
     The batch chi square measures the short-term randomness while the
cumulative chi square measures the long-term randomness.   The RNG may pass
one chi square test but not the other.   For example, if the RNG produced a
slight systematic excess of high integers, it could still pass the chi
square test for each batch, but the cumulative chi square would gradually
become larger and larger, showing a greater and greater variation from true
randomness.
 
2.  Mean and Standard Deviation
 
     Random numbers (0 to 1) should average about 0.5000, with a standard
deviation of about 0.2887.   For each batch of 1000 random numbers the
program displays the batch number, the batch mean, the batch standard
deviation, the cumulative mean of all batches, and the cumulative standard
deviation of all batches.
 
3.  Correlation Test
 
     Random numbers (0 to 1) should have a serial correlation coefficient
of approximately zero.   That is, each new number should not depend upon
the previous number in any regular way.   A positive correlation shows a
tendency to run, while a negative correlation shows a tendency to
alternate.
 
     For each batch of 1000 sequential pairs of random numbers, the program
shows the batch number, the correlation coefficient, the batch Z score, and
the 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, the
count of each of the nine possible pairs, the batch chi square, and the
cumulative chi square of all batches taken together.
 
5.  Run Test   (DF = 4)
 
     Random numbers (0 to 1) should have `runs up,' where each succeeding
number 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 next
random number was smaller), one, two, three, and four or more, then the
batch 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 is
less than 0.5 in a predictable manner.   For this test, the number of
random numbers needed to give 1000 `gaps' will also vary a little.
 
     For each batch of 1000 `gaps,' the program shows the batch number, the
count for the batch, the number of `gaps' of length zero (succeeding random
numbers were less than 0.5), one through four, and five or more, then the
batch 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 six
different possible permutations.   That is, LMH, LHM, MLH, MHL, HLM, and
HML are equally likely, where L = low, M = medium, and H = high.
 
     For each batch of 1000 sets of permutations, the program shows the
batch number, the bin count for each permutation, the batch chi square, and
the cumulative chi square for all batches.
 
8.  Mean Square Successive Difference Test
 
     Successive random numbers (0 to 1) should have a predictable average
squared difference.   This routine calculates both eta and the Z score.
Eta is the sum of the squared differences of the successive iterations
divided by the sum of the square of the differences between each iteration
and the mean.   The Z score is then calculated from eta.   A positive Z
score indicates long trends, while a negative Z score indicates short
oscillations.
 
     For each batch of 1000 differences, the program shows the batch
number, 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 predictable
number 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 matches
for 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 the
cumulative chi square.
 
10.  Coupon Collector's Test
 
     For random integers (1 through 10), to get at least one of each
integer takes an average of 29.29 samples.   For each batch of 100 complete
sets, the program shows the batch number, the average number of samples for
the 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-ball
machine), the distribution of balls in the eleven bins at the bottom is
predictable (see Table 4).   Since there will be so few balls in the end
bins, these are combined with the next-to-end bins for the chi square
calculation.
 
     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 cumulative
chi square.
 
12.  Maximum of `T' Test   (DF = 3)
 
     For three successive random integers (1 through 10), each should be
the largest 28.5% of the time, with no maximum (duplication of the high
integer) 14.5% of the time.
 
     For each batch of 1000 sets of three random integers, the program
shows the batch number, the number of `no maximum,' the number of times
that 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 random
walk started at X = 0 and Y = 0, on the average it should take 472 random
steps (right or left, up or down) to hit a wall, with a standard deviation
of 260 steps.   After 25 walks, the average number of steps should be
between 420 and 524.   Also, X and Y, plus and minus walls should be hit
equally often.
 
     For each random walk, the program shows the walk number, the X and Y
values after a wall is hit, the number of steps required for this walk, the
average number of steps for all walks, and (after the first walk) the
standard deviation of the number of steps.   Note that this test is quite
variable.
 
14.  Looping Test
 
     A random number generator, producing floating point random numbers
between zero and one, should not reproduce its original random number, nor
any random number produced later.   If it does, it will loop (produce the
same sequence of numbers again).   This test checks the random numbers for
duplication at user-specified intervals.   Up to 100 random numbers may be
stored and checked against all new random numbers for duplication.   If a
duplication occurs, the sequence numbers of the duplicate random numbers
are displayed.   Otherwise, only the number of random numbers that have
been checked is displayed.
 
     Note that if one duplication is found, many others will follow, as the
same loop of numbers is generated again and again.   To test one million
random 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-DOS
Basic and compiled Basic random number generators?   First, a digression:
 
     While Table 2 shows only six columns for two interpreted Basics and
four compiled Basics, more were tested.   BASICA version 1.19 and GWBasics
versions 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 compiled
Basics, BASCOM and QuickBasic 1 gave identical numeric results, as shown in
the `BASCOM' column.   In addition, the QuickBasics 3, 4, 4.5, plus Basic
Compiler 7 (sometimes called QuickBasic Extended) and Visual Basic for DOS
all 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 in
Table 2, that depends:   If you are using QB64 or the QBasic interpreted
Basic, it passed all tests, at least up to one million numbers (10 million
for QB64), though there are 4 tests for which the results were `questionable'.
 
     However, random number generators in all other Basics failed at least
one test.   The worst was the old BASCOM compiler, which failed three
tests, Run, Gap, and Poker.   The best was the Turbo Basic compiler, which
failed only the Looping test.   If you need `only' 2.5 million or so random
numbers, Turbo Basic would be best.   QuickBasic 4 and its clones failed
two tests and was questionable on two others.   The Microway Basic compiler
failed only one test, and was questionable on one other test.
 
     The two tests that were most difficult for the random number
generators to pass were the Poker test, failed by three RNGs, and the
Pascal's Triangle test, failed by two RNGs.
 
     One other interesting finding was the wide difference in how fast the
RNGs in the various Basic compilers worked.   Table 3 shows that Quick
Basic 4 was far faster than seven of the other compilers, over ten times as
fast as QuickBasic 4.5 or Turbo Basic.   The Microway compiler was nearly
as fast as QuickBasic 4.
 
     Since most of the MS-DOS Basics were from Microsoft, the variety of
results obtained is quite surprising.   I would have expected them to find
a good random number generator and to then stick with it.   But no, after
getting a pretty good RNG for QBasic, they had a couple of very poor ones
for their BASCOM and QuickBasic 1 compilers.   Then for QuickBasic 3 and
later, they apparently use the same, not too good RNG routine.   Odd, to
say the least.
 
     I have heard that more computer time has been spent in testing RNGs
than in actually using them.   After evaluating these RNGs, I can believe
it!
 
Other Systems:
 
     The `RNG' program is correct for running on most Basic interpreters or
for compiling on the most compilers.   One minor change was needed for it
run on the old BASCOM compiler:   Delete all CLS statements.   The BASCOM
compiler does not support CLS.   Also another minor change was used for
running the program with the Microway compiler:   Add a statement to retain
certain variables in the numeric coprocessor's registers.   This speeds up
operation by about 12.5%.
 
If Your RND Fails the Test:
 
     If the RND function in your BASIC (or compiled BASIC) is not good
enough for your requirements, must you give up?   No, you can still
simulate the RND function with a few BASIC statements.   The best form
of RNG is called a `mixed congruential' RNG.   To implement this form of
RNG, you simply select a seed number `R' (which may be zero), multiply
it by a selected number, add a different selected number, then save only
the decimal portion of the resulting number for both the random number
and for the seed number for the next random number.   That sounds rather
complicated, but it's really quite simple:   Just use something like
P = 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 greater
the precision of calculation, the better.   You might try P = 9821,
Q = 0.211327 and R = 0;  I tried these in the subroutine above in a
compiled BASIC version of the RNG test program, and after one million
numbers, it passed all but the Maximum of `T' test.   Of course, I used
double precision for P, Q, and R to give enough significant digits in
the subroutine calculations.   Experiment, or check the references.   I
can't guarantee that these values will work for you, as each BASIC does
its calculations a little differently.   Good luck!
 
Conclusions:
 
     The subject of random numbers may seem simple, but it actually gets
heavily into statistics.   By the use of these RNG tests, you may check
whether your BASIC's RND function is good enough for you to base serious
programs on it, or whether it can be used only for selecting a starting
point for a game.   Even if your RND fails, you may be able to replace
it 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 CLS
20 PRINT"-  -  -  -  -  -  Random Number Generator (RNG) Tests  -  -  -  -  -  -
30 PRINT
40  '              (c)  by Delmer D. Hinrichs  1997
50  '                      2116 S. E. 377th Ave.
60  '                      Washougal, Wash. 98671
70  '
80  '   Written in MS-DOS GWBASIC for IBM PC/XT/AT/PS2 and compatibles
90  '
100 DEFINT I-N :S=.25 :H=.5 :T=.75 :B=0
110 IF II=0 THEN DIM A(100), C(10), J(10), P(10) :II=1
120 FOR I=0 TO 10 :C(I)=0 :J(I)=0 :P(I)=0 :NEXT I
130 FOR I=0 TO 100 :A(I)=0 :NEXT I
140 A=RND(-1)
150 PRINT"Select the test that you want from the following list" :PRINT
160 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" :PRINT
310 PRINT TAB(15); :INPUT"Which Test No."; A$ :A=VAL(A$) :IF A$="0" THEN END
320 ON A GOTO 350,590,800,1040,1290,1570,1840,2170,2410,2750,2970,3290,3580,3810
330 PRINT"Your Entry `"; A$; "' is Illegal.  Try Again"; :GOTO 310
340  '
350 CLS
360 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$="" :WEND
400 FOR J=1 TO 1000
410   IR=INT(RND*10)+1 :J(IR)=J(IR)+1
420   NEXT J
430 X=0 :C=0 :B=B+1
440 FOR I=1 TO 10
450   X=X+(J(I)-100)^2/100
460   C(I)=C(I)+J(I)
470   C=C+(C(I)-100*B)^2/(100*B)
480   NEXT I
490 PRINT USING"#####  ###  ###  ###  ###  ###  ###  ###  ###  ###  ###   ###.## ###.##"; B; J(1); J(2); J(3); J(4); J(5); J(6); J(7); J(8); J(9); J(10); X; C
500 FOR I=1 TO 10 :J(I)=0 :NEXT I
510 IF INKEY$<>"" THEN PRINT :GOTO 540
520 IF B<1000 GOTO 570
530   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 10
570 GOTO 400
580  '
590 CLS
600 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=0
640 R1=0 :R2=0
650 FOR I=1 TO 1000
660   R=RND :R1=R1+R :R2=R2+R*R
670   NEXT I
680 B=B+1 :T=1000*B :C=C+R1 :D=D+R2
690 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, CD
720 IF INKEY$<>"" THEN PRINT :GOTO 750
730 IF B<1000 GOTO 780
740   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 10
780 GOTO 640
790  '
800 CLS
810 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=0
850 X1=0 :X2=0 :Y1=0 :Y2=0 :XY=0 :X=RND
860 FOR I=1 TO 1000
870   Y=X :X=RND
880   X1=X1+X :X2=X2+X*X :Y1=Y1+Y :Y2=Y2+Y*Y :XY=XY+X*Y
890   NEXT I
900 B=B+1 :C1=C1+X1 :C2=C2+X2 :D1=D1+Y1 :D2=D2+Y2 :CD=CD+XY
910 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, CZ
960 IF INKEY$<>"" THEN PRINT :GOTO 990
970 IF B<1000 GOTO 1020
980   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 10
1020 GOTO 850
1030  '
1040 CLS
1050 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$="" :WEND
1090 FOR I=1 TO 1000
1100   K=3*INT(3*RND) + INT(3*RND)
1110   J(K)=J(K)+1
1120   NEXT I
1130 B=B+1 :X=0 :C=0
1140 FOR I=0 TO 8
1150   X=X+(J(I)-111)^2/111
1160   C(I)=C(I)+J(I)
1170   C=C+(C(I)-111*B)^2/(111*B)
1180   NEXT I
1190 PRINT USING"#####    ###  ###  ###  ###  ###  ###  ###  ###  ###   ###.## ###.##"; B, J(0), J(1), J(2), J(3), J(4), J(5), J(6), J(7), J(8), X, C
1200 FOR I=0 TO 8 :J(I)=0 :NEXT I
1210 IF INKEY$<>"" THEN PRINT :GOTO 1240
1220 IF B<1000 GOTO 1270
1230   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 10
1270 GOTO 1090
1280  '
1290 CLS
1300 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.333
1340 WHILE INKEY$="" :WEND :M=0
1350 FOR I=1 TO 1000
1360   D=RND :K=0
1370   R=RND :IF R>D THEN D=R :K=K+1 :GOTO 1370
1380   IF K>4 THEN K=4
1390   J(K)=J(K)+1 :M=M+K+1
1400   NEXT I
1410 B=B+1 :X=0 :C=0
1420 FOR I=0 TO 4
1430   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 I
1470 PRINT USING"#####    ####        ###  ###  ###  ##  ##      ##.##    ##.##"; B, M+1000, J(0), J(1), J(2), J(3), J(4), X, C
1480 FOR I=0 TO 4 :J(I)=0 :NEXT I
1490 IF INKEY$<>"" THEN PRINT :GOTO 1520
1500 IF B<1000 GOTO 1550
1510   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 10
1550 M=0 :GOTO 1350
1560  '
1570 CLS
1580 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=0
1630 FOR I=1 TO 1000
1640   K=K+1 :IF RND>H THEN L=L+1 :GOTO 1640
1650   IF L>5 THEN L=5
1660   J(L)=J(L)+1 :L=0
1670   NEXT I
1680 B=B+1 :X=0 :C=0
1690 FOR I=0 TO 5
1700   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 I
1740 PRINT USING"#####   ####      ###  ###  ### ### ### ###     ##.##    ##.##"; B, K, J(0), J(1), J(2), J(3), J(4), J(5), X, C
1750 FOR I=0 TO 5 :J(I)=0 :NEXT I
1760 IF INKEY$<>"" THEN PRINT :GOTO 1790
1770 IF B<1000 GOTO 1820
1780   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 10
1820 K=0 :GOTO 1630
1830  '
1840 CLS
1850 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=0
1890 WHILE INKEY$="" :WEND
1900 FOR I=1 TO 1000
1910   D=RND :E=RND :F=RND
1920   IF D>E GOTO 1950
1930   IF F>D GOTO 1980
1940   K=2 :GOTO 1990
1950   IF F>D THEN K=3 :GOTO 1990
1960   IF E>F THEN K=0 ELSE K=1
1970   GOTO 1990
1980   IF E>F THEN K=4 :ELSE K=5
1990   J(K)=J(K)+1
2000   NEXT I
2010 B=B+1 :X=0 :C=0
2020 FOR I=0 TO 5
2030   X=X+(J(I)-P)^2/P
2040   C(I)=C(I)+J(I)
2050   C=C+(C(I)-B*P)^2/(B*P)
2060   NEXT I
2070 PRINT USING"#####            ###  ###  ###  ###  ###  ###     ##.##    ##.##"; B, J(0), J(1), J(2), J(3), J(4), J(5), X, C
2080 FOR I=0 TO 5 :J(I)=0 :NEXT I :K=0
2090 IF INKEY$<>"" THEN PRINT :GOTO 2120
2100 IF B<1000 GOTO 2150
2110   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 10
2150 GOTO 1900
2160  '
2170 CLS
2180 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=0
2220 X=RND :D2=0 :X1=0 :X2=0
2230 FOR I=1 TO 1000
2240   Y=X :X=RND
2250   D2=D2+(Y-X)*(Y-X) :X1=X1+X :X2=X2+X*X
2260   NEXT I
2270 B=B+1 :W=1000*B :CD=CD+D2 :C1=C1+X1 :C2=C2+X2
2280 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, CZ
2330 IF INKEY$<>"" THEN PRINT :GOTO 2360
2340 IF B<1000 GOTO 2390
2350   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 10
2390 GOTO 2220
2400  '
2410 CLS
2420 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)=37
2460 WHILE INKEY$="" :WEND
2470 FOR I=1 TO 1000 :K=0
2480   J1=INT(RND*10)+1 :J2=INT(RND*10)+1 :J3=INT(RND*10)+1 :J4=INT(RND*10)+1
2490   IF J1=J2 THEN K=1
2500   IF J1=J3 THEN K=K+1
2510   IF J1=J4 THEN K=K+1
2520   IF J2=J3 THEN K=K+1
2530   IF J2=J4 THEN K=K+1
2540   IF J3=J4 THEN K=K+1
2550   IF K>4 THEN K=4
2560   J(K)=J(K)+1
2570   NEXT I
2580 B=B+1
2590 J(3)=J(3)+J(4) :X=0 :C=0
2600 FOR I=0 TO 3
2610   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 I
2650 PRINT USING"#####            ###  ### ### ### ##             ##.##   ##.##"; B, J(0), J(1), J(2), J(3), J(4), X, C
2660 FOR I=0 TO 4 :J(I)=0 :NEXT I
2670 IF INKEY$<>"" THEN PRINT :GOTO 2700
2680 IF B<1000 GOTO 2730
2690   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 10
2730 GOTO 2470
2740  '
2750 CLS
2760 PRINT"(10)  *   *   *   *   Coupon Collector's Test for RNGs   *   *   *   *"
2770 PRINT, "Average", "Cum.Avg"
2780 PRINT"Ideal", 29.29, 29.29
2790 WHILE INKEY$="" :WEND :C=0 :M=0
2800 FOR I=1 TO 100
2810   FOR L=1 TO 10
2820     K=INT(RND*10)+1 :M=M+1 :IF J(K)=1 GOTO 2820
2830     J(K)=1
2840     NEXT L
2850   FOR K=1 TO 10 :J(K)=0 :NEXT K
2860   NEXT I
2870 B=B+1 :C=C+M
2880 PRINT USING"#####         ###.##        ###.##"; B, M/100, C/(100*B)
2890 IF INKEY$<>"" THEN PRINT :GOTO 2920
2900 IF B<1000 GOTO 2950
2910   PRINT"1000 Trials"
2920   PRINT, "Average", "Cum.Avg"
2930   PRINT"Ideal", 29.29, 29.29
2940   WHILE INKEY$="" :WEND :GOTO 10
2950 M=0 :GOTO 2800
2960  '
2970 CLS
2980 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)=252
3020 P(6)=210 :P(7)=120 :P(8)=45 :P(9)=11
3030 WHILE INKEY$="" :WEND
3040 FOR I=1 TO 1024 :K=10
3050   FOR L=1 TO 10
3060     R=RND
3070     IF R>H THEN K=K+1 ELSE K=K-1
3080     NEXT L :J(K/2)=J(K/2)+1
3090   NEXT I
3100 B=B+1
3110 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=0
3140 FOR I=1 TO 9
3150   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 I
3190 PRINT USING"#####   ##  ##  ##  ###  ###  ###  ###  ###  ##  ## ##    ##.##   ##.##"; B, J0, J1, J(2), J(3), J(4), J(5), J(6), J(7), J(8), J9, JT, X, C
3200 FOR I=0 TO 10 :J(I)=0 :NEXT I
3210 IF INKEY$<>"" THEN PRINT :GOTO 3240
3220 IF B<1000 GOTO 3270
3230   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 10
3270 GOTO 3040
3280  '
3290 CLS
3300 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)=285
3340 WHILE INKEY$="" :WEND
3350 FOR I=1 TO 1000
3360   J=INT(RND*10)+1 :K=INT(RND*10)+1 :L=INT(RND*10)+1
3370   IF J>K THEN IF J>L THEN J(1)=J(1)+1
3380   IF K>J THEN IF K>L THEN J(2)=J(2)+1
3390   IF L>J THEN IF L>K THEN J(3)=J(3)+1
3400   NEXT I
3410 J(0)=1000-J(1)-J(2)-J(3)
3420 B=B+1 :X=0 :C=0
3430 FOR I=0 TO 3
3440   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 I
3480 PRINT USING"#####            ###  ###  ###  ###              ##.##   ##.##"; B, J(0), J(1), J(2), J(3), X, C
3490 FOR I=0 TO 3 :J(I)=0 :NEXT I
3500 IF INKEY$<>"" THEN PRINT :GOTO 3530
3510 IF B<1000 GOTO 3560
3520   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 10
3560 GOTO 3350
3570  '
3580 CLS
3590 :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=0
3630 D=0 :J=0 :K=0
3640 R=RND :D=D+1 :IF R>H GOTO 3680
3650 IF R>S THEN J=J+1 ELSE J=J-1
3660 IF ABS(J)<20 GOTO 3640
3670 GOTO 3700
3680 IF R>T THEN K=K+1 ELSE K=K-1
3690 IF ABS(K)<20 GOTO 3640
3700 B=B+1 :C=C+D :C2=C2+D*D
3710 IF B>1 THEN SD=SQR((C2-C*C/B)/(B-1))
3720 PRINT USING"#####          ###   ###        ####   ####.##     ###.##"; B, J, K, D, C/B, SD
3730 IF INKEY$<>"" THEN PRINT :GOTO 3760
3740 IF B<1000 GOTO 3790
3750   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 10
3790 GOTO 3630
3800  '
3810 CLS
3820 PRINT"(14)   *   *   *   *   *   *   RNG Looping Test   *   *   *   *   *   *"
3830 R=RND :A(0)=R
3840 PRINT :INPUT"Test `RND' at Intervals of (10000 for 1 Million RNDs) "; W
3850 IF W>0 THEN IF W=INT(W) GOTO 3870
3860 PRINT"Your Entry "; W;" is Illegal.  Try Again" :GOTO 3840
3870 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 W
3900     R=RND
3910     FOR K=0 TO J
3920       IF R=A(K) THEN PRINT :PRINT :PRINT"Looping Test FAILED.  RND #"; P+U;           "Duplicates"; K*W :END
3930       NEXT K
3940     NEXT U
3950   A(I)=R
3960   NEXT I
3970 PRINT W*100 :PRINT :PRINT"Looping Test Passed."
3980 WHILE INKEY$="" :WEND :GOTO 10
3990 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.Chi
Ideal  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
    • qbasicmichael.com
    • Email
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.

That page links to:
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 >>>
    • Pete's Qbasic Site
    • Email
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
Download Q-Basics Code Demo: Q-Basics.zip
Download QB64 BAT, IconAdder and VBS shortcuts: QB64BAT.zip
Download QB64 DLL files in a ZIP: Program64.zip

mcalkins

  • Hero Member
  • *****
  • Posts: 1269
    • qbasicmichael.com
    • Email
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/

  • Print