A number of applications in mathematics and computer science require collections of random numbers over some range. For instance, one may wish to take a political poll by choosing telephone numbers at random. A list of the possible exchanges could be made, and a household or business in one of them selected by generating a random number between 1 and 9999 for the last four digits of the call.
A random number in some range is chosen in such a way that every number is equally likely to be selected on each occasion.
In most cases, including the one just mentioned, a sequence of random numbers, rather than a single number is required. Not only ought each number to be random, but so should the order in which they occur in the sequence. A great deal of work has been done on the properties of such collections of numbers and it is possible to obtain carefully worked out tables of such numbers that satisfy the most stringent criteria of randomness.
What is desired here is to fashion such a sequence using the computer. One difficulty of doing this is the fact that is that once a given piece of code has been written for a computer, the same sequence of random numbers will be generated each time. The numbers produced may well pass various tests of randomness, but they will not be truly random in at least one important sense--they are predetermined.
Random numbers generated by formula in a pre-determined sequence are called pseudo-random numbers.
Some sequences of pseudo-random numbers are better than others. Suppose one were working over the range [0 ..99] and used the sequence 0, 1, 2, 3, .. 99 or the sequence 0, 2, 4, 6, .. 98, 1, 3, 5, ...99. In both cases, each number occurs only once (is equally represented), but in neither is any given number equally likely to be at any given position. Moreover, if the second did not start over (mod 99), only half the numbers would occur at all. A random number generator must be able to do better than that.
As already noted, it is difficult to generate a sequence of truly random numbers. However, the following function might do for many non-critical applications. It does generate pseudo-random numbers--and purists might argue that it does not do so very well, but it will serve for some purposes.
DEFINITION MODULE Generator;
(* by R. Sutcliffe revised 1993 04 06 *)
PROCEDURE Random ( ) : REAL;
(* Returns a random number 0 .. 1 *)
END Generator.
IMPLEMENTATION MODULE Generator;
(* by R. Sutcliffe revised 1993 04 06 *)
FROM RealMath IMPORT
exp, ln, pi;
VAR
seed : REAL;
(* global variable maintained between calls to Random *)
PROCEDURE Random ( ) : REAL;
(* Returns a random number 0 .. 1 *)
VAR
temp : REAL;
BEGIN
temp := seed + pi; (* scramble up digits of seed *)
temp := exp (5.0 * ln (temp)); (* scramble them some more *)
seed := temp - FLOAT (TRUNC (temp));
(* remove whole number part *)
RETURN seed;
END Random;
BEGIN (* initialize seed in body of module*)
seed := 4.0;
END Generator.
Notice that in true Modula-2 style the value of seed is maintained, hidden away from the main program inside the library module Generator until the next time Random is called. A brief test program follows:
MODULE TestGenerator;
(* by R. Sutcliffe revised 1993 04 06 *)
IMPORT STextIO;
IMPORT SRealIO;
IMPORT Generator;
IMPORT RedirStdIO; (* non-standard *)
VAR
count : CARDINAL;
BEGIN
RedirStdIO.OpenOutput;
FOR count := 0 TO 19
DO
IF count MOD 4 = 0 (* write them four per line *)
THEN
STextIO.WriteLn;
END;
SRealIO.WriteFixed (Generator.Random (), 10,15);
END; (* if *)
RedirStdIO.CloseOutput;
END TestGenerator.
0.0000000000 0.0198669434 0.8190307617 0.5805053711 0.3963012695 0.2736206055 0.6108093262 0.9552612305 0.1240234375 0.3860473633 0.2877807617 0.3228149414 0.0510864258 0.7236328125 0.7287597656 0.4659423828 0.0164794922 0.1307373047 0.2195129395 0.9540710449
These numbers look somewhat random, but the same sequence is obtained each time this program is run as things stand. This library module could be improved with the addition of:
PROCEDURE Randomize (newSeed : REAL); (* reset the random number seed to a user supplied one *)
The most common technique for producing pseudo-random cardinal numbers, that satisfy most tests for randomness is called the linear congruential method. It generates a new random number by multiplying the previous one by a constant a, adding one, and take the remainder modulo (MOD) a second constant m. The result is a new random number that lies between zero and m - 1. That is,
r := (r * a + 1) MOD m
Any initial value for r can be used as the seed to start things off and, as before, this seed will in fact determine the entire sequence.
Now, a close examination of this formula reveals a potentially serious problem, for if one assumes that r and a are of type CARDINAL, then unless both are kept very small, or the range for CARDINAL is large, the calculations are frequently going to overflow the data type, possibly on the very first attempt to obtain r * a.
Suppose one were generating random numbers in the range [0 .. 9999]. To prevent overflows, say when the previous value of r had been close to 9999, the value for a would have to be less than six, but this is not too practical. Indeed, without going into the details here, it has been found that the best results are obtained if a has one digit less than m, ends with the digits "21" where x is even, and otherwise has no particular pattern in its digits. What the calculation needs to do, then, is (r * a) MOD 10 000 but without any possible overflow when (r * a) is computed. To accomplish this, break the number into two parts. Since 10 000 is the square of 100, it is convenient to write
r = 100 * r1 + r0
a = 100 * 1 + 0
and express the multiplication as
r * a = (100 * r1 + r0) * (100 * 1 + 0)
= (10 000 * r1 * 1) + (100 * (r1 * 0 + r0 * 1)) + (r0 * 0)
In this answer, one is now interested only in the third term and the rightmost two significant digits of the second term as the first one is already more than 10 000.
7245 * 6175
= (100 * 72 + 45) + (100 * 61 + 75)
= 10000 (72 * 61) + 100 (72 * 75 + 45 * 61) + 45 * 75
= 10000 (4392) + 100 (8145) + 3375 (* break up the second term *)
= 10000 (4392 + 81) + 4500 + 3375
= 10000 (4473) + 7875
Had the last number also overflowed 10 000, one step would be added to carry the one to the first term. Note, for later use, that this is a way to multiply two numbers in the range 0 .. 9999 and keep track of all the digits in the answer--in this case, store 7875 in one CARDINAL, and 4473 in another. (99992 is too large to fit into the cardinal type in many implementations. Even if not, a similar argument could be applied to 9999999, if that were the largest such number representable in the cardinal type.)
Once the whole process in finished in this situation, however, the focus is on only in the rightmost four digits (7875). Putting these observations together in a problem refinement and writing the code yields:
Write a procedure to multiply two numbers in the range [0 .. 9999] and retain only the last four digits.
Break the numbers into two parts of the form 100 * 1 + 0.
for a number n
1 = n DIV 100
0 = n MOD 100
Multiply these as binomials as shown above.
The resulting terms are:
1 = r1 * 1 (understood to be times 10000)
2 = r1 * 0 + r0 * 1 (times 100)
3 = r0 * 0
Add together the parts affecting the rightmost digits.
answer =100 (2 MOD 100) + 3
Discard any portion over 10000.
answer = answer MOD 10000
Here is the code for this procedure:
PROCEDURE Multiply (mplier, mcand : CARDINAL) : CARDINAL;
(* Pre: mplier and mcand are both in the range [0 .. 9999]
Post: Returns (mplier * mcand) mod 10000 *)
TYPE
parts = ARRAY [0 .. 1] OF CARDINAL;
VAR
x, y : parts;
BEGIN
x[1] := mplier DIV 100; (* break numbers into two parts *)
x[0] := mplier MOD 100;
y[1] := mcand DIV 100;
y[0] := mcand MOD 100;
RETURN
(((x[1] * y[0] + x[0] * y[1]) MOD 100) * 100
+ x[0] * y[0]) MOD 10000;
END Multiply;
With that, it is possible to re-write the procedure Random for cardinals.
PROCEDURE Random ( ) : CARDINAL; CONST a = 73421; BEGIN seed := (Multiply (seed, a) + 1) MOD 10000; RETURN seed; END Random;
Finally, there are several ways to initialize r for the first time in the body of the library module containing this procedure. It could be read in from the keyboard, or an attempt could be made to generate it in some truly random fashion. So the body could look like this:
BEGIN
WriteString (" What is the random number seed");
ReadCard (r);
WriteLn;
END Randomizer;
Some versions of Modula-2 provide a procedure BusyRead. This procedure polls the keyboard for a character and immediately returns either CHR(0) if nothing has yet been typed, or the character typed if one was. So, perhaps randomizing the seed could be done like this:
PROCEDURE RandomizeK ();
VAR
ch : CHAR;
BEGIN
STextIO.WriteString (" Randomizing the seed number ");
STextIO.WriteLn;
STextIO.WriteString (" Please press a key ");
REPEAT
BusyRead (ch);
seed := Random ();
UNTIL ch # CHR (0)
END RandomizeK;
The idea here is that seed begins with whatever number happens to be in the memory at the time the procedure is entered. This depends on the sequence of events since start up of the computer. Random is called repeatedly until the user presses a key. Since everyone has a different reaction time and this loop would execute many times before the keypress, the value of the initial seed and thus of the whole sequence would be more genuinely random.
Still another method is available to those who have an intimate knowledge of the workings of the particular hardware on which their randomizer is running. Many systems change a particular memory location at regular (or irregular) intervals, and if a method of accessing such a location is known, its contents may be interpreted as a cardinal, and used for the seed.
NOTE: Forcible re-interpretation of the data in some memory location, initially named by a variable of some other type is done using newName := CAST (newTypeName, variableName); where CAST is imported from the module SYSTEM. Use of any procedures or types in the module SYSTEM implies that the code is not portable to some other computing system.
A popular choice is a memory location representing the value held for a clock of some kind. The latter may be a real-time clock, whose contents are a representation of the time and (possibly) the date. It could be a pseudo-clock, sometimes called a tick count, which is a location automatically incremented by the operating system as part of various other tasks it performs. The contents of the tick count location, where applicable, are predictable for a given time and sequence of events after system start up, but as both of these may vary considerably before the randomizing function is called, the results may be quite satisfactory.
Note also that all references to 10000 and 100 in the above discussion could be replaced by some other numbers, say m and n with n2 = m as long as m is less than half of the largest possible CARDINAL. A suitable expansion is possible if the nonstandard LONGCARD is available or the type CARDINAL has a sufficiently large maximum value. It might be best to write the algorithm with references to 10 000 and 100 replaced by constants so that these can easily be changed.
It is also worth noting that the last digits of a pseudo-random number sequence cycle through a regular pattern (test this!) and that sequences of random numbers over some specified range that are produced from some generic generator should therefore use the first digits for the numbers in the sequence.
For this reason, most generators return a REAL between 0 and 1 (say, by dividing by 10 000 in this case). Given this assumption, a program would then acquire a random number from such a generator--say in the range [0 .. j - 1] by executing the statement:
num := TRUNC (j * RealRandom () )
The procedure RealRandom might then be:
PROCEDURE RealRandom () : REAL; CONST m = 10 000; (* easy to change this way *) BEGIN seed := (Multiply (seed, a) + 1) MOD m; RETURN FLOAT (seed) / FLOAT (m) END RealRandom;
Here is a library module based on some of the ideas in this section. The implementation was written specifically for the Macintosh computer and uses a system utility to access the memory location containing the tick count, so the code is not portable to other systems.
DEFINITION MODULE Randoms;
(* Written by Mark Harder
revised by R. Sutcliffe 1993 04 15 *)
(* Rnd and Random use the linear congruential generator method *)
PROCEDURE Randomize ( );
(* Pre: none
Post: Sets the seed for random Rnd based on tick count. *)
PROCEDURE SetSeed (seed : CARDINAL);
(* Pre: none
Post: Sets the seed to specified seed value. *)
PROCEDURE Rnd ( ) : CARDINAL;
(* Pre: none
Post: Returns a pseudorandom number of the range of cardinal *)
PROCEDURE RndInRange (low, high : CARDINAL) : CARDINAL;
(* Pre: none
Post: Returns a pseudorandom cardinal number in the range low .. high *)
PROCEDURE Random ( ) : REAL;
(* Pre: none
Post: Returns a pseudorandom real number in the range of [0 .. 1) *)
END Randoms.
IMPLEMENTATION MODULE Randoms;
(* Written by Mark Harder
for the Macintosh computer (system dependent)
revised by R. Sutcliffe 1993 04 15 *)
FROM SYSTEM IMPORT
CAST; (* means this is non-portable *)
FROM Events IMPORT (* Mac system only *)
TickCount;
FROM RealMath IMPORT
ln, sqrt, pi;
CONST
a = 94621; (* fits pattern of digits and is prime *)
maxcard = MAX (CARDINAL);
cardmodulus = FLOAT (maxcard) + 1.0;
(* can't express cardinal modulus in cardinal *)
VAR
gSeed : CARDINAL;
PROCEDURE Multiply (mplier, mcand : CARDINAL) : CARDINAL;
(* Returns (x * y) mod maxcard *)
TYPE
parts = ARRAY [0 .. 1] OF CARDINAL;
VAR
x, y : parts;
rtmxcard : CARDINAL;
temp : REAL;
BEGIN
rtmxcard := TRUNC (sqrt (cardmodulus) + 0.01);
(* but can its square root *)
x[1] := mplier DIV rtmxcard;(* break numbers into two parts *)
x[0] := mplier MOD rtmxcard;
y[1] := mcand DIV rtmxcard;
y[0] := mcand MOD rtmxcard;
temp := FLOAT (((x[1] * y[0] + x[0] * y[1])
MOD rtmxcard) * rtmxcard + x[0] * y[0]);
IF temp > cardmodulus
THEN
temp := temp - cardmodulus; (* will only have to do once *)
END;
RETURN TRUNC (temp);
END Multiply;
PROCEDURE Randomize ();
(* You have to know this cast is appropriate. *)
BEGIN
gSeed := CAST (CARDINAL, TickCount());
END Randomize;
PROCEDURE SetSeed (seed : CARDINAL);
BEGIN
gSeed := seed;
END SetSeed;
PROCEDURE Rnd () : CARDINAL;
(* pseudorandom numbers in the range of cardinal *)
BEGIN
(* linear congruential generator. *)
gSeed := Multiply (a, gSeed);
IF gSeed < maxcard (* now do the add 1 and mod maxcard part *)
THEN
INC (gSeed);
ELSE
gSeed := 0;
END;
RETURN (gSeed)
END Rnd;
PROCEDURE RndInRange (low, high : CARDINAL) : CARDINAL;
VAR
range, result : CARDINAL;
BEGIN
range := high - low;
result := (Rnd () MOD (range + 1) ); (* in range 0 .. range *)
INC (result, low); (* adjust range *)
RETURN (result);
END RndInRange;
PROCEDURE Random () : REAL;
VAR
temp : REAL;
BEGIN
RETURN FLOAT (Rnd ())/ FLOAT (cardmodulus);
END Random;
BEGIN
Randomize ();
END Randoms.
Note the changes to the multiplication procedure to make it operate over the entire range of cardinals. There are other ways of doing this, including converting both numbers to reals, performing a real multiplication, stripping the amount over maxcard and converting back to a cardinal. You will be asked to produce this and some other refinements of this material in the exercises at the end of this chapter. With the proper imports, a client program could itself employ the keyboard response method of randomizing the seed by changing the line in the procedure RandomizeK earlier in this section from
seed := Random ();
to
Randoms.Randomize ();
It is also worth observing at this point that some versions of Modula-2 provide random number generator functions provided in a library module. As these are not in any sense standard items, the names, syntax, and location of such items will vary from one implementation to another and they are likely to be machine dependent. Check the documentation for details.