DataMuseum.dk

top - download

    Length: 4021 (0xfb5)
    Types: TextFile
    Names: »B«

└─⟦85b835f43⟧ Bits:30000549 8mm tape, Rational 1000, Xlib rev 6.00
    └─ ⟦0c20f784e⟧ »DATA« 
        └─⟦1abbe589f⟧ 
            └─⟦306851c02⟧ 
                └─⟦this⟧ 

TextFile

with Unchecked_Deallocation;

with Xlbt_Arithmetic;  
use Xlbt_Arithmetic;

package body Ran1_Package is
------------------------------------------------------------------------------
-- Ran1 returns a uniform random deviate between 0.0 and 1.0.  Uses linear
-- congruential methods using three portable not-as-good random number
-- generators.  Near infinite period with no sensible sequential correlation.
--
-- Based upon:  Ran1, pg 196;
-- Numerical Recipes: The Art Of Scientific Computing
-- W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling
-- Cambridge University Press, 1986
------------------------------------------------------------------------------

    type Float_Array is array (S_Natural range <>) of Float;

    M1  : constant       := 259200;  
    Ia1 : constant       := 7141;  
    Ic1 : constant       := 54773;  
    Rm1 : constant Float := 1.0 / Float (M1);

    M2  : constant       := 134456;  
    Ia2 : constant       := 8121;  
    Ic2 : constant       := 28411;  
    Rm2 : constant Float := 1.0 / Float (M2);

    M3  : constant := 243000;  
    Ia3 : constant := 4561;  
    Ic3 : constant := 51349;

    type Ran1_Data_Rec is  
        record  
            Ix1 : S_Natural;  
            Ix2 : S_Natural;  
            Ix3 : S_Natural;  
            R   : Float_Array (1 .. 97);  
        end record;

    procedure Free_Ran1_Data is new Unchecked_Deallocation (Ran1_Data_Rec,  
                                                            Ran1_Data);

--\f

    function Ran1_Initialize (Iseed : S_Natural) return Ran1_Data is
------------------------------------------------------------------------------
--  Iseed   - Specifies the initial seed value to be passed to the RAN()
--              during our initialization.
-- Called to create and initialize a Ran1_Data area for use with Ran1.
------------------------------------------------------------------------------
        Data : Ran1_Data := new Ran1_Data_Rec;  
    begin

----Seed the first generator.

        Data.Ix1 := (Ic1 - Iseed) mod M1;  
        Data.Ix1 := (Ia1 * Data.Ix1 + Ic1) mod M1;

----Use the first to seed the second and third.

        Data.Ix2 := Data.Ix1 mod M2;  
        Data.Ix1 := (Ia1 * Data.Ix1 + Ic1) mod M1;  
        Data.Ix3 := Data.Ix1 mod M3;

----Fill the table with sequential uniform deviates generated by the first two
--  generators.

        for I in Data.R'Range loop  
            Data.Ix1   := (Ia1 * Data.Ix1 + Ic1) mod M1;  
            Data.Ix2   := (Ia2 * Data.Ix2 + Ic2) mod M2;  
            Data.R (I) := (Float (Data.Ix1) + Float (Data.Ix2) * Rm2) * Rm1;  
        end loop;  
        return Data;

    end Ran1_Initialize;

--\f

    procedure Ran1_Free (Idum : in out Ran1_Data) is
------------------------------------------------------------------------------
--  Idum        - Specifies the data area to free.
-- Called to free up the Ran1_Data area allocated by Ran1_Initialize.  It is
-- OK to call this with a null.
------------------------------------------------------------------------------
    begin

        Free_Ran1_Data (Idum);

    end Ran1_Free;

--\f

    function Ran1 (Idum : Ran1_Data) return Float is
------------------------------------------------------------------------------
--  Idum    - Specifies the Ran1_Data to use in calculating the next random
--              value.
-- Called to calculate the next uniform random deviate between 0.0 and 1.0.
------------------------------------------------------------------------------
        J : S_Long;  
        F : Float;  
    begin

        Idum.Ix1 := (Ia1 * Idum.Ix1 + Ic1) mod M1;  
        Idum.Ix2 := (Ia2 * Idum.Ix2 + Ic2) mod M2;  
        Idum.Ix3 := (Ia3 * Idum.Ix3 + Ic3) mod M3;  
        J        := 1 + (Idum.R'Last * Idum.Ix3) / M3;  
        if J > Idum.R'Last then  
            J := J - 1;  
        end if;  
        F          := Idum.R (J);  
        Idum.R (J) := (Float (Idum.Ix1) + Float (Idum.Ix2) * Rm2) * Rm1;  
        return F;

    end Ran1;

--\f

end Ran1_Package;