A random number generator

Today I neede a random number generator that needs to be random and scalable. And whatever I had used before wasn't good enough anymore. So I took my trusty old 'Programming in Oberon' by Wirth and Reiser (Martin, not Hans) and converted the one listed on page 12. I converted it to CARDINAL and INTEGER. Here's the definition module:

```DEFINITION MODULE Random;

PROCEDURE  nr (base  : CARDINAL) : CARDINAL;

PROCEDURE  initSeed (seed	: CARDINAL);

END Random.
```
Could have been bigger and harder to understand I guess. With 'initSeed' you can enter a new seed for the random number generator. With 'nr' you get a new random number between '0' and 'base'. Here's how we donstructed the lot:
```IMPLEMENTATION MODULE Random;

IMPORT	InOut, SysLib;

CONST	a	= 16807;
m	= 2147483647;	(*   2^31 -1	*)
q	= 127773;	(*   m DIV a	*)
r	= 2836;		(*   m MOD a	*)

VAR	z, t	: INTEGER;

PROCEDURE  nr (base  : CARDINAL) : CARDINAL;

VAR	gamma			: INTEGER;
result, factor		: CARDINAL;

BEGIN
gamma := a * (z MOD q) - r * (z DIV q);
factor := m DIV base;
IF  gamma > 0  THEN  z := gamma  ELSE  z := gamma + m  END;
result := CARDINAL (z) DIV factor;
RETURN result
END nr;

PROCEDURE  initSeed (seed	: CARDINAL);

BEGIN
z := seed
END initSeed;

BEGIN
SysLib.time (t);
z := CARDINAL (t)
END Random.
```
When the library is initialised, the system time is used to set the seed for the generator. From then on, new values are calculated on the fly.

The 'factor' is used for scaling. I first had a simpler solution:

```   result := CARDINAL (z) MOD base;
```
but that produced random numbers that were relatively close together. Now, with the factor, things look much better.

Running the random number generator

Below is a short test program to check the random number generator:

```MODULE tln;

IMPORT  InOut, Random;

VAR	n	: CARDINAL;

BEGIN
FOR  n := 1  TO  100  DO
InOut.WriteCard (Random.nr (1000), 10);
InOut.WriteLn
END
END tln.
```
It produced the following list:

Col 1 Col 2 Col 3 Col 4
```  232
510
975
513
736
66
991
457
216
83
191
307
626
665
90
686
253
875
246
773
933
718
755
387
995
```
```  976
173
13
803
548
608
190
888
224
161
961
158
655
115
607
21
834
994
581
537
190
218
637
766
956
```
```  268
712
146
819
132
104
826
98
80
230
178
101
44
548
538
921
513
861
966
478
120
843
371
620
716
```
```   88
95
596
192
367
443
89
527
301
294
838
625
516
119
620
607
63
925
101
772
716
866
229
118
463
```

Don't show the accompanying graph of random data to an environmentalist. He'll see proof of global warming in it... Still, it's the output of the random number generator...

Checking the generator

Nice the table is, it doesn't really show the quality of the random numbers generated. The line graph is excellent proof for the environmentallists, but it's not for us, engineers. The sample was just too small. So I adapted the tln program as follows:

```
MODULE tln;

IMPORT  InOut, Random;

VAR	n	: CARDINAL;

BEGIN
FOR  n := 1  TO  10000  DO
InOut.WriteCard (Random.nr (1024), 10);
InOut.WriteLn
END
END tln.
```
It now generates 10000 numbers in the range between 0 and 1023 (inclusive). And I invented a simple yet effective way to see in a glance if the generator is sound or not. Seeing is believing. Period. To spoil your fun... I won't say it was good or not. But here's the source for the visual checker program:
```
MODULE chkrnd;

IMPORT	InOut;

VAR	field		: ARRAY [0..15] OF ARRAY [0..63] OF CHAR;
dot, n, m	: CARDINAL;

BEGIN
FOR  n := 0  TO  15  DO
FOR  m := 0  TO  63  DO
field [n, m] := ' '
END
END;
LOOP
IF  InOut.Done () = FALSE  THEN  EXIT  END;
m := dot MOD 64;
n := dot DIV 64;
field [n, m] := 'X'
END;
InOut.WriteLn;
FOR  n := 0  TO  15  DO  InOut.WriteString (field [n]); InOut.WriteLn  END;
InOut.WriteBf
END chkrnd.
```
• I create an array of 16 lines of 64 tokens
• The array is cleared to all spaces
• I enter a LOOP
• I read a random number
• If nothing was read, we're done
• I calculate the horizontal position
• I calculate the line number
• I fill the table entry witha big 'X'
• I print the array
Simple yet effective, methinks.

