Random numbers

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.

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

Col 1

Col 2

Col 3

Col 4

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...

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.

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'
- Next one please!
I print the array

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.

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.

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....

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

Page equipped with GoogleBuster technology