--                              -*- Mode: Ada -*- 
-- lecuyer.a --- 
-- Author          : William M. Graham
-- Created On      : Mon Feb 11 08:34:03 1991
-- Last Modified By: William M. Graham
-- Last Modified On: Wed Mar 27 17:14:15 1991
-- Update Count    : 3
-- Status          : Available
-- 

package LECUYER is

-- This package provides a pseudo random number generator that produces
-- FLOAT numbers that are uniformly distributed between 0.0 and 1.0 .
--
-- The random number generator can also be initialized with two positive
-- integer seeds.  The only requirement of these seeds is that they fall
-- in the specified ranges 1..2147483562 and 1..2147483398 .
--
-- This algorithm was developed by Pierre L'ecuyer.  The ADA implementation
-- of this algorithm is taken from the code example given in L'ecuyer's
-- EFFICIENT AND PORTABLE COMBINED RANDOM NUMBER GENERATORS paper in the
-- Communications of the ACM, June 1988, pp. 742-774.  This code is 
-- appropriate for any 32 bit or larger machine.  The period for this 
-- generator is reported to be ~ 2.3 * 10**18 .
--
-- Coded by William N. Graham, Applied Research Labs, 2-11-91.


        subtype SEED_RANGE_1 is INTEGER range 1..2147483562;
        subtype SEED_RANGE_2 is INTEGER range 1..2147483398;

        procedure INITIALIZE(SEED1 : in SEED_RANGE_1;  -- First Seed.
                             SEED2 : in SEED_RANGE_2); -- Second Seed.
        -- Initialize the random number generator with two seeds.

        function UNIFORM return FLOAT;
        -- Random number generator will return a number between 0.0 and 1.0 .
        -- This random number generator combines two Multiplicative
        -- Linear Congruential Generators (MLCG's).  A MLCG can be
        -- characterized as F(S) := (A*S) mod M.

end LECUYER;


package body LECUYER is

        -- F(S) := (A*S + C) mod M
        -- For a Multiplicative Linear Congruential Generator (MLCG),
        -- C will be 0.

        M1 : constant POSITIVE := 2147483563;    -- First Modulus.
        A1 : constant POSITIVE := 40014;         -- First Multiplier.
        Q1 : constant POSITIVE := 53668;         -- First Quotient.
        R1 : constant POSITIVE := 12211;         -- First Remainder.
        M2 : constant POSITIVE := 2147483399;    -- Second Modulus.
        A2 : constant POSITIVE := 40692;         -- Second Multiplier.
        Q2 : constant POSITIVE := 52774;         -- Second Quotient.
        R2 : constant POSITIVE := 3791;          -- Second Remainder.

        S1 : INTEGER := 1;  -- The first current state variable.
        S2 : INTEGER := M2 / 2;  -- The second current state variable.

        M1_MINUS_1 : constant POSITIVE := M1 - 1;    -- M1 - 1 .
        NORMALIZE  : constant FLOAT := 4.656613E-10; -- To normalize result.

        procedure INITIALIZE(SEED1 : in SEED_RANGE_1;     -- First Seed.
                             SEED2 : in SEED_RANGE_2) is  -- Second Seed.
        -- Initialize the random number generator with two seeds.

        begin

            S1 := SEED1;  -- Assign first seed to first state variable.
            S2 := SEED2;  -- Assign second seed to second state variable.

        end INITIALIZE;


        function UNIFORM return FLOAT is
        -- Random number generator will return a number between 0.0 and 1.0 .
        -- This random number generator combines two Multiplicative
        -- Linear Congruential Generators.

        Z : INTEGER;  -- Random integer result.
        K : INTEGER;  -- Quotient.

        begin

            K := S1 / Q1;
            S1 := A1 * (S1 - K * Q1) - K * R1;
            if S1 < 0 then S1 := S1 + M1; end if;

            K := S2 / Q2;
            S2 := A2 * (S2 - K * Q2) - K * R2;
            if S2 < 0 then S2 := S2 + M2; end if;

            Z := S1 - S2;
            if Z < 1 then Z := Z + M1_MINUS_1; end if;

            return( FLOAT(Z) * NORMALIZE );

        end UNIFORM;


end LECUYER;