Sample size Distribution

# 10000

This doesn't mean the generator is true. For that, we need to make an average and standard deviation calculation (which is left for tomorrow). Still, it shows that (given a fair sample size) all values are used at least one.

• After 100 samples, not all 1000 vales have been covered, but that's not more than logical of course.
• After 1000 samples, the majority of numbers was used at least once. Still, a lot of gaps are present, indicating that the sequence IS random.
• After 5000 samples, we use a 400% oversampling and that's clearly visible. Only very few empty spot. All randomly spaced out. A good sign!
• After 10000 samples, all values have been used at least once. Our generator is good, to say the least.
Time to see if Linux has a statistics package from the command line.

Command line statistix

Linux does have a wc but nothing similar that handles statistics. But now it has: Stx. Stx takes its arguments from StdIn and pits out averga and standard deviation on StdOut.

```MODULE Stx;

IMPORT	InOut, MathLib, NumConv, Strings;

VAR	sx, sx2, n	: REAL;
num		: Strings.String;
val		: LONGCARD;
ok		: BOOLEAN;

PROCEDURE  Sum (val	: INTEGER);

VAR	Val		: REAL;

BEGIN
Val := MathLib.real (val);
sx := sx + Val;
sx2 := sx2 + Val * Val
END Sum;

PROCEDURE  Avg;

VAR	avg		: REAL;

BEGIN
avg := sx / n;
InOut.WriteString ("Average = ");
InOut.WriteReal (avg, 8, 2);
InOut.WriteLn
END Avg;

PROCEDURE Std;

VAR	sigma		: REAL;

BEGIN
sigma := MathLib.sqrt ((sx2 - (sx * sx) / n) / (n - 1.0) );
InOut.WriteString ("Sigma = ");
InOut.WriteReal (sigma, 8, 2);
InOut.WriteLn;
END Std;

BEGIN
n := 0.0;
sx := 0.0;
sx2 := 0.0;
LOOP
IF  NOT InOut.Done ()  THEN  EXIT  END;
NumConv.Str2Num (val, 10, num, ok);
IF  NOT ok  THEN
InOut.WriteString ("Not a number! Check ");
InOut.WriteString (num);
InOut.WriteLn;
HALT
END;
Sum (INTEGER (val));
n := n + 1.0
END;
InOut.WriteString ("Number of samples : ");
InOut.WriteReal (n, 8, 0);
InOut.WriteLn;
Avg;
Std;
InOut.WriteString ("Done. Have a nice day.");
InOut.WriteLn;
InOut.WriteBf
END Stx.
```
The first version was written using CARDINALs and INTEGERs where possible, but from sample sizes exceeding 5000 you ran into problems (one of the squares overflowed into a negative number...) so I REALed the full source.
```
jan@Beryllium:~/modula/lib\$ tln >File
jan@Beryllium:~/modula/lib\$ chkrnd <File

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
jan@Beryllium:~/modula/lib\$ Stx <File
Number of samples :    10000
Average =   510.23
Sigma =   296.64
Done. Have a nice day.
jan@Beryllium:~/modula/lib\$
```
Well, that doesn't look too bad, all in all. All 10.000 values are used. The average is 510 (with a maximum of 1023, this is spot on the middle of the population. The standard deviation (for n-1) is 300, which looks a bit high to me, but it shows the amount of chaos and chaos is what we wanted in the first place. Also, this is by no means a gaussian curve, so the standard deviation as such is rather meaningless. What we would need here, is a frequency chart. See how far we get here....

