program gentst(gentstp,data,genhisp,output);
(* generate a bunch of numbers, each the sum of random numbers
   to test genhis.  Tom Schneider *)

label 1; (* end of program *)

const
(* begin module version *)
version = 3.13; (* of gentst.p 1994 sep 5
origin 1985 dec 19 tds *)
(* end module version *)

(* begin module describe.gentst *)
(*
name
   gentst: test random generator

synopsis
   gentst(gentstp:in, data: out, output: out) 

files 
   gentstp: parameter file controlling the program.
      Three numbers, one per line:
         seed: random seed to start the process
         total: the number of numbers to generate
         components: the number of random numbers between 0 and 1 to add
            together to generate the total
   data: the input file for genhis.  this is a set of numbers which
      should have gaussian distribution if the random number generator
      is a reasonable one.  It will be N(0,1), a normal distribution
      with mean 0 and standard deviation 1.
   genhisp: control file for genhis
   output: messages to the user

description 
   test of a random number generator by creating a gaussian
   distribution of numbers for plotting by genhis

example
  seed := 0.5;
  total := 10000;
  components:= 100;

see also
   tstrnd.p, genhis.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.gentst *)

var
(* begin module gentst.var *)
   gentstp,data, genhisp: text; (* output files *)
(* end module gentst.var *)

(* 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 gentst.themain *)
procedure themain(var data, genhisp: text);
(* the main procedure of gentst *)
const
  xnumber = 20; (* number of points on x axis to plot *)
var
  total: integer; (* total numbers to make *)
  components: integer; (* each number is made of this many random numbers *)
  indexfortotal,indexforcomponents: integer; (* indices *)
  seed: real; (* the current random number seed *)
  sum: real; (* the sum of the seeds *)
begin
  writeln(output,'gentst ',version:4:2);

  reset(gentstp);
  readln(gentstp,seed);
  readln(gentstp,total);
  readln(gentstp,components);

  rewrite(genhisp);
  (* set x scale to fit a page: *)
  writeln(genhisp,'x n ', xnumber:1);
  writeln(genhisp,'p g'); (* plot gaussian *)

  rewrite(data);
  writeln(data,'* gentst ',version:4:2);
  writeln(data,'* ',seed:10:9, ' random seed');
  writeln(data,'* ',total:1,' numbers, each the sum of');
  writeln(data,'* ',components:1, ' random numbers');

  for indexfortotal :=1 to total do begin
    sum := 0;
    for indexforcomponents := 1 to components do begin
       random(seed);
       sum := sum + seed;
    end;
    writeln(data,(sum-components/2.0) * sqrt(12.0/components));
  end
end;
(* end module gentst.themain *)

begin
   themain(data,genhisp);
1:end.