home *** CD-ROM | disk | FTP | other *** search
/ Oakland CPM Archive / oakcpm.iso / sigm / vol131 / rndknuth.lib < prev    next >
Encoding:
Text File  |  1984-04-29  |  6.0 KB  |  182 lines
  1. {
  2. KNUTH Random number generator (subtractive method) :
  3.  
  4. Definitions and globals required by this library :
  5. CONST    bigeven#    = MAXINT - 1;
  6. TYPE    intarray    = ARRAY[1..55] OF INTEGER;
  7.     byte        = 0..255;
  8. VAR    randarray    : intarray;
  9.       randindex,
  10.     seed        : INTEGER;
  11.  
  12. Contents of this library :
  13. PROCEDURE initknuth(VAR randarray : intarray;seed : INTEGER);
  14. FUNCTION  rndknuth(VAR randarray : intarray) : byte;
  15. FUNCTION  random#i  : INTEGER;
  16. FUNCTION  random#r  : REAL;
  17. FUNCTION  random#sr : REAL;
  18. FUNCTION  random#n  : REAL;
  19. Comment:
  20.     The following routines setup the KNUTH subtractive method
  21. random number generator , as described by :
  22.  Knuth,D.E.;The art of computer programming (2nd Edition);
  23.          Vol 2 : Seminumerical Algorithms;pp.170-171
  24.          Reading,Ma.;Addison-Wesley Publ.Co. (1981)
  25. Anyone contemplating the use of random numbers is strongly advised to
  26. read Knuth's analysis of the problems involved.
  27. Many of the random number generators widely used have serious faults.
  28. Many of the commonly recommended algorithms , e.g.some of the linear
  29. congruence method routines ,MAY not work properly when coded in a 
  30. higher level language such as PASCAL/Z,and some are totally misbehaved.
  31. The linear congruence technic is dependent on both the proper selection 
  32. of parameters (modulus , multiplier , and increment ) AND the method
  33. and range of integer arithmetic of the machine and language used.It is
  34. often desirable therefore to have alternative generators available , 
  35. using differring algorithms .
  36.  
  37. The Knuth generator has the following features :
  38.     1 : It was designed to be written in a high level language ( the
  39. original language was FORTRAN ) and transportable.
  40.     2 : The numbers generated are in a Fibonacci-like sequence but
  41. corrected for the actual nonrandomness of older Fibonacci method generators.
  42.  
  43. Note Knuth's extreme caution in accepting any results using Monte Carlo
  44. programs dependent on random number generators not explicitly tested for
  45. that application.Such programs should be tried with at least 2 differrent
  46. random number generators,and all random number generators "will fail in
  47. at least one application."
  48.  
  49. Note also that you dont have to turn off Pascal/Z's error checking to
  50. use this method.
  51.  
  52. Three random number generators are included below : the first (randomi)
  53. returns a positive INTEGER in the range 0 .. (MAXINT - 1) , the second
  54. (randomr) returns a positive REAL in the range 0.0 .. 1.0,and the third
  55. (randomsr) returns a signed REAL in the range -1.0 .. +1.0.
  56. Preferably only one of these should be used in the one program.
  57.  
  58. A fourth generator (randomn) returns a random REAL numberfrom a NORMAL 
  59. distribution with mean = zero and variance = 1.This FUNCTION requires
  60. FUNCTION randomr .
  61.  
  62. All 4 generators require the FUNCTION rndknuth,the PROCEDURE initknuth
  63. and a collection of global definitions as listed above.
  64.  
  65. Modified for Pascal/Z by G.M.Acland : University of Pennsylvania,1982.
  66. }
  67.  
  68. FUNCTION rndknuth(VAR randarray : intarray) : byte;
  69. {
  70. comment : fills the array "randarray" with 55 pseudo random INTEGERS in
  71.       the range 0..bigeven#. Knuth originally specified 10^9 for
  72.       bigeven# . For Pascal/Z the best number = MAXINT - 1.
  73.       Requires the following definitions globally :
  74.  CONST    bigeven#    = MAXINT - 1;
  75.  TYPE    "intarray"    = ARRAY[1..55] OF INTEGER;
  76.     "byte"        = 0..255;
  77.  VAR    "randarray"    : "intarray";
  78.        Returns the value 1 ( for reinitializing index to "randarray").
  79. }
  80. VAR
  81.     i,j,k    :  INTEGER;
  82. BEGIN
  83.  FOR i := 1 TO 55 DO
  84.   BEGIN
  85.     k := I + 31;
  86.     IF k > 55 THEN k := k - 55;
  87.     j := randarray[i] - randarray[k];
  88.     IF j < 0 THEN j := j + bigeven#;
  89.     randarray[i] := j
  90.   END;
  91.  rndknuth := 1;
  92. END; { of : FUNCTION rndknuth }
  93.  
  94. PROCEDURE initknuth(VAR randarray : intarray;seed : INTEGER);
  95. {
  96. comment : Initializes randarray.Has the same requirements as rndknuth ,
  97.       which FUNCTION  is called by initknuth,plus the input value :
  98.       "seed" : this may be a zero,one or any other positive
  99.       INTEGER value.A useful technic when you want to use a "random"
  100.       seed is to create an integer from the time of day , if you have
  101.       it available to your computer.
  102.  e.g.: procedure initseed;
  103.        BEGIN
  104.         timestring := '  :  :  ';
  105.         seed := 0;
  106.         time(timestring);
  107.         FOR i := 1 TO 8 DO seed := seed + ORD(timestring[i]);
  108.        END;
  109. }
  110. VAR
  111.     i,ii,j,k : INTEGER;
  112. BEGIN
  113.  randarray[55]    := seed;
  114.  j        := seed;
  115.  k        := 1;
  116.  FOR i := 1 TO 54 DO
  117.   BEGIN
  118.    ii := (21 * i) MOD 55;
  119.    randarray[ii] := k;
  120.    k := j - k;
  121.    IF k < 0 THEN k := k + bigeven#;
  122.    j := randarray[ii]
  123.   END;
  124.  i := rndknuth(randarray);
  125.  i := rndknuth(randarray);
  126.  i := rndknuth(randarray);
  127. END; { of : PROCEDURE initknuth }
  128.  
  129. FUNCTION random#r : REAL;
  130. {
  131. comment : Returns a REAL pseudo random number in the range 0.0 .. 1.0.
  132.       Requires the definitions needed by rndknuth and  initknuth
  133.       plus the following global :
  134.   VAR    randindex    : INTEGER;
  135. }
  136. BEGIN
  137.  randindex := randindex + 1;
  138.  IF randindex > 55 THEN randindex := rndknuth(randarray);
  139.  random#r  := randarray[randindex]/bigeven#;
  140. END; { of : FUNCTION random#r }
  141.  
  142. FUNCTION random#i : INTEGER;
  143. {
  144. comment : Returns an INTEGER pseudo random number in the range 0.. bigeven#.
  145.       Requires the definitions needed by rndknuth and  initknuth
  146.       plus the following global :
  147.   VAR    randindex    : INTEGER;
  148. }
  149. BEGIN
  150.  randindex := randindex + 1;
  151.  IF randindex > 55 THEN randindex := rndknuth(randarray);
  152.  random#i  := randarray[randindex];
  153. END; { of : FUNCTION random#i }
  154.  
  155. FUNCTION random#sr : REAL;
  156. {
  157. comment : Returns a REAL pseudo random number in the range -1.0 .. +1.0
  158.       Requires the definitions needed by rndknuth and  initknuth
  159.       plus the following global :
  160.   VAR    randindex    : INTEGER;
  161. }
  162. BEGIN
  163.  randindex := randindex + 1;
  164.  IF randindex > 55 THEN randindex := rndknuth(randarray);
  165.  random#sr  := 2.0 * (0.5 - ( randarray[randindex]/bigeven# ));
  166. END; { of : FUNCTION random#sr }
  167.  
  168. FUNCTION random#n : REAL;
  169. {
  170. comment : Returns a REAL number that is randomly selected from a normally
  171. distributed population whose mean is zero and variance (and standard dev.)
  172. is 1.0.
  173. }
  174. VAR
  175.     n    : INTEGER;
  176.     total    : REAL;
  177. BEGIN
  178.  total := -6.0;
  179.  FOR n := 1 TO 12 DO total := total + random#r;
  180.  random#n := total;
  181. END; { of : function randomn }
  182.