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.

Arandom numberin 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 calledpseudo-randomnumbers.

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.

DEFINITIONMODULEGenerator; (* by R. Sutcliffe revised 1993 04 06 *)PROCEDURERandom ( ) :REAL; (* Returns a random number 0 .. 1 *)ENDGenerator.IMPLEMENTATIONMODULEGenerator; (* by R. Sutcliffe revised 1993 04 06 *)FROMRealMathIMPORTexp, ln, pi;VARseed :REAL; (* global variable maintained between calls to Random *)PROCEDURERandom ( ) :REAL; (* Returns a random number 0 .. 1 *)VARtemp :REAL;BEGINtemp := 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 *)RETURNseed;ENDRandom;BEGIN(* initialize seed in body of module*) seed := 4.0;ENDGenerator.

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:

MODULETestGenerator; (* by R. Sutcliffe revised 1993 04 06 *)IMPORTSTextIO;IMPORTSRealIO;IMPORTGenerator;IMPORTRedirStdIO; (* non-standard *)VARcount :CARDINAL;BEGINRedirStdIO.OpenOutput;FORcount := 0TO19DOIFcountMOD4 = 0 (* write them four per line *)THENSTextIO.WriteLn;END; SRealIO.WriteFixed (Generator.Random (), 10,15);END; (* if *) RedirStdIO.CloseOutput;ENDTestGenerator.

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:

PROCEDURERandomize (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)MODm

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 * r_{1} + r_{0}

a = 100 * 1 + 0

and express the multiplication as

r * a = (100 * r_{1} + r_{0}) * (100 * 1 + 0)

= (10 000 * r_{1} * 1) + (100 * (r_{1} * 0 + r_{0} * 1)) + (r_{0} * 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. (9999^{2} 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 = r_{1} * 1 (understood to be times 10000)

2 = r_{1} * 0 + r_{0} * 1 (times 100)

3 = r_{0} * 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:

PROCEDUREMultiply (mplier, mcand :CARDINAL) :CARDINAL; (* Pre: mplier and mcand are both in the range [0 .. 9999] Post: Returns (mplier * mcand) mod 10000 *)TYPEparts =ARRAY[0 .. 1]OFCARDINAL;VARx, y : parts;BEGINx[1] := mplierDIV100; (* break numbers into two parts *) x[0] := mplierMOD100; y[1] := mcandDIV100; y[0] := mcandMOD100;RETURN(((x[1] * y[0] + x[0] * y[1])MOD100) * 100 + x[0] * y[0])MOD10000;ENDMultiply;

With that, it is possible to re-write the procedure Random for cardinals.

PROCEDURERandom ( ) :CARDINAL;CONSTa = 73421;BEGINseed := (Multiply (seed, a) + 1)MOD10000;RETURNseed;ENDRandom;

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:

BEGINWriteString (" What is the random number seed"); ReadCard (r); WriteLn;ENDRandomizer;

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:

PROCEDURERandomizeK ();VARch :CHAR;BEGINSTextIO.WriteString (" Randomizing the seed number "); STextIO.WriteLn; STextIO.WriteString (" Please press a key ");REPEATBusyRead (ch); seed := Random ();UNTILch #CHR(0)ENDRandomizeK;

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 n^{2} = 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:

PROCEDURERealRandom () :REAL;CONSTm = 10 000; (* easy to change this way *)BEGINseed := (Multiply (seed, a) + 1)MODm;RETURNFLOAT(seed) /FLOAT(m)ENDRealRandom;

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.

DEFINITIONMODULERandoms; (* Written by Mark Harder revised by R. Sutcliffe 1993 04 15 *) (* Rnd and Random use the linear congruential generator method *)PROCEDURERandomize ( ); (* Pre: none Post: Sets the seed for random Rnd based on tick count. *)PROCEDURESetSeed (seed :CARDINAL); (* Pre: none Post: Sets the seed to specified seed value. *)PROCEDURERnd ( ) :CARDINAL; (* Pre: none Post: Returns a pseudorandom number of the range of cardinal *)PROCEDURERndInRange (low, high :CARDINAL) :CARDINAL; (* Pre: none Post: Returns a pseudorandom cardinal number in the range low .. high *)PROCEDURERandom ( ) :REAL; (* Pre: none Post: Returns a pseudorandom real number in the range of [0 .. 1) *)ENDRandoms.IMPLEMENTATIONMODULERandoms; (* Written by Mark Harder for the Macintosh computer (system dependent) revised by R. Sutcliffe 1993 04 15 *)FROMSYSTEMIMPORTCAST; (* means this is non-portable *)FROMEventsIMPORT(* Mac system only *) TickCount;FROMRealMathIMPORTln, sqrt, pi;CONSTa = 94621; (* fits pattern of digits and is prime *) maxcard =MAX(CARDINAL); cardmodulus =FLOAT(maxcard) + 1.0; (* can't express cardinal modulus in cardinal *)VARgSeed :CARDINAL;PROCEDUREMultiply (mplier, mcand :CARDINAL) :CARDINAL; (* Returns (x * y) mod maxcard *)TYPEparts =ARRAY[0 .. 1]OFCARDINAL;VARx, y : parts; rtmxcard :CARDINAL; temp :REAL;BEGINrtmxcard :=TRUNC(sqrt (cardmodulus) + 0.01); (* but can its square root *) x[1] := mplierDIVrtmxcard;(* break numbers into two parts *) x[0] := mplierMODrtmxcard; y[1] := mcandDIVrtmxcard; y[0] := mcandMODrtmxcard; temp :=FLOAT(((x[1] * y[0] + x[0] * y[1])MODrtmxcard) * rtmxcard + x[0] * y[0]);IFtemp > cardmodulusTHENtemp := temp - cardmodulus; (* will only have to do once *)END;RETURNTRUNC(temp);ENDMultiply;PROCEDURERandomize (); (* You have to know this cast is appropriate. *)BEGINgSeed := CAST (CARDINAL, TickCount());ENDRandomize;PROCEDURESetSeed (seed :CARDINAL);BEGINgSeed := seed;ENDSetSeed;PROCEDURERnd () :CARDINAL; (* pseudorandom numbers in the range of cardinal *)BEGIN(* linear congruential generator. *) gSeed := Multiply (a, gSeed);IFgSeed < maxcard (* now do the add 1 and mod maxcard part *)THENINC(gSeed);ELSEgSeed := 0;END;RETURN(gSeed)ENDRnd;PROCEDURERndInRange (low, high :CARDINAL) :CARDINAL;VARrange, result :CARDINAL;BEGINrange := high - low; result := (Rnd ()MOD(range + 1) ); (* in range 0 .. range *)INC(result, low); (* adjust range *)RETURN(result);ENDRndInRange;PROCEDURERandom () :REAL;VARtemp :REAL;BEGINRETURNFLOAT(Rnd ())/FLOAT(cardmodulus);ENDRandom;BEGINRandomize ();ENDRandoms.

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.