home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / gnat-2.06-src.tgz / tar.out / fsf / gnat / ada / a-nudira.adb < prev    next >
Text File  |  1996-09-28  |  8KB  |  256 lines
  1. ------------------------------------------------------------------------------
  2. --                                                                          --
  3. --                         GNAT RUNTIME COMPONENTS                          --
  4. --                                                                          --
  5. --         A D A . N U M E R I C S . D I S C R E T E _ R A N D O M          --
  6. --                                                                          --
  7. --                                 B o d y                                  --
  8. --                                                                          --
  9. --                            $Revision: 1.9 $                              --
  10. --                                                                          --
  11. --           Copyright (c) 1992,1993,1994 NYU, All Rights Reserved          --
  12. --                                                                          --
  13. -- The GNAT library is free software; you can redistribute it and/or modify --
  14. -- it under terms of the GNU Library General Public License as published by --
  15. -- the Free Software  Foundation; either version 2, or (at your option) any --
  16. -- later version.  The GNAT library is distributed in the hope that it will --
  17. -- be useful, but WITHOUT ANY WARRANTY;  without even  the implied warranty --
  18. -- of MERCHANTABILITY  or  FITNESS FOR  A PARTICULAR PURPOSE.  See the  GNU --
  19. -- Library  General  Public  License for  more  details.  You  should  have --
  20. -- received  a copy of the GNU  Library  General Public License  along with --
  21. -- the GNAT library;  see the file  COPYING.LIB.  If not, write to the Free --
  22. -- Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.        --
  23. --                                                                          --
  24. ------------------------------------------------------------------------------
  25.  
  26. with Ada.Calendar;
  27. package body Ada.Numerics.Discrete_Random is
  28.  
  29.    -----------------------
  30.    -- Local Subprograms --
  31.    -----------------------
  32.  
  33.    function Square_Mod_N (X, N : Int_32) return Int_32;
  34.    pragma Inline (Square_Mod_N);
  35.  
  36.    procedure Test_Math (N : Int_32);
  37.  
  38.    ------------
  39.    -- Create --
  40.    ------------
  41.  
  42.    function Create return Pointer is
  43.    begin
  44.       return new State' (X1 => 2999 ** 2, X2 => 1439 ** 2,
  45.         P => 94_833_359, Q => 47_416_679, FP => 94_833_359.0,
  46.         Offset =>  Floating (Result_Subtype'Pos (Result_Subtype'First)) + 0.5,
  47.         Scale  => (Floating (Result_Subtype'Pos (Result_Subtype'Last)) -
  48.                    Floating (Result_Subtype'Pos (Result_Subtype'First)) + 1.0)
  49.                    / (94_833_359.0 * 47_416_679.0));
  50.    end Create;
  51.  
  52.    -----------
  53.    -- Image --
  54.    -----------
  55.  
  56.    function Image (Of_State : State) return String is
  57.    begin
  58.       return Int_32'Image (Of_State.X1) & ',' & Int_32'Image (Of_State.X2)
  59.              & ',' &
  60.              Int_32'Image (Of_State.P)  & ',' & Int_32'Image (Of_State.Q);
  61.    end Image;
  62.  
  63.    ------------
  64.    -- Random --
  65.    ------------
  66.  
  67.    function Random (Gen : Generator) return Result_Subtype is
  68.       Temp : Int_32;
  69.    begin
  70.       Gen.Point.X1 := Square_Mod_N (Gen.Point.X1, Gen.Point.P);
  71.       Gen.Point.X2 := Square_Mod_N (Gen.Point.X2, Gen.Point.Q);
  72.       Temp := Gen.Point.X2 - Gen.Point.X1;
  73.  
  74.       if Temp < 0 then
  75.          Temp := Temp + Gen.Point.Q;
  76.       end if;
  77.  
  78.       if Temp < 0 then
  79.          Temp := Temp + Gen.Point.Q;
  80.       end if;
  81.  
  82.       --  duplication is not an error, it is a loop unwinding.
  83.  
  84.       Temp := Int_32 (Gen.Point.Offset + (Floating (Temp) *
  85.          Floating (Gen.Point.P) + Floating (Gen.Point.X1)) * Gen.Point.Scale);
  86.  
  87.       if Temp > Int_32 (Result_Subtype'Pos (Result_Subtype'Last)) then
  88.          return Result_Subtype'First;
  89.       else
  90.          return Result_Subtype'Val (Temp);
  91.       end if;
  92.  
  93.       --  Pathological, but there do exist cases where the rounding implicit
  94.       --  in calculating the scale factor will cause rounding to 'Last + 1.
  95.       --  In those cases, returning 'First results in the least bias.
  96.    end Random;
  97.  
  98.    -----------
  99.    -- Reset --
  100.    -----------
  101.  
  102.    procedure Reset (Gen : Generator; Initiator : Integer) is
  103.       X1, X2, P, Q : Int_32;
  104.    begin
  105.       P := 94_833_359;
  106.       Q := 47_416_679;
  107.       X1 := 2 + Int_32 (Initiator) rem (P - 3);
  108.       X2 := 2 + Int_32 (Initiator) rem (Q - 3);
  109.  
  110.       for I in 1 .. 5 loop
  111.          X1 := Square_Mod_N (X1, P);
  112.          X2 := Square_Mod_N (X2, Q);
  113.       end loop;
  114.  
  115.       --  eliminate effects of small Initiators.
  116.  
  117.       Gen.Point.all :=
  118.         (X1 => X1, X2 => X2, P => P, Q => Q, FP => Floating (P),
  119.         Offset =>  Floating (Result_Subtype'Pos (Result_Subtype'First)) + 0.5,
  120.         Scale  => (Floating (Result_Subtype'Pos (Result_Subtype'Last)) -
  121.                    Floating (Result_Subtype'Pos (Result_Subtype'First)) + 1.0)
  122.                   / (Floating (P) * Floating (Q)));
  123.    end Reset;
  124.  
  125.    -----------
  126.    -- Reset --
  127.    -----------
  128.  
  129.    procedure Reset (Gen : Generator) is
  130.       use Ada.Calendar;
  131.       Now : Time := Clock;
  132.       X1, X2, P, Q : Int_32;
  133.    begin
  134.       P := 94_833_359;
  135.       Q := 47_416_679;
  136.       X1 := Int_32 (Year (Now)) * 12 * 31 +
  137.             Int_32 (Month (Now) * 31) + Int_32 (Day (Now));
  138.       X2 := Int_32  (Seconds (Now) * Duration (1000.0));
  139.       X1 := 2 + X1 rem (P - 3);
  140.       X2 := 2 + X2 rem (Q - 3);
  141.  
  142.       for I in 1 .. 5 loop
  143.          X1 := Square_Mod_N (X1, P);
  144.          X2 := Square_Mod_N (X2, Q);
  145.       end loop;
  146.  
  147.       --  less justification here, but eliminate visible effects of same day
  148.       --  starts.
  149.  
  150.       Gen.Point.all :=
  151.         (X1 => X1, X2 => X2, P => P, Q => Q, FP => Floating (P),
  152.          Offset =>  Floating (Result_Subtype'Pos (Result_Subtype'First)) + 0.5,
  153.          Scale  => (Floating (Result_Subtype'Pos (Result_Subtype'Last)) -
  154.                     Floating (Result_Subtype'Pos (Result_Subtype'First)) + 1.0)
  155.                   / (Floating (P) * Floating (Q)));
  156.  
  157.       --  Why save the actual assignments to the end?  To insure to the
  158.       --  greatest extent possible that an exception won't leave the generator
  159.       --  inconsistent.
  160.  
  161.    end Reset;
  162.  
  163.    -----------
  164.    -- Reset --
  165.    -----------
  166.  
  167.    procedure Reset (Gen : Generator; From_State : State) is
  168.    begin
  169.       Gen.Point.all := From_State;
  170.    end Reset;
  171.  
  172.    ----------
  173.    -- Save --
  174.    ----------
  175.  
  176.    procedure Save (Gen : Generator; To_State : out State) is
  177.    begin
  178.       To_State := Gen.Point.all;
  179.    end Save;
  180.  
  181.    ------------------
  182.    -- Square_Mod_N --
  183.    ------------------
  184.  
  185.    function Square_Mod_N (X, N : Int_32) return Int_32 is
  186.       Temp : Floating := Floating (X) * Floating (X);
  187.       Div  : Int_32   := Int_32 (Temp / Floating (N));
  188.    begin
  189.       Div := Int_32 (Temp - Floating (Div) * Floating (N));
  190.       if Div < 0 then
  191.          return Div + N;
  192.       else
  193.          return Div;
  194.       end if;
  195.    end Square_Mod_N;
  196.  
  197.  
  198.    ---------------
  199.    -- Test_Math --
  200.    ---------------
  201.  
  202.    procedure Test_Math (N : Int_32) is
  203.    begin
  204.       if Square_Mod_N (N - 1, N) /= 1 then
  205.          raise Program_error;
  206.       end if;
  207.    end Test_Math;
  208.  
  209.    -----------
  210.    -- Value --
  211.    -----------
  212.  
  213.    function Value (Coded_State : String) return State is
  214.       Start, Stop : Positive := Coded_State'First;
  215.       Out_State : State;
  216.  
  217.    begin
  218.       while Coded_State (Stop) /= ',' loop
  219.          Stop := Stop + 1;
  220.       end loop;
  221.  
  222.       Out_State.X1 := Int_32'Value (Coded_State (Start .. Stop - 1));
  223.       Start := Stop + 1;
  224.  
  225.       loop
  226.          Stop := Stop + 1;
  227.          exit when Coded_State (Stop) = ',';
  228.       end loop;
  229.  
  230.       Out_State.X2 := Int_32'Value (Coded_State (Start .. Stop - 1));
  231.       Out_State.Q  := Int_32'Value (Coded_State
  232.                                        (Stop + 1 .. Coded_State'Last));
  233.       Out_State.P  := Out_State.Q * 2 + 1;
  234.       Out_State.FP := Floating (Out_State.P);
  235.       Out_State.Scale :=
  236.           (Floating (Result_Subtype'Pos (Result_Subtype'Last)) -
  237.            Floating (Result_Subtype'Pos (Result_Subtype'First)) + 1.0)
  238.        /  (Floating (Out_State.P) * Floating (Out_State.Q));
  239.  
  240.       --  now do SOME sanity checks.
  241.  
  242.       if Out_State.Q < 31
  243.         or else Out_State.X1 not in 2 .. Out_State.P - 1
  244.         or else Out_State.X2 not in 2 .. Out_State.Q - 1
  245.       then
  246.          raise Constraint_Error;
  247.       end if;
  248.  
  249.       Test_Math (Out_State.P);
  250.       return Out_State;
  251.    end Value;
  252.  
  253. begin
  254.    Test_Math (94_833_359);
  255. end Ada.Numerics.Discrete_Random;
  256.