!pr0
Random Numbers for Applesoft...............Bob Sander-Cederlof

The RND function in Applesoft is faulty, and many periodicals have loudly proclaimed its faults.  "Call APPLE", Jan 83, pages 29-34, tells them in "RND is Fatally Flawed", and presents an alternative routine which can be called with the USR function.

First, the flaws:  1) the initialization code fails to preset all five bytes of the seed value (only the first four of five are loaded); 2) the RND code uses a poor algorithm, and depends on "tweaks" to make the numbers more random; 3) the RND code does not properly implement the algorithm it appears to be aiming at.

BAD INITIALIZATION.  The initialization code is at $F150 in the Applesoft ROMs.  This loop moves the CHRGET subroutine down to $B1-C8, and is also supposed to copy the random number seed into $C9-CD.  The last byte does not get copied, due to a bug.  Changing $F151 from $1C to $1D would fix it.  Most of us don't really care about this bug, because we are trying to get random numbers for games and the like, and the more random the better:  not copying the last byte could make the numbers generated a little more random from one run to the next.  However, some applications in simulation programs require REPEATABLE sequences of random numbers, so the effect of model changes can be seen independent of the random number generator.

POOR ALGORITHM.  Most generators use an algorithm which makes the next random number by multiplying the previous one by a constant, and adding another constant.  The result is reduced by dividing by a third constant and saving the remainder as the next random number.  More on this later.  The proper choice of the three constants is critical.  I am not sure whether the Applesoft authors just made poor choices, or whether the bugs mentioned below drove them to tweaking.  Tweaking the generated value is often thought to produce even more random results.  In fact, according to authorities like Donald Knuth, they almost always ruin the generator.  Applesoft tweaks the generated value by reversing the middle two bytes of the 32-bit value.  Guess what: it ruins the generator, assuming it was good to start with.

BUGGY ALGORITHM.  The congruency algorithm described in words above will only work properly when integer arithmetic is used.  Applesoft uses floating point arithmetic.  Further, Applesoft arithmetic routines expect five-byte operands.  For some reason the constants used in RND are only four bytes long each.  It appears that the exponents may have been omitted, in the expectation that integer arithmetic was going to be used.  You can see the code for RND at $EFAE.

If you want to see some non-random features using RND, type in and RUN the following program:

       10 HGR:HCOLOR=3
       20 X=RND(1)*280:Y=RND(1)*160
       30 HPLOT X,Y
       40 GO TO 20
You will see the Hi-Res screen being sprinkled with dots.  After about seven minutes, but long before the screen is full, new dots stop appearing.  RND has looped, and is replotting the same sequence of numbers.  Another test disclosed that the repetition starts at the 37,758th "random" number.

Mathematicians have developed many sophisticated tests for random number generators, but Applesoft fails even these simple ones!  Depending on the starting value, you can get the Applesoft generator in a loop.  You never get anywhere near the theoretically possible 4 billion different values.

The Call APPLE article proposes a new algorithm.  It comes with impressive claims and credentials, but I have not found it to be better than a properly implemented congruential algorithm.  The algorithm multiplies the previous seed by 8192, and takes the remainder after dividing by 67099547.  This is a congruency algorithm:

       X(n+1) = ( a * X(n) + c ) mod m

       with a=8192, c=0, m=67099547

I re-implemented the Call APPLE algorithm, and my listing follows.  The Call APPLE version would not quite fit in page 3, but mine does with a little room to spare.  I also dug into some other references and came up with another algorithm, from Knuth.  It is also a congruency, but with a=314159269, c=907633386, and m=2^32.  This turns out to be easier to compute, and according to Knuth it should be "better".  "Better" is in quotes because it is really hard to pin down what are the most important properties.  Anyway this one should have very good characteristics.

The RND function does three different things, depending on the argument.  You write something like R=RND(X).  If X=0, you get the same number as the previous use of RND produced.  If X<0, the absolute value of X becomes the new seed value.  This allows you to control the sequence when you wish, and also to randomize it somewhat by using a "random" seed.  If X>0, you get the next random number.  The value will always be a positive number less than 1.  If you want to generate a number in a range, you multiply by the width of the range and add the starting value.  For example, to generate a random integer between 1 and 10:

       R = INT( RND(1)*10 ) + 1

The programs I have written build a little on the options available with RND.  They all begin with a little routine which hooks in the USR vector.  After executing this, you can write R=USR(X), in other words substitute USR(X) anywhere you would have used RND(X).  But I have added, following the Call APPLE article, the option to automatically generate integers in a range based at 0.  If 0<X<2, you will get the next random fraction.  If X is 2 or greater than 2, you will get a random integer between 0 and X-1.  Thus you can make a random integer between 1 and 10 like this:

       R = USR(10) + 1

as well as with:

       R = INT (USR(1)*10) + 1

I wrote a third program which makes a 16-bit random value.  This one uses the seed at $4E and $4F which the Apple increments continuously whenever the standard monitor input loop is waiting for an input keystroke.  Integer BASIC uses this seed, and as a result is quite valuable in writing games.  My new program gives you all the options stated above, and is significantly quicker than any of the others.  It uses a=19125, c=13843, and m=2^16 in a standard congruency algorithm.

