home *** CD-ROM | disk | FTP | other *** search
/ Oakland CPM Archive / oakcpm.iso / sigm / vol148 / rndknuth.lib < prev    next >
Encoding:
Text File  |  1984-04-29  |  6.1 KB  |  184 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(seed : INTEGER);
  14. FUNCTION  rndknuth  : 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. Four 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 the above 3 should be used in the one program.)
  57. The fourth generator (randomn) returns a random REAL numberfrom a NORMAL 
  58. distribution with mean = zero and variance = 1.This FUNCTION requires
  59. FUNCTION randomr .
  60. All 4 generators require the FUNCTION rndknuth,the PROCEDURE initknuth
  61. and a collection of global definitions as listed above.
  62.  
  63. Modified for Pascal/Z by G.M.Acland : University of Pennsylvania,1982.
  64. }
  65.  
  66. FUNCTION rndknuth : byte;
  67. {
  68. comment : fills the array "randarray" with 55 pseudo random INTEGERS in
  69.       the range 0..bigeven#. Knuth originally specified 10^9 for
  70.       bigeven# . For Pascal/Z the best number = MAXINT - 1.
  71.       Requires the following definitions globally :
  72.  CONST    bigeven#    = MAXINT - 1;
  73.  TYPE    "intarray"    = ARRAY[1..55] OF INTEGER;
  74.     "byte"        = 0..255;
  75.  VAR    "randarray"    : "intarray";
  76.        Returns the value 1 ( for reinitializing index to "randarray").
  77. }
  78. VAR
  79.     i,j,k    :  INTEGER;
  80. BEGIN
  81.  FOR i := 1 TO 55 DO
  82.   BEGIN
  83.     k := I + 31;
  84.     IF k > 55 THEN k := k - 55;
  85.     j := randarray[i] - randarray[k];
  86.     IF j < 0 THEN j := j + bigeven#;
  87.     randarray[i] := j
  88.   END;
  89.  rndknuth := 1;
  90. END; { of : FUNCTION rndknuth }
  91.  
  92. PROCEDURE initknuth(seed : INTEGER);
  93. {
  94. comment : Initializes the GLOBAL randarray. Has the same requirements as
  95.           rndknuth, which FUNCTION  is called by initknuth, plus the input
  96.           value : "seed" : this may be a zero,one or any other positive
  97.       INTEGER value.A useful technic when you want to use a "random"
  98.       seed is to create an integer from the time of day , if you have
  99.       it available to your computer.
  100.  e.g.: procedure initseed;
  101.        BEGIN
  102.         timestring := '  :  :  ';
  103.         seed := 0;
  104.         time(timestring);
  105.         FOR i := 1 TO 8 DO seed := seed + ORD(timestring[i]);
  106.        END;
  107.          Note that the GLOBAL randindex is initialized within this routine.
  108.     This was not initialized in the previous version distributed thru
  109.     PZUG.
  110. }
  111. VAR
  112.     i,ii,j,k : INTEGER;
  113. BEGIN
  114.  randindex    := 1;
  115.  randarray[55]    := seed;
  116.  j        := seed;
  117.  k        := 1;
  118.  FOR i := 1 TO 54 DO
  119.   BEGIN
  120.    ii := (21 * i) MOD 55;
  121.    randarray[ii] := k;
  122.    k := j - k;
  123.    IF k < 0 THEN k := k + bigeven#;
  124.    j := randarray[ii]
  125.   END;
  126.  i := rndknuth;
  127.  i := rndknuth;
  128.  i := rndknuth;
  129. END; { of : PROCEDURE initknuth }
  130.  
  131. FUNCTION random#r : REAL;
  132. {
  133. comment : Returns a REAL pseudo random number in the range 0.0 .. 1.0.
  134.       Requires the definitions needed by rndknuth and  initknuth
  135.       plus the following global :
  136.   VAR    randindex    : INTEGER;
  137. }
  138. BEGIN
  139.  randindex := randindex + 1;
  140.  IF randindex > 55 THEN randindex := rndknuth;
  141.  random#r  := randarray[randindex]/bigeven#;
  142. END; { of : FUNCTION random#r }
  143.  
  144. FUNCTION random#i : INTEGER;
  145. {
  146. comment : Returns an INTEGER pseudo random number in the range 0.. bigeven#.
  147.       Requires the definitions needed by rndknuth and  initknuth
  148.       plus the following global :
  149.   VAR    randindex    : INTEGER;
  150. }
  151. BEGIN
  152.  randindex := randindex + 1;
  153.  IF randindex > 55 THEN randindex := rndknuth;
  154.  random#i  := randarray[randindex];
  155. END; { of : FUNCTION random#i }
  156.  
  157. FUNCTION random#sr : REAL;
  158. {
  159. comment : Returns a REAL pseudo random number in the range -1.0 .. +1.0
  160.       Requires the definitions needed by rndknuth and  initknuth
  161.       plus the following global :
  162.   VAR    randindex    : INTEGER;
  163. }
  164. BEGIN
  165.  randindex := randindex + 1;
  166.  IF randindex > 55 THEN randindex := rndknuth;
  167.  random#sr  := 2.0 * (0.5 - ( randarray[randindex]/bigeven# ));
  168. END; { of : FUNCTION random#sr }
  169.  
  170. FUNCTION random#n : REAL;
  171. {
  172. comment : Returns a REAL number that is randomly selected from a normally
  173. distributed population whose mean is zero and variance (and standard dev.)
  174. is 1.0.
  175. }
  176. VAR
  177.     n    : INTEGER;
  178.     total    : REAL;
  179. BEGIN
  180.  total := -6.0;
  181.  FOR n := 1 TO 12 DO total := total + random#r;
  182.  random#n := total;
  183. END; { of : function randomn }
  184.