From:	SCOVCB::GALVIN_FX    31-JAN-1989 09:02
To:	ARISIA::EVERHART
Subj:	

package body RANDOM is

use UTILITY.REAL_IO; use UTILITY.INT_IO; 

   SEED_MAX : float := 1.0E+08;
--  The exponent is the number of digits in the pseudorandom number.

   TRC_OPTION, RPT_OPTION, DMP_OPTION : UTILITY.OPTION := UTILITY.OFF;

-- Random Stream Creating Operators :
-- ----------------------------------

   -- Random stream creating procedures provide seeds with a default
   -- value of 0, which invokes the random deviate functions to
   -- re-assigning the zero seed with a system-determined seed value.

   procedure CREATE_UNIFORM_FLOAT
            (FLOAT_STREAM : out UNIFORM_FLOAT_RANDOM;
             TITLE : in string;
             MIN_VALUE, MAX_VALUE : in float;
             SEED : in positive := 365437) is

         UFR : UNIFORM_FLOAT_RANDOM;
      begin
         UFR := new UNIFORM_FLOAT_REC;
         UFR.TITLE := UTILITY.STRING_MATCH (TITLE, 8);
         UFR.SEED := float (SEED) / SEED_MAX;
         UFR.MIN_VALUE := MIN_VALUE;
         UFR.DIFF := MAX_VALUE - MIN_VALUE;
         FLOAT_STREAM := UFR;
         if TRC_OPTION = UTILITY.ON then
            UTILITY.OUTPUT_TRACE (SIMULATION.CURRENT_SIMULATION_TIME, 
                            SIMULATION.CURRENT_PROCESS_NAME, 
                            "UNIFORM FLOAT STREAM CREATED", UTILITY.NO_LF);
            put ("     TITLE = "); put (TITLE); 
            put ("     SEED = "); put (SEED); 
            put ("     MIN_VALUE = "); put (MIN_VALUE); 
            put ("     MAX_VALUE = "); put (MAX_VALUE); 
            new_line;
         end if;
      end CREATE_UNIFORM_FLOAT;


   procedure CREATE_UNIFORM_INTEGER
            (INTEGER_STREAM : out UNIFORM_INTEGER_RANDOM;
             TITLE : in string;
             MIN_VALUE, MAX_VALUE : in integer;
             SEED : in positive := 365437) is

         UIR : UNIFORM_INTEGER_RANDOM;
      begin
         UIR := new UNIFORM_INTEGER_REC;
         UIR.TITLE := UTILITY.STRING_MATCH (TITLE, 8);
         UIR.SEED := float (SEED) / SEED_MAX;
         UIR.MIN_VALUE := MIN_VALUE;
         UIR.DIFF := float (MAX_VALUE - MIN_VALUE);
         INTEGER_STREAM := UIR;
         if TRC_OPTION = UTILITY.ON then
            UTILITY.OUTPUT_TRACE (SIMULATION.CURRENT_SIMULATION_TIME, 
                            SIMULATION.CURRENT_PROCESS_NAME, 
                            "UNIFORM INTEGER STREAM CREATED", UTILITY.NO_LF);
            put ("     TITLE = "); put (TITLE); 
            put ("     SEED = "); put (SEED); 
            put ("     MIN_VALUE = "); put (MIN_VALUE); 
            put ("     MAX_VALUE = "); put (MAX_VALUE); 
            new_line;
         end if;
      end CREATE_UNIFORM_INTEGER;


   procedure CREATE_UNIFORM_BOOLEAN
            (BOOLEAN_STREAM : out UNIFORM_BOOLEAN_RANDOM;
             TITLE : in string;
             PROPORTION_TRUE : in float;
             SEED : in positive := 365437) is

         UBR : UNIFORM_BOOLEAN_RANDOM;
      begin
         UBR := new BOOLEAN_REC;
         UBR.TITLE := UTILITY.STRING_MATCH (TITLE, 8);
         UBR.SEED := float (SEED) / SEED_MAX;
         UBR.PROPORTION_TRUE := PROPORTION_TRUE;
         BOOLEAN_STREAM := UBR;
         if TRC_OPTION = UTILITY.ON then
            UTILITY.OUTPUT_TRACE (SIMULATION.CURRENT_SIMULATION_TIME, 
                            SIMULATION.CURRENT_PROCESS_NAME, 
                            "BOOLEAN STREAM CREATED", UTILITY.NO_LF);
            put ("     TITLE = "); put (TITLE); 
            put ("     SEED = "); put (SEED); 
            put ("     PROPORTION_TRUE = "); put (PROPORTION_TRUE); 
            new_line;
         end if;
      end CREATE_UNIFORM_BOOLEAN;


   procedure CREATE_NORMAL
            (NORMAL_STREAM : out NORMAL_RANDOM;
             TITLE : in string;
             MEAN, STD_DEV, MIN_VALUE, MAX_VALUE : in float;
             SEED : in positive := 365437) is

         NR : NORMAL_RANDOM;
      begin
         NR := new NORMAL_REC;
         NR.TITLE := UTILITY.STRING_MATCH (TITLE, 8);
         NR.SEED := float (SEED) / SEED_MAX;
         NR.MEAN := MEAN;
         NR.STD_DEV := STD_DEV;
         NR.MIN_VALUE := MIN_VALUE;
         NR.MAX_VALUE := MAX_VALUE;
         NORMAL_STREAM := NR;
         if TRC_OPTION = UTILITY.ON then
            UTILITY.OUTPUT_TRACE (SIMULATION.CURRENT_SIMULATION_TIME, 
                            SIMULATION.CURRENT_PROCESS_NAME, 
                            "NORMAL STREAM CREATED", UTILITY.NO_LF);
            put ("     TITLE = "); put (TITLE); 
            put ("     SEED = "); put (SEED); 
            put ("     MEAN = "); put (MEAN); 
            put ("     STD_DEV = "); put (STD_DEV); 
            put ("     MIN_VALUE = "); put (MIN_VALUE); 
            put ("     MAX_VALUE = "); put (MAX_VALUE); 
            new_line;
         end if;
      end CREATE_NORMAL;


   procedure CREATE_POISSON
            (POISSON_STREAM : out POISSON_RANDOM;
             TITLE : in string;
             MEAN, MIN_VALUE, MAX_VALUE : in integer;
             SEED : in positive := 365437) is

         PR : POISSON_RANDOM;
      begin
         PR := new POISSON_REC;
         PR.TITLE := UTILITY.STRING_MATCH (TITLE, 8);
         PR.SEED := FLOAT(SEED) / SEED_MAX;
         PR.MEAN := MEAN;
         PR.MIN_VALUE := MIN_VALUE;
         PR.MAX_VALUE := MAX_VALUE;
         POISSON_STREAM := PR;
         if TRC_OPTION = UTILITY.ON then
            UTILITY.OUTPUT_TRACE (SIMULATION.CURRENT_SIMULATION_TIME, 
                            SIMULATION.CURRENT_PROCESS_NAME, 
                            "POISSON STREAM CREATED", UTILITY.NO_LF);
            put ("     TITLE = "); put (TITLE); 
            put ("     SEED = "); put (SEED); 
            put ("     MEAN = "); put (MEAN); 
            put ("     MIN_VALUE = "); put (MIN_VALUE); 
            put ("     MAX_VALUE = "); put (MAX_VALUE); 
            new_line;
         end if;
      end CREATE_POISSON;


   procedure CREATE_EXPONENTIAL
            (EXPONENTIAL_STREAM : out EXPONENTIAL_RANDOM;
             TITLE : in string;
             MEAN, MIN_VALUE, MAX_VALUE : in float;
             SEED : in positive := 365437) is

         ER : EXPONENTIAL_RANDOM;
      begin
         ER := new EXPONENTIAL_REC;
         ER.TITLE := UTILITY.STRING_MATCH (TITLE, 8);
         ER.SEED := float (SEED) / SEED_MAX;
         ER.MEAN := MEAN;
         ER.MIN_VALUE := MIN_VALUE;
         ER.MAX_VALUE := MAX_VALUE;
         EXPONENTIAL_STREAM := ER;
         if TRC_OPTION = UTILITY.ON then
            UTILITY.OUTPUT_TRACE (SIMULATION.CURRENT_SIMULATION_TIME, 
                            SIMULATION.CURRENT_PROCESS_NAME, 
                            "EXPONENTIAL STREAM CREATED", UTILITY.NO_LF);
            put ("     TITLE = "); put (TITLE); 
            put ("     SEED = "); put (SEED); 
            put ("     MEAN = "); put (MEAN); 
            put ("     MIN_VALUE = "); put (MIN_VALUE); 
            put ("     MAX_VALUE = "); put (MAX_VALUE); 
            new_line;
         end if;
      end CREATE_EXPONENTIAL;


   procedure CREATE_ERLANG
            (ERLANG_STREAM : out ERLANG_RANDOM;
             TITLE : in string;
             NUM_OF_DEVIATES : in positive;
             MEAN, MIN_VALUE, MAX_VALUE : in float;
             SEED : in positive := 365437) is

         ER : ERLANG_RANDOM;
      begin
         ER := new ERLANG_REC;
         ER.TITLE := UTILITY.STRING_MATCH (TITLE, 8);
         ER.SEED := float (SEED) / SEED_MAX;
         ER.NUMBER_OF_DEVIATES := NUM_OF_DEVIATES;
         ER.MEAN := MEAN;
         ER.MIN_VALUE := MIN_VALUE;
         ER.MAX_VALUE := MAX_VALUE;
         ERLANG_STREAM := ER;
         if TRC_OPTION = UTILITY.ON then
            UTILITY.OUTPUT_TRACE (SIMULATION.CURRENT_SIMULATION_TIME, 
                            SIMULATION.CURRENT_PROCESS_NAME, 
                            "ERLANG STREAM CREATED", UTILITY.NO_LF);
            put ("     TITLE = "); put (TITLE); 
            put ("     SEED = "); put (SEED); 
            put ("     NUM_OF_DEVIATES = "); put (NUM_OF_DEVIATES); 
            put ("     MEAN = "); put (MEAN); 
            put ("     MIN_VALUE = "); put (MIN_VALUE); 
            put ("     MAX_VALUE = "); put (MAX_VALUE); 
            new_line;
         end if;
      end CREATE_ERLANG;



