--------------------------------------------------------------------------
-- The following is an implementation of a "universal" random number    --
-- generator algorithm developed by Dr. George Marsaglia of the         --
-- Supercomputer Computations Research Institute (SCRI) at Florida      -- 
-- State University.  This generator has a period of ~2**144 and has    --
-- been tailored for reproducibility in all CPU's with at least         --
-- 16 bit integer arithmetic and 24 bit floating point.  This algorithm --
-- does not generate random numbers < 2**-24.  At the end of this file  --
-- you will find a self test program that checks generated results      --
-- against known expected results and reports any inaccuracies.         --
--                                                                      --
-- Further references: "Toward a Universal Random Number Generator",    --
-- appearing in the Journal of The American Statistical Association.    --
--                                                                      --
-- This code appeared in the March/April publication of SIGAda's        --
-- Ada Letters and is considered public domain.  PCK                    --
--------------------------------------------------------------------------

package U_Rand is
    M1 : constant := 179 ;
    M2 : constant := M1 - 10 ;

    subtype SEED_RANGE_1 is integer range 1..M1-1 ;
    subtype SEED_RANGE_2 is integer range 1..M2-1 ;

    Default_I : constant SEED_RANGE_1 := 12 ;
    Default_J : constant SEED_RANGE_1 := 34 ;
    Default_K : constant SEED_RANGE_1 := 56 ;
    Default_L : constant SEED_RANGE_1 := 78 ;

    procedure Start(New_I : SEED_RANGE_1 := Default_I ;
                    New_J : SEED_RANGE_1 := Default_J ;
                    New_K : SEED_RANGE_1 := Default_K ;
                    New_L : SEED_RANGE_2 := Default_L) ;

    function Next return FLOAT ;
end U_Rand ;

package body U_Rand is
    M3     : constant := 97 ;
    Init_C : constant := 362436.0/16777216.0 ;
    CD     : constant := 7654321.0/16777216.0 ;
    CM     : constant := 16777213.0/16777216.0 ;

    subtype RANGE_1 is INTEGER range 0..M1-1 ;
    subtype RANGE_2 is INTEGER range 0..M2-1 ;
    subtype RANGE_3 is INTEGER range 1..M3 ;

    I, J, K : RANGE_1 ;
    NI, NJ : INTEGER ;
    L : RANGE_2 ;
    C : FLOAT ;
    U : array(RANGE_3) of FLOAT ;

    procedure Start(New_I : SEED_RANGE_1 := Default_I ;
                    New_J : SEED_RANGE_1 := Default_J ;
                    New_K : SEED_RANGE_1 := Default_K ;
                    New_L : SEED_RANGE_2 := Default_L) is
        S, T : FLOAT ;
        M : RANGE_1 ;
    begin
        I := New_I ; J := New_J ; K := New_K ; L := New_L ;
        NI := RANGE_3'last ;
        NJ := (RANGE_3'last/3) + 1 ;
        C := Init_C ;

        for II in RANGE_3 loop 
            S := 0.0 ; T := 0.5 ;
            for JJ in 1..24 loop 
                M := (((J * I) mod M1) * K) mod M1 ;
                I := J ; J := K ; K := M ;
                L := (53 * L + 1) mod M2 ;
                if ((L * M) mod 64) >= 32 then
                    S := S + T ;
                end if ;
                T := 0.5 * T ;
            end loop ;
            U(II) := S ;
        end loop ;
    end Start ;

    function Next return FLOAT is
        Temp : FLOAT ;
    begin
        Temp := U(NI) - U(NJ) ;
        if Temp < 0.0 then
            Temp := Temp + 1.0 ;
        end if ;
        U(NI) := Temp ;
        NI := NI - 1 ;
        if NI = 0 then
            NI := RANGE_3'last ;
        end if ;
        NJ := NJ - 1 ;
        if NJ = 0 then
            NJ := RANGE_3'last ;
        end if ;
        C := C - CD ;
        if C < 0.0 then
            C := C + CM ;
        end if ;
        Temp := Temp - C ;
        if Temp < 0.0 then
            Temp := Temp + 1.0 ;
        end if ;
        return Temp ;
    end Next ;

begin
    Start ; -- initialize table U
end U_Rand ;
------------------------------------------------------------------------
-- test program follows
--
--
--with Text_IO ; use Text_IO ;
--with U_Rand ;
--procedure Test_Urand is
--    package Float_IO is new Text_IO.Float_IO(FLOAT) ;
--    use Float_IO ;
--    Test_File : Text_IO.File_Type ;
--    File_Name : STRING(1..21) := "URand_Test_Output.dat" ;
--    subtype A_RESULT is STRING(1..11) ;
--    type RESULTS_TYPE is array(1..6) of A_RESULT ;
--    Actual_Results : RESULTS_TYPE ;
--    EOS : NATURAL ;
--    Results_OK : BOOLEAN := True ;
--    Expected_Results : RESULTS_TYPE := (1=>" 6533892.00",
--                                        2=>"14220222.00",
--                                        3=>" 7275067.00",
--                                        4=>" 6172232.00",
--                                        5=>" 8354498.00",
--                                        6=>"10633180.00") ;

--    Rnum : FLOAT ;
--begin
--    for I in 1..20000 loop
--        Rnum := U_Rand.Next ;
--    end loop ;

--    Text_IO.Create(File=>Test_File, Name=>File_Name) ;
--    for I in 1..6 loop
--        Put(File=>Test_File,
--            Item=>(2.0**24)*U_Rand.Next,
--            Fore=>8, Aft=>2, Exp=>0) ; 
--        New_Line(Test_File) ;
--    end loop ;

--    Text_IO.Close(File=>Test_File) ;
--    Text_IO.Open(File=>Test_File, Mode=>In_File, Name=>File_Name) ;
--    Put("Expected Results    ") ; Put_Line("Actual Results") ; 

--    for I in 1..6 loop
--        Get_Line(File=>Test_File, Item=>Actual_Results(I), Last=>EOS) ;
--        Skip_Line(Test_File) ;
--        Put(Expected_Results(I)) ; Put("         ") ;
--        Put_Line(Actual_Results(I)) ;
--        if Actual_Results(I) /= Expected_Results(I) then
--            Results_OK := False ;
--        end if ;
--    end loop ;

--    New_Line(2) ; 
--    if not Results_OK then
--        Put_Line("!! CAUTION!! Random Number Generator inconsistent on this inplementation") ;
--    else
--        Put_Line("Random Number Generator Consistent") ;
--    end if ;
--end Test_Urand ;