If you are seriously interested in random numbers, you need to read and study Donald Knuth.  Volume 2 of his series "The Art of Computer Programming" is called "Seminumerical Algorithms".  Chapter 3, pages 1-160, is all about random numbers.  (There is only one other chapter in this volume, all about arithmetic in nearly 300 pages!)  Knuth started the series back in the 60's, with the goal of seven volumes covering most of what programmers do.  He finished the first three by 1972, went back and revised the first one, and then evidently got sidetracked into typesetting (several books around a typesetting language he calls "Tex").

Speaking of being sidetracked...!

Knuth ends his chapter with a list of four rules for selecting a, c, and m for congruency algorithms.  Let me summarize those rules here:

!lm+5
1.  The number m is conveniently taken as the word size.  In Applesoft, the floating point mantissa is 32 bits; hence, I chose m=2^32.

2.  If m is a power of 2 (and mine is), pick "a" so that "a mod 8 = 5".  This, together with the rules on choosing c below, ensure that all m values will produced before the series repeats.

3.  Pick "a" between m/100 and m-sqrt(m).  The binary digits should NOT have a simple, regular pattern.  Knuth recommends taking some haphazard constant, such as a=3131492621.

4.  "c" should be odd, and preferable near "m*.2113248654".
!lm-5

Now for the program listings.

The first listing is for my rendition of Call APPLE's algorithm.  Lines 1220-1280 link in the USR vector.  Lines 1370-1450 branch according to the value of the argument of the USR function.  If the argument is negative, lines 1550-1620 set up its absolute value as the new seed.  If the argument is zero, the old seed is used without change, lines 1420-1450.  If positive non-zero, lines 1470-1490 set up the argument as the RANGE.

Lines 1640-1690 calculate the new seed, which will be 8192 times the old seed, modulo 67099547.  8192 is 2^13, so we can multiply be 13 left shifts.  After each shift, if the result is bigger than 67099547, we subtract that value and keep the remainder.  The final result will be some number smaller than 67099547.

Lines 1700-1770 save the new seed, and then divide it by 67099547 to get a fraction for the USR function result.  Lines 1780-1860 check the initial argument to see if you wanted a fraction between 0 and 1, or an integer between 0 and arg-1.  If the latter, the fraction is multiplied by the range and reduced to an integer.

The subroutine named MODULO subtracts 67099547 from the seed value if it would leave a positive remainder, and then renormalizes the result into floating point.

Line 2270 defines the initial seed after loading the program to be 1.0.  If you want some other seed, change this line or be sure to seed it with R=USR(-seed) in your Applesoft program.



<<<<listing of S.USRND S-C>>>>>
!np
The second listing is for my 32-bit algorithm based on Knuth's rules.  Again, lines 1210-1270 set up the USR linkage.  Lines 1360-1400 decide what kind of argument has been used.  If negative, lines 1470-1590 prepare a new seed value.  If zero, the previous value is re-used.  If positive, the argument is the range.

In this version the seed is maintained as a 32-bit integer.  Lines 1470-1590 convert from the floating point form of the argument in FAC to the integer form in SEED.  If the argument happens to be bigger than 2^32, I simply force the exponent to 2^32.

Lines 1600-1690 form the next seed by multiplying by 314159269 and adding 907633386.  The calculation is done in a somewhat tricky way.  Essentially it involves loading 907633386 into the product register, and then adding the partial products of 314159269*seed to that register.  The tricks allow me to do all that with a minimum of program and variable space, and I hope with plenty of speed.  I understood it all this morning, but it is starting to get hazy now.  If you really need a detailed explanation, call me some day.  The modulo 2^32 part is automatic, because bits beyond 32 are thrown away.

Lines 1700-1780 load the seed value into FAC and convert it to a floating point fraction.

Lines 1790-1870 check the range requested.  If less than 2, the fraction is returned as the USR result.  If 2 or more, the fraction is multiplied by the range and integerized.



<<<S.RANDOM KNUTH here>>>
!np
The third listing is cut down from the second one, to produce a 16-bit random number.  The code is very similar to the program above, so I will not describe it line-by-line.  If you want an optimized version of this, the multiply especially could be shortened.


<<<S.RANDOM KEYIN>>>
!np
What do you do if you want even more randomness than you can get from one generator?  You can use two together.  The best way (for greatest randomness) is to use one to select values from a table produced by the other.  First generate, say 50 or 100, random values with one generator.  The generate a random value with the second generator and use it to pick one of the 50 or 100 values.  That picked value is the first number to use.  Then replace the picked value with a new value from the first generator.  Pick another value randomly using the second generator, and so on.  This is analogous to two people working together.  The first person picks a bowlful at random from the universe.  The second person picks items one at a time from the bowl.  The first person keeps randomly picking from the universe to replace the items removed from the bowl by the second person.

You could use the 16-bit generator to pick values from a "bowl" kept full by my 32-bit generator.

Now back to those tests mentioned at the beginning.  I am happy to report that all three of the algorithms listed above completely fill the hi-res screen, no holes left, eventually.

By the way, the August 1981 AAL contained an article about the Integer BASIC RND function, and how to use it from assembly language.
