-- Author: William N. Graham
--         Applied Research Labs
--         University of Texas
--         Austin, Texas
--
-- Date:   2-27-91

with TEXT_IO;
use  TEXT_IO;

package MITCHELL_MOORE_KNUTH_RANDOM is

    subtype SEED_INDEX is INTEGER range 1..100;

    procedure INITIALIZE(SI : in SEED_INDEX);
    -- Initialize array of initial seeds.

    function RANDOM return FLOAT;
    -- Return another random number between 0 and 1.

    function CHI_SQUARE( N : POSITIVE;    -- Number of iterations.
                         R : POSITIVE;    -- 0 <= RANDOM_INT < R .
                         S : SEED_INDEX)  -- Initial seed index.
                                        return FLOAT;
    -- Compute a chi-square statistic for the random number generator.
    -- This function is from Robert Sedgewick's ALGORITHMS, 1984.

end MITCHELL_MOORE_KNUTH_RANDOM;


package body MITCHELL_MOORE_KNUTH_RANDOM is

    -- Number of initial seeds required in an array.
    MAX_ARRAY_SIZE : constant POSITIVE := 55;

    -- Last index in the array.
    MAX_INDEX : constant NATURAL := MAX_ARRAY_SIZE - 1;

    -- Modulus number, should be large and even.
    MODULUS : INTEGER;

    -- Modulus inverse = 1 / modulus, this will be a float.
    MODULUS_INV : FLOAT;

    -- Array of initial seeds and previous sequence numbers generated.
    -- The initial seeds must not all be even.
    X : array (0..MAX_INDEX) of NATURAL;

    -- Array index of first number in previous sequence.
    -- This is the index of the oldest number in the sequence.
    FIRST : NATURAL := 0;


    procedure INITIALIZE( SI : in SEED_INDEX ) is

        -- Modulus used for linear method.
        M : constant INTEGER := 100000000;

        -- Square root of M.
        M1 : constant INTEGER := 10000;

        -- Multiplier used for linear method.
        -- This should be of length 1 digit less than M and end in the
        -- sequence ...821 .
        B : constant INTEGER := 31415821;

        -- Multiply two large integers without overflow.
        -- p*q = (10000*p1 + p0) * (10000*q1 + q0)
        --     = 100000000*p1*q1 + 10000*(p1*q0+p0*q1) + p0*q0
        function MULT( P : INTEGER;
                       Q : INTEGER ) return INTEGER is

            P1 : INTEGER;
            P0 : INTEGER;
            Q1 : INTEGER;
            Q0 : INTEGER;

        begin

            P1 := P / M1;
            P0 := P mod M1;
            Q1 := Q / M1;
            Q0 := Q mod M1;

            return( (((P0*Q1+P1*Q0) mod M1) * M1 + P0 * Q0) mod M );

        end MULT;
            
        
        -- Linear Congruential random integer.
        -- N := (N*B+1) mod M.
        function LINEAR_RANDOM( N : INTEGER ) return INTEGER is

        begin

            return( (MULT(N,B) + 1) mod M );

        end LINEAR_RANDOM;

    begin

        -- Half of max integer will ensure against overflow.
        MODULUS := INTEGER'LAST / 2;
        
        -- Ensure MODULUS is even.
        if (MODULUS mod 2) /= 0 then  -- MODULUS Odd.
           MODULUS := MODULUS - 1;  -- Make even.
        end if;

        -- Compute inverse of modulus.
        MODULUS_INV := 1.0 / FLOAT(MODULUS);

        -- The initial seeds must not all be even.
        -- Initialize array with seeds by a weak method.
        -- Use the linear congruential method to initialize array.

        -- Set first seed from seed index.
        X(0) := (SI * 2) + 1;  -- Ensure first seed is odd.

        -- Invoke linear generator to fill out array.
        for I in 0..MAX_INDEX-1 loop  -- Until array filled.
            X(I+1) := LINEAR_RANDOM( X(I) );
        end loop;

    end INITIALIZE;


    function RANDOM return FLOAT is

        -- Will return a random Float number that is greater than
        -- or equal to 0.0 and less than 1.0 .

        -- The next index of the oldest number in the sequence.
        -- This is the index of the first addend.
        FIRST_PLUS_1 : NATURAL;

        -- This is the index of the second addend.
        FIRST_PLUS_32 : NATURAL;

        RANDOM_INT : INTEGER;  -- New integer generated;

    begin

        -- Compute next index of the oldest number in the sequence.
        -- This is also the index of the first addend.
        FIRST_PLUS_1 := FIRST + 1;
        -- Check for index out of bounds.
        if FIRST_PLUS_1 > MAX_INDEX then  -- Adjust.
           FIRST_PLUS_1 := FIRST_PLUS_1 - MAX_ARRAY_SIZE;
        end if;

        -- Compute index of the second addend.
        FIRST_PLUS_32 := FIRST + 32;
        -- Check for index out of bounds.
        if FIRST_PLUS_32 > MAX_INDEX then  -- Adjust.
           FIRST_PLUS_32 := FIRST_PLUS_32 - MAX_ARRAY_SIZE;
        end if;

        -- Compute new random integer.  
        RANDOM_INT := X(FIRST_PLUS_1) + X(FIRST_PLUS_32);
        -- Check for integer out of bounds.
        if RANDOM_INT >= MODULUS then  -- Adjust.
           RANDOM_INT := RANDOM_INT - MODULUS;
        end if;

        X(FIRST) := RANDOM_INT;  -- Assign new integer to the array.

        -- Set new index of the oldest number in the sequence.
        FIRST := FIRST_PLUS_1;

        -- Normalize random number and return it.
        return( FLOAT(RANDOM_INT) * MODULUS_INV );

    end RANDOM;


    function CHI_SQUARE( N : POSITIVE;    -- Number of iterations.
                         R : POSITIVE;    -- 0 <= RANDOM_INT < R .
                         S : SEED_INDEX)  -- Initial seed index.
                                        return FLOAT is

    -- Compute a chi-square statistic for the random number generator.
    -- This function is from Robert Sedgewick's ALGORITHMS, 1984.

    -- Maximum R value.
    RMAX : constant NATURAL := 2000;

    -- Array of frequencies.
    FREQ : array (0..RMAX) of NATURAL;
    
    -- Random Integer generated for the test, 0 <= RANDOM_INT < R .
    RANDOM_INT : NATURAL;

    -- Sum of the squares of the frequencies.
    TOTAL : NATURAL;

    begin

        -- Check for R in range.
        if R > RMAX then  -- R is to large.
           put("R TOO LARGE");
           return(0.0);
        end if;

        -- Initialize random number generator.
        INITIALIZE(S);

        -- Initialize frequency array.
        for J in 0..RMAX loop
            FREQ(J) := 0;
        end loop;

        -- Collect frequencies.
        for J in 1..N loop
            -- Generate random integer, 0 <= RANDOM_INT < R .
            RANDOM_INT := INTEGER(RANDOM * FLOAT(R));        
            FREQ(RANDOM_INT) := FREQ(RANDOM_INT) + 1;
        end loop;

        -- Compute sum of squares.
        TOTAL := 0;  -- Initialize sum of squares.
        for J in 0..R-1 loop
            TOTAL := TOTAL + FREQ(J) * FREQ(J);
        end loop;

        -- Compute chi-square.
        return( (FLOAT(R * TOTAL) / FLOAT(N)) - FLOAT(N) );

    end CHI_SQUARE;
    
begin

    INITIALIZE(1);  -- Initialize array with seeds.

end MITCHELL_MOORE_KNUTH_RANDOM;