-- Sampling Operators :
-- --------------------

   --  BASIC PSEUDORANDOM NUMBER GENERATOR from
   --  SIMULATION:  PRINCIPLES and METHODS by
   --  Wayne Graybeal and Udo W. Pooch, sec. 4.2.5, p. 84.

   procedure RND (SEED : in out float) is
         CONST : float := 291.0;
         TERM  : float;
      begin
         TERM :=  CONST * SEED;
         SEED :=  TERM - float (integer (TERM + 0.5) - 1);
      end RND;


   function UNIFORM_FLOAT
           (FLOAT_STREAM : in UNIFORM_FLOAT_RANDOM) return float is
         RANDOM : float;
      begin
         if FLOAT_STREAM = null then
            raise RANDOM_VARIABLE_DOES_NOT_EXIST;
         end if;
         RND (FLOAT_STREAM.SEED);
         RANDOM := FLOAT_STREAM.SEED * FLOAT_STREAM.DIFF + 
            FLOAT_STREAM.MIN_VALUE;
         if TRC_OPTION = UTILITY.ON then
            UTILITY.OUTPUT_TRACE (SIMULATION.CURRENT_SIMULATION_TIME, 
                            SIMULATION.CURRENT_PROCESS_NAME, 
                            "UNIFORM FLOAT RANDOM", UTILITY.NO_LF);
            put (" computed for "); put (FLOAT_STREAM.TITLE); put (" = ");
            put (RANDOM);
            new_line;
         end if;
         return RANDOM;
      end UNIFORM_FLOAT;


   function UNIFORM_INTEGER
           (INTEGER_STREAM : in UNIFORM_INTEGER_RANDOM) return integer is
         RANDOM : integer;

      begin
         if INTEGER_STREAM = null then
            raise RANDOM_VARIABLE_DOES_NOT_EXIST;
         end if;
         RND (INTEGER_STREAM.SEED);
         RANDOM := integer (INTEGER_STREAM.SEED * INTEGER_STREAM.DIFF)
                            + INTEGER_STREAM.MIN_VALUE;
         if TRC_OPTION = UTILITY.ON then
            UTILITY.OUTPUT_TRACE (SIMULATION.CURRENT_SIMULATION_TIME, 
                            SIMULATION.CURRENT_PROCESS_NAME, 
                            "UNIFORM INTEGER RANDOM", UTILITY.NO_LF);
            put (" computed for "); put (INTEGER_STREAM.TITLE); put (" = ");
            put (RANDOM);
            new_line;
         end if;
         return RANDOM;
      end UNIFORM_INTEGER;


   function UNIFORM_BOOLEAN
           (BOOLEAN_STREAM : in UNIFORM_BOOLEAN_RANDOM) return boolean is
         RANDOM : float;
         LOGICAL : boolean;
      begin
         if BOOLEAN_STREAM = null then
            raise RANDOM_VARIABLE_DOES_NOT_EXIST;
         end if;
         RND (BOOLEAN_STREAM.SEED);
         RANDOM := BOOLEAN_STREAM.SEED;
         if RANDOM >= BOOLEAN_STREAM.PROPORTION_TRUE then
            LOGICAL := TRUE;
         else
            LOGICAL := FALSE;
         end if;
         if TRC_OPTION = UTILITY.ON then
            UTILITY.OUTPUT_TRACE (SIMULATION.CURRENT_SIMULATION_TIME, 
                            SIMULATION.CURRENT_PROCESS_NAME, 
                            "BOOLEAN RANDOM", UTILITY.NO_LF);
            put (" computed for "); put (BOOLEAN_STREAM.TITLE); put (" = ");
            if LOGICAL = TRUE then
               put ("TRUE");
            else
               put ("FALSE");
            end if;
            new_line;
         end if;
         return LOGICAL;
      end UNIFORM_BOOLEAN;


   function NORMAL (NORMAL_STREAM : in NORMAL_RANDOM) return float is
      --
      -- Using 'Polar Method', compute a normal deviate :
      --
      --   {Repeat}
      --    1. Generate two random deviates, U_1 & U_2, uniformly
      --       distributed between 0 & 1.
      --    2. Translate U_1 & U_2 to yield two deviates uniformly
      --       distributed bwteen  -1  and  +1 :
      --          Let   U_1  <--  2 * U_1  -  1   and
      --                U_2  <--  2 * U_2  -  1.
      --    3. Let   S  <--  U_1 ** 2  +  U_2 ** 2.
      --    4. If  S >= 1 then jump back to step 1.
      --   {Until  S < 1}
      --    5. Now two potential normal deviates with
      --       mean 0  and  standard deviation 1 are :
      --          X_1  =  U_1  *  square_root_of (-2 * ln (S) / S)  and
      --          X_2  =  U_2  *  square_root_of (-2 * ln (S) / S).
      --       
         NORMAL_RAND_NO, SUM_OF_SQUARES, 
         UNIFORM_RAND_NO_1, UNIFORM_RAND_NO_2 : float;

      begin
         --
         if NORMAL_STREAM = null then
            raise RANDOM_VARIABLE_DOES_NOT_EXIST;
         end if;
         SUM_OF_SQUARES := 1.0;
         while  SUM_OF_SQUARES >= 1.0  loop
            RND (NORMAL_STREAM.SEED);
            UNIFORM_RAND_NO_1 := NORMAL_STREAM.SEED;
            RND (NORMAL_STREAM.SEED);
            UNIFORM_RAND_NO_2 := NORMAL_STREAM.SEED;
            UNIFORM_RAND_NO_1 := UNIFORM_RAND_NO_1 + UNIFORM_RAND_NO_1 - 1.0;
            UNIFORM_RAND_NO_2 := UNIFORM_RAND_NO_2 + UNIFORM_RAND_NO_2 - 1.0;
            SUM_OF_SQUARES := UNIFORM_RAND_NO_1 * UNIFORM_RAND_NO_1
                            + UNIFORM_RAND_NO_2 * UNIFORM_RAND_NO_2;
         end loop;  -- SUM_OF_SQUARES < 1
         --
         NORMAL_RAND_NO := UNIFORM_RAND_NO_1
            * ((-2.0 * MATH.LN (SUM_OF_SQUARES) / SUM_OF_SQUARES) ** 0.5);
         -- Transform normal deviate to correspond
         -- with actual distribution parameters :
         NORMAL_RAND_NO := NORMAL_RAND_NO * NORMAL_STREAM.STD_DEV
                         + NORMAL_STREAM.MEAN;
         --
         -- Assure that normal deviate computed is
         -- within the limits of possible values :
         if NORMAL_RAND_NO < NORMAL_STREAM.MIN_VALUE then
            NORMAL_RAND_NO := NORMAL_STREAM.MIN_VALUE;
         elsif NORMAL_RAND_NO > NORMAL_STREAM.MAX_VALUE then
            NORMAL_RAND_NO := NORMAL_STREAM.MAX_VALUE;
         end if;
         --
         --
         if TRC_OPTION = UTILITY.ON then
            UTILITY.OUTPUT_TRACE (SIMULATION.CURRENT_SIMULATION_TIME, 
                            SIMULATION.CURRENT_PROCESS_NAME, 
                            "NORMAL RANDOM", UTILITY.NO_LF);
            put (" computed for "); put (NORMAL_STREAM.TITLE); put (" = ");
            put (NORMAL_RAND_NO);
            new_line;
         end if;
         return NORMAL_RAND_NO;
      end NORMAL;


   function POISSON
           (POISSON_STREAM : in POISSON_RANDOM) return integer is
      --
      -- {Repeat}
      --  1. For e = 2.71828..., compute exponential part & inititialize :
      --        p <-- e ** (- poisson mean),  q <-- 0,  &  N <-- 0.
      --  2. Get unifrom random deviate :  q <-- q * {random number}.
      --  3. If  q >= p  then
      --        increment N by 1, &
      --        jump back to step 2;
      --     else if  q < p  then
      --        {random deviate with poisson distribution} <-- N;
      --     end if;
      --  {Until  q < p}
      --
         POISSON_RAND_NO : integer;
         EXPON_PART, RAND_NO, RANDOM_PRODUCT : float;
      --
      begin
         --
         if POISSON_STREAM = null then
            raise RANDOM_VARIABLE_DOES_NOT_EXIST;
         end if;
         EXPON_PART := MATH.EXP (- float (POISSON_STREAM.MEAN));
         POISSON_RAND_NO := 0;
         RANDOM_PRODUCT := 1.0;
         while  EXPON_PART <= RANDOM_PRODUCT  loop
            RND (POISSON_STREAM.SEED);
            RANDOM_PRODUCT := RANDOM_PRODUCT * POISSON_STREAM.SEED;
            POISSON_RAND_NO := POISSON_RAND_NO + 1;
         end loop;
         --
         -- Transform poisson deviate to correspond with
         -- actual distribution required :
         POISSON_RAND_NO := POISSON_RAND_NO + 
                            integer (POISSON_STREAM.MIN_VALUE);
         if POISSON_RAND_NO > POISSON_STREAM.MAX_VALUE then
            POISSON_RAND_NO := POISSON_STREAM.MAX_VALUE;
         end if;
         --
         --
         if TRC_OPTION = UTILITY.ON then
            UTILITY.OUTPUT_TRACE (SIMULATION.CURRENT_SIMULATION_TIME, 
                            SIMULATION.CURRENT_PROCESS_NAME, 
                            "POISSON RANDOM", UTILITY.NO_LF);
            put (" computed for "); put (POISSON_STREAM.TITLE); put (" = ");
            put (POISSON_RAND_NO);
            new_line;
         end if;
         return POISSON_RAND_NO;
      end POISSON;


   function EXPONENTIAL
           (EXPONENTIAL_STREAM : in EXPONENTIAL_RANDOM) return float is
         RANDOM : float;
      begin
         if EXPONENTIAL_STREAM = null then
            raise RANDOM_VARIABLE_DOES_NOT_EXIST;
         end if;
         RND (EXPONENTIAL_STREAM.SEED);
         RANDOM := - EXPONENTIAL_STREAM.MEAN * 
                   MATH.LN (EXPONENTIAL_STREAM.SEED);
         if TRC_OPTION = UTILITY.ON then
            UTILITY.OUTPUT_TRACE (SIMULATION.CURRENT_SIMULATION_TIME, 
                            SIMULATION.CURRENT_PROCESS_NAME, 
                            "EXPONENTIAL RANDOM", UTILITY.NO_LF);
            put (" computed for "); put (EXPONENTIAL_STREAM.TITLE); 
            put (" = "); put (RANDOM);
            new_line;
         end if;
         return RANDOM;
      end EXPONENTIAL;


   function ERLANG (ERLANG_STREAM : in ERLANG_RANDOM) return float is
