{
KNUTH Random number generator (subtractive method) :

Definitions and globals required by this library :
CONST   bigeven#        = MAXINT - 1;
TYPE    intarray        = ARRAY[1..55] OF INTEGER;
       byte            = 0..255;
VAR     randarray       : intarray;
       randindex,
       seed            : INTEGER;

Contents of this library :
PROCEDURE initknuth(seed : INTEGER);
FUNCTION  rndknuth  : byte;
FUNCTION  random#i  : INTEGER;
FUNCTION  random#r  : REAL;
FUNCTION  random#sr : REAL;
FUNCTION  random#n  : REAL;
Comment:
       The following routines setup the KNUTH subtractive method
random number generator , as described by :
Knuth,D.E.;The art of computer programming (2nd Edition);
            Vol 2 : Seminumerical Algorithms;pp.170-171
            Reading,Ma.;Addison-Wesley Publ.Co. (1981)
Anyone contemplating the use of random numbers is strongly advised to
read Knuth's analysis of the problems involved.
Many of the random number generators widely used have serious faults.
Many of the commonly recommended algorithms , e.g.some of the linear
congruence method routines ,MAY not work properly when coded in a
higher level language such as PASCAL/Z,and some are totally misbehaved.
The linear congruence technic is dependent on both the proper selection
of parameters (modulus , multiplier , and increment ) AND the method
and range of integer arithmetic of the machine and language used.It is
often desirable therefore to have alternative generators available ,
using differring algorithms .

The Knuth generator has the following features :
       1 : It was designed to be written in a high level language ( the
original language was FORTRAN ) and transportable.
       2 : The numbers generated are in a Fibonacci-like sequence but
corrected for the actual nonrandomness of older Fibonacci method generators.

Note Knuth's extreme caution in accepting any results using Monte Carlo
programs dependent on random number generators not explicitly tested for
that application.Such programs should be tried with at least 2 differrent
random number generators,and all random number generators "will fail in
at least one application."

Note also that you dont have to turn off Pascal/Z's error checking to
use this method.

Four random number generators are included below : the first (randomi)
returns a positive INTEGER in the range 0 .. (MAXINT - 1) , the second
(randomr) returns a positive REAL in the range 0.0 .. 1.0,and the third
(randomsr) returns a signed REAL in the range -1.0 .. +1.0.
(Preferably only one of the above 3 should be used in the one program.)
The fourth generator (randomn) returns a random REAL numberfrom a NORMAL
distribution with mean = zero and variance = 1.This FUNCTION requires
FUNCTION randomr .
All 4 generators require the FUNCTION rndknuth,the PROCEDURE initknuth
and a collection of global definitions as listed above.

Modified for Pascal/Z by G.M.Acland : University of Pennsylvania,1982.
}

FUNCTION rndknuth : byte;
{
comment : fills the array "randarray" with 55 pseudo random INTEGERS in
         the range 0..bigeven#. Knuth originally specified 10^9 for
         bigeven# . For Pascal/Z the best number = MAXINT - 1.
         Requires the following definitions globally :
CONST  bigeven#        = MAXINT - 1;
TYPE   "intarray"      = ARRAY[1..55] OF INTEGER;
       "byte"          = 0..255;
VAR    "randarray"     : "intarray";
          Returns the value 1 ( for reinitializing index to "randarray").
}
VAR
       i,j,k   :  INTEGER;
BEGIN
FOR i := 1 TO 55 DO
 BEGIN
   k := I + 31;
   IF k > 55 THEN k := k - 55;
   j := randarray[i] - randarray[k];
   IF j < 0 THEN j := j + bigeven#;
   randarray[i] := j
 END;
rndknuth := 1;
END; { of : FUNCTION rndknuth }

PROCEDURE initknuth(seed : INTEGER);
{
comment : Initializes the GLOBAL randarray. Has the same requirements as
         rndknuth, which FUNCTION  is called by initknuth, plus the input
         value : "seed" : this may be a zero,one or any other positive
         INTEGER value.A useful technic when you want to use a "random"
         seed is to create an integer from the time of day , if you have
         it available to your computer.
e.g.: procedure initseed;
      BEGIN
       timestring := '  :  :  ';
       seed := 0;
       time(timestring);
       FOR i := 1 TO 8 DO seed := seed + ORD(timestring[i]);
      END;
        Note that the GLOBAL randindex is initialized within this routine.
       This was not initialized in the previous version distributed thru
       PZUG.
}
VAR
       i,ii,j,k : INTEGER;
BEGIN
randindex      := 1;
randarray[55]  := seed;
j              := seed;
k              := 1;
FOR i := 1 TO 54 DO
 BEGIN
  ii := (21 * i) MOD 55;
  randarray[ii] := k;
  k := j - k;
  IF k < 0 THEN k := k + bigeven#;
  j := randarray[ii]
 END;
i := rndknuth;
i := rndknuth;
i := rndknuth;
END; { of : PROCEDURE initknuth }

FUNCTION random#r : REAL;
{
comment : Returns a REAL pseudo random number in the range 0.0 .. 1.0.
         Requires the definitions needed by rndknuth and  initknuth
         plus the following global :
 VAR   randindex       : INTEGER;
}
BEGIN
randindex := randindex + 1;
IF randindex > 55 THEN randindex := rndknuth;
random#r  := randarray[randindex]/bigeven#;
END; { of : FUNCTION random#r }

FUNCTION random#i : INTEGER;
{
comment : Returns an INTEGER pseudo random number in the range 0.. bigeven#.
         Requires the definitions needed by rndknuth and  initknuth
         plus the following global :
 VAR   randindex       : INTEGER;
}
BEGIN
randindex := randindex + 1;
IF randindex > 55 THEN randindex := rndknuth;
random#i  := randarray[randindex];
END; { of : FUNCTION random#i }

FUNCTION random#sr : REAL;
{
comment : Returns a REAL pseudo random number in the range -1.0 .. +1.0
         Requires the definitions needed by rndknuth and  initknuth
         plus the following global :
 VAR   randindex       : INTEGER;
}
BEGIN
randindex := randindex + 1;
IF randindex > 55 THEN randindex := rndknuth;
random#sr  := 2.0 * (0.5 - ( randarray[randindex]/bigeven# ));
END; { of : FUNCTION random#sr }

FUNCTION random#n : REAL;
{
comment : Returns a REAL number that is randomly selected from a normally
distributed population whose mean is zero and variance (and standard dev.)
is 1.0.
}
VAR
       n       : INTEGER;
       total   : REAL;
BEGIN
total := -6.0;
FOR n := 1 TO 12 DO total := total + random#r;
random#n := total;
END; { of : function randomn }