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
      InOut.ReadCard (dot);
      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.
   
Simple yet effective, methinks.

Sample size Distribution

100

1000

5000

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.

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
      InOut.ReadString (num);
      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
      InOut.ReadCard (dot);
      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