program tstrnd(output); 
(* test random generator 
   thomas schneider copyright (c) 1988
*)
  
label 1; (* end of the program *) 
  
const 
(* begin module version *)
version = 1.07; (* of tstrnd.p 2000 May 26
origin 1988 October 17 or earlier *)
(* end module version *)
  
(* begin module describe.tstrnd *)
(*
name
      tstrnd: test random generator

synopsis
      tstrnd(output: out)

files 
      output: The version of tstrnd is printed.  Successful compilation 
         and running of the program indicates that the modules are correct. 

description 
      Test of a compiler and system independent random number generator.

see also
      {gaussian random number generator: } gentst.p

author
      Thomas D. Schneider

bugs
      none known
  
technical notes 
      The constant n in procedure randomtest determines how many times  
      the random number generator will be in a series of tests.  If n 
      is small, the the test will be poor, if it is large then the test may 
      take a long time. 

*)
(* end module describe.tstrnd *)
  
(* begin module tstrnd.const *) 
      debugging = false; (* for debugging purposes *) 
(* end module tstrnd.const *) 
  
(* begin module halt *) 
procedure halt; 
(* stop the program.  the procedure performs a goto to the end of the   
   program.  you must have a label: 
      label 1;
   declared, and also the end of the program must have this label:  
      1: end. 
   examples are in the module libraries.
   this is the only goto in the delila system. *) 
begin 
      writeln(output,' program halt.'); 
      goto 1
end;  
(* end module halt version = 'delmod 6.44 82 aug 29 tds/gds'; *)
  
(* begin module random *) 
procedure random(var seed: real); 
(* random generator 2.  version = 1.00; 1986 December 31
   Test this routine with the program tstrnd.
   written by David Masternarde *)

(* This random number generator is based on a shift register with a single bit
of feedback, as described in Electronics for Neurobiologists, by Brown, Maxfield
and Moraff, MIT press 1973, referencing Random Process Simulation and
Measurement by Korn, McGraw-Hill 1966.  The random seed rand, a number between 0
and 1 exclusive, is converted to an integer between 1 and 2**23-1, inclusive.
This 23-bit number is shifted right one bit and the output of the last (23rd)
bit and the 9th bit are added modulo 2 (exclusive orred) and fed back into the
new first bit.  This is done between 4 and 11 times, depending on the last 3
bits of the original number.  The result is converted back to a real number
between 0 and 1 from which the 23 bit integer can be recovered on the next call.
The 23-bit shift register goes through all 2**23-1 values before repeating;
the repetition frequency of this algorithm could be less or greater depending
on the seed, because of the random number of multiple shifts per call. *)
(* powers of 2 *)
const pow14=16384; pow15=32768; pow22= 4194304; pow23=8388608;
var iseed, (* integer shift register *)
    i, nrep: integer;  (* index, number of times to do shift *)
begin (* random *)  
      iseed := round(seed*pow23);  (* convert to 23 bit number *)
      if (iseed<1) or (iseed>=pow23) then iseed := 1;
      nrep := 4 + iseed mod 8;  (* do it 4 to 11 times based on last 3 bits *)
      for i:= 1 to nrep do
         (* if last bit and 9th bit are equal, feed back a 0, otherwise a 1 *)
	 if( odd(iseed) = ((iseed mod pow15) >= pow14))
	    then iseed := iseed div 2
	    else iseed := (iseed div 2) + pow22;
      seed := iseed/pow23;
end; (* random *) 
(* end module random *) 
  
(* begin module random.test *)  
procedure randomtest(var table: text);  
(* this procedure tests the random procedure by calling it many 
times and printing a table representing the results,  the mean and the  
standard deviation for a series of tests.  this random testing procedure
should be run every time random is used on a new computer system. *)
const 
      n = 10000; (* the number of calls to random per test *) 
      tests = 20; (* the number of tests to apply *)  
      seed = 4.0; (* the seed for generation *) 
var 
      test: integer; (* the number of the current test *) 
      call: integer; (* the number of the current call *) 
      x, (* the random variable *)
      sumx, (* the sum of x during one test *)
      sumxsqd, (* the sum of x squared *) 
      mean, (* the mean of x in one test *) 
      std: (* the standard deviation of x in one test *)  
         real;  
begin (* randomtest *)  
      writeln(table); 
      (* initialize things *) 
      x := seed; (* start x.  note that this is not the first random number *)
      writeln(table, ' test of the random number generator:');  
      writeln(table, ' std. dev is 1/sqrt(12); found by integrating:');
      writeln(table, ' from 0 to 1 of (x- mean)^2 dx, where the mean is 0.5');
      writeln(table, ' ');
      writeln(table, ' ');
      writeln(table); 
      writeln(table,'     mean     std dev'); 
      writeln(table,'    (0.5000)  (0.2887)  (expected values)'); 
      writeln(table,' ===================== ', n:1, ' calls per test'); 
  
      for test := 1 to tests do begin 
         (* initialize values for one test *) 
         sumx := 0; 
         sumxsqd := 0;
  
         for call := 1 to n do begin
            random(x); (* generate the next random number *)
            sumx := sumx + x; 
            sumxsqd := sumxsqd + x*x
         end; 
  
         mean := sumx/n;
         std := sqrt( (n/(n-1)) * (sumxsqd/n - mean * mean) );
         writeln(table,' ', mean:10:4, std:10:4)
      end 
end; (* randomtest *) 
(* end module random.test *)  
  
begin (* tstrnd *)  
(* begin module tstrnd.main *)  
      writeln(output, ' test new random generator ', version:4:2);
  
      randomtest(output); 
(* end module tstrnd.main *)  
1: end. (* tstrnd *)