Frequency diagram

So it was time to make a frequency lister, similar to the chkrnd lister:

```MODULE chkFrq;

IMPORT	InOut;

VAR	field		: ARRAY [0..15] OF ARRAY [0..63] OF CHAR;
dot, n, m	: CARDINAL;

BEGIN
FOR  n := 0  TO  15  DO
FOR  m := 0  TO  63  DO
field [n, m] := 'A'
END
END;
LOOP
IF  InOut.Done () = FALSE  THEN  EXIT  END;
m := dot MOD 64;
n := dot DIV 64;
IF  field [n, m] < 'Z'  THEN  INC (field [n, m])  END
END;
InOut.WriteLn;
FOR  n := 0  TO  15  DO  InOut.WriteString (field [n]); InOut.WriteLn  END;
InOut.WriteBf
END chkFrq.
```
The array is now filled with A's and if a value is found, the letter is incremented. After 10000 samples, this is what it looks like:
```
jan@Beryllium:~/modula/lib\$ chkFrq <File

NLOJKKJKJMIILJMLFLKHLKHTHIHDJJMINLILIOHGHJHMJFMMDLGQPJHFJLKKJLKO
PFJKLOLFLLGPKJMNHFLHOMJJMHKKOIGIJKLMOLLIMSHJGRIHJOGJKGGNFJIOUIKP
EIJNLLMJDLIKMLOLFHPGNHIQHKQNLKJKMQFONKHHLHLHFIVKNQJMNKLILFPOMJHF
FTOHONMJFKIPFOJLINTQMKPQEIPQGIQIKJIRJKHQJLNHGENQGIJKKJMKIJKGJJLL
KIMKIFIGKIKKJQNFLGIIFRHMOKJEFIMNHLDDOPVMEGFKOMMLHIJNMEGFKKDHKFLF
KFFINPHKJNIHMJJOHGNIKIILHLIHMHNMJGIJIHGIHJFKNKLICINPOGOMJOMOGLML
KNLHJJFHICKGKIOLIEKLLQJIFMLMJMKKLHKIOMIITKLKKFKJJLMNHMKJLLGIFOIG
MPLGFMFJGIJNNIINKJOOMHGJEJKONHLGHHLFKHGKIHHHNGIKLKKEPNMMLHHGFOMK
HKPJDHMKOOGIJLMHIMIJLJKJLHIIFKPIKILHKINGLKGKQNNHJKOJJGENKEHNJMIP
NLKFJOIHLJJPMKHJDJKKOKIILNMPJPHIIJNIKLLJFKJIIELMHLIMHFJLILJHHLHG
JKGEPMERLQLNJHKJNJKKILFMHOKLMGKKHLLPLFGKKNKNMJKQFNJKJQGFINQKKNML
JJNKKOMOKLHKHKJIOLJNJNNMMLNIJOKLSHNMJGOLPHLLHMKJJGKIMJEPREPKILLI
KPKKKHNILGJMOHJLRILGNKKFOPFJJKHOHHKHHNNJHHHIMHNIMIIKMILGKMIIHGLH
INJILMMGKJJRFILJGFQIHJMILIHGDJIQHMNIHGJFLPDHNNMKJJNFMNGIJQOIIIKG
GPGMIHMJIEOQFFKPMHIOJCKFIILMLJFMNHIKHOLKRKJIHJLLJFILLGIOILJHIKLJ
LKRINLKJNMIFNNJJKJMHLIVMNGILIONMGFJKNNGLFJTNNKHMLGJFGHOLLMMJGIHK
```
One C, 10 D's and the rest is more frequent than 4. The result is graphic but not easily comprehended. Perhaps I need to create a file that must be processed by Kspread. See how far we get here.

Kspread crashes when it needs to make graphs on this amount of data...

Page created 22 June 2009 and