--
-- An Erlang function is defined by the parameters ALPHA and BETA.
-- The NUMBER_OF_DEVIATES is ALPHA.  The MEAN is (ALPHA * BETA).
--
-- Erlang(x) = e**(-x/BETA) * x**(ALPHA-1) / ((ALPHA-1)! * BETA**ALPHA)
--
-- y = 1-exp(-x/BETA)  x=-BETA ln y
--
         RANDOM : float := 1.0;
      begin
         if ERLANG_STREAM = null then
            raise RANDOM_VARIABLE_DOES_NOT_EXIST;
         end if;
         for i in 1..ERLANG_STREAM.NUMBER_OF_DEVIATES loop
            RND (ERLANG_STREAM.SEED);
            RANDOM := RANDOM * ERLANG_STREAM.SEED;
         end loop;
         RANDOM := - ERLANG_STREAM.MEAN / 
                   float (ERLANG_STREAM.NUMBER_OF_DEVIATES) * MATH.LN (RANDOM);
         if RANDOM < ERLANG_STREAM.MIN_VALUE then
            RANDOM := ERLANG_STREAM.MIN_VALUE;
         elsif RANDOM > ERLANG_STREAM.MAX_VALUE then
            RANDOM := ERLANG_STREAM.MAX_VALUE;
         end if;
         if TRC_OPTION = UTILITY.ON then
            UTILITY.OUTPUT_TRACE (SIMULATION.CURRENT_SIMULATION_TIME, 
                            SIMULATION.CURRENT_PROCESS_NAME, 
                            "ERLANG RANDOM", UTILITY.NO_LF);
            put (" computed for "); put (ERLANG_STREAM.TITLE); put (" = ");
            put (RANDOM);
            new_line;
         end if;
         return RANDOM;
      end ERLANG;


-- Trace Operator :
-- ------------------------------------

   procedure SET_RANDOM_TRACE_OPTION (TRACE_OPTION : in UTILITY.OPTION) is
      begin
         TRC_OPTION := TRACE_OPTION;
         if TRC_OPTION = UTILITY.ON then
            UTILITY.OUTPUT_TRACE (SIMULATION.CURRENT_SIMULATION_TIME, 
                            SIMULATION.CURRENT_PROCESS_NAME, 
                            "SET RANDOM TRACE OPTION ON", UTILITY.LF);
         end if;
      end;

end RANDOM;
