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-nuflra.adb < prev    next >
Text File  |  1996-09-28  |  9KB  |  299 lines
  1. ------------------------------------------------------------------------------
  2. --                                                                          --
  3. --                         GNAT RUNTIME COMPONENTS                          --
  4. --                                                                          --
  5. --            A D A . N U M E R I C S . F L O A T _ R A N D O M             --
  6. --                                                                          --
  7. --                                 B o d y                                  --
  8. --                                                                          --
  9. --                            $Revision: 1.11 $                             --
  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.Float_Random is
  28.  
  29.    -----------------------
  30.    -- Local Subprograms --
  31.    -----------------------
  32.  
  33.    procedure Euclid (P, Q : in Int_32; X, Y : out Int_32; GCD : out Int_32);
  34.    function  Euclid (P, Q : Int_32) return Int_32;
  35.    function Square_Mod_N (X, N : Int_32) return Int_32;
  36.  
  37.    pragma Inline (Square_Mod_N);
  38.    procedure Test_Math (N : in Int_32);
  39.  
  40.    ------------
  41.    -- Create --
  42.    ------------
  43.  
  44.    function Create return Pointer is
  45.    begin
  46.       return new State' (Initial_State);
  47.    end Create;
  48.  
  49.    ------------
  50.    -- Euclid --
  51.    ------------
  52.  
  53.    procedure Euclid (P, Q : in Int_32; X, Y : out Int_32; GCD : out Int_32) is
  54.  
  55.       XT : Int_32 := 1;
  56.       YT : Int_32 := 0;
  57.  
  58.       procedure Recur
  59.         (P,  Q  : in Int_32;                 --  a (i-1), a (i)
  60.          X,  Y  : in Int_32;                 --  x (i),   y (i)
  61.          XP, YP : in out Int_32;             --  x (i-1), y (i-1)
  62.          GCD    : out Int_32);
  63.  
  64.       procedure Recur
  65.         (P,  Q  : in Int_32;
  66.          X,  Y  : in Int_32;
  67.          XP, YP : in out Int_32;
  68.          GCD    : out Int_32)
  69.       is
  70.          Quo : Int_32  := P / Q;             --  q <-- |_ a (i-1) / a (i) _|
  71.          XT  : Int_32 := X;                  --  x (i)
  72.          YT  : Int_32 := Y;                  --  y (i)
  73.  
  74.       begin
  75.          if P rem Q = 0 then                 --  while does not divide
  76.             GCD := Q;
  77.             XP  := X;
  78.             YP  := Y;
  79.          else
  80.             Recur (Q, P - Q * Quo, XP - Quo * X, YP - Quo * Y, XT, YT, Quo);
  81.             --  a (i) <== a (i)
  82.             --  a (i+1) <-- a (i-1) - q*a (i)
  83.             --  x (i+1) <-- x (i-1) - q*x (i)
  84.             --  y (i+1) <-- y (i-1) - q*y (i)
  85.             --  x (i) <== x (i)
  86.             --  y (i) <== y (i)
  87.             XP  := XT;
  88.             YP  := YT;
  89.             GCD := Quo;
  90.          end if;
  91.       end Recur;
  92.  
  93.    begin
  94.       Recur (P, Q, 0, 1, XT, YT, GCD);
  95.       X := XT;
  96.       Y := YT;
  97.    end Euclid;
  98.  
  99.    ------------
  100.    -- Euclid --
  101.    ------------
  102.  
  103.    function Euclid (P, Q : Int_32) return Int_32 is
  104.       X, Y, GCD : Int_32;
  105.    begin
  106.       Euclid (P, Q, X, Y, GCD);
  107.       return X;
  108.    end Euclid;
  109.  
  110.    -----------
  111.    -- Image --
  112.    -----------
  113.  
  114.    function Image (Of_State : State) return String is
  115.    begin
  116.       return Int_32'Image (Of_State.X1) & ',' & Int_32'Image (Of_State.X2)
  117.              & ',' &
  118.              Int_32'Image (Of_State.P)  & ',' & Int_32'Image (Of_State.Q)
  119.              & ',' &
  120.              Int_32'Image (Of_State.X);
  121.    end Image;
  122.  
  123.    ------------
  124.    -- Random --
  125.    ------------
  126.  
  127.    function Random  (Gen : Generator) return  Uniformly_Distributed is
  128.       X1, X2, Temp : Int_32;
  129.  
  130.    begin
  131.       Gen.Point.X1 := Square_Mod_N (Gen.Point.X1,  Gen.Point.P);
  132.       Gen.Point.X2 := Square_Mod_N (Gen.Point.X2,  Gen.Point.Q);
  133.       return
  134.         Float ((Longer_Float (((Gen.Point.X2 - Gen.Point.X1) * Gen.Point.X)
  135.         mod Gen.Point.Q) * Longer_Float (Gen.Point.P)
  136.         + Longer_Float (Gen.Point.X1)) * Gen.Point.Scale);
  137.    end Random;
  138.  
  139.    -----------
  140.    -- Reset --
  141.    -----------
  142.  
  143.    procedure Reset (Gen : in Generator; Initiator : in Integer) is
  144.       X1, X2, P, Q : Int_32;
  145.  
  146.    begin
  147.       P := Initial_State.P;
  148.       Q := Initial_State.Q;
  149.       X1 := 2 + Int_32 (Initiator) rem  (P - 3);
  150.       X2 := 2 + Int_32 (Initiator) rem  (Q - 3);
  151.       for I in 1 .. 5 loop
  152.          X1 := Square_Mod_N (X1, P);
  153.          X2 := Square_Mod_N (X2, Q);
  154.       end loop;
  155.       --  eliminate effects of small Initiators.
  156.       Gen.Point.all :=
  157.         (X1 => X1, X2 => X2, P => P, Q => Q, X => Euclid (P, Q),
  158.          Scale => 1.0 / (Longer_Float (P) * Longer_Float (Q)));
  159.    end Reset;
  160.  
  161.    -----------
  162.    -- Reset --
  163.    -----------
  164.  
  165.    procedure Reset (Gen : Generator; From_State : State) is
  166.    begin
  167.       Gen.Point.all := From_State;
  168.    end Reset;
  169.  
  170.    -----------
  171.    -- Reset --
  172.    -----------
  173.  
  174.    procedure Reset (Gen : Generator) is
  175.       use Ada.Calendar;
  176.       Now : Time := Clock;
  177.       X1, X2, P, Q : Int_32;
  178.    begin
  179.       P := Initial_State.P;
  180.       Q := Initial_State.Q;
  181.       X1 := Int_32 (Year (Now)) * 12 * 31 +
  182.             Int_32 (Month (Now)) * 31 +
  183.             Int_32 (Day (Now));
  184.       X2 := Int_32 (Seconds (Now) * Duration (1000.0));
  185.       --  X2 := Int_32 (1000 * Seconds (Now)); raises Constraint_Error
  186.       X1 := 2 + X1 rem  (P - 3);
  187.       X2 := 2 + X2 rem  (Q - 3);
  188.       for I in 1 .. 5 loop
  189.          X1 := Square_Mod_N (X1, P);
  190.          X2 := Square_Mod_N (X2, Q);
  191.       end loop;
  192.       --  less justification here, but eliminate visible effects of same day
  193.       --  starts.
  194.       Gen.Point.all :=
  195.         (X1 => X1, X2 => X2, P => P, Q => Q,
  196.          X => Euclid (P, Q),
  197.          Scale => 1.0 / (Longer_Float (P) * Longer_Float (Q)));
  198.  
  199.       --  Why save the actual assignments to the end?  To insure to the
  200.       --  greatest extent possible that an exception won't leave the generator
  201.       --  inconsistant.
  202.  
  203.    end Reset;
  204.  
  205.    ----------
  206.    -- Save --
  207.    ----------
  208.  
  209.    procedure Save (Gen : in Generator; To_State : out State) is
  210.    begin
  211.       To_State := Gen.Point.all;
  212.    end Save;
  213.  
  214.    ------------------
  215.    -- Square_Mod_N --
  216.    ------------------
  217.  
  218.    --  don't mess with working code, but I want a function that returns X.
  219.    --  Note rem is used below instead of mod where both arguments are known
  220.    --  positive. This is faster on some hardware.
  221.  
  222.    function Square_Mod_N (X, N : Int_32) return Int_32 is
  223.       subtype LF is Longer_Float;
  224.       Temp : LF  := LF (X) * LF (X);
  225.       Div : Int_32 := Int_32 (Temp / LF (N));
  226.    begin
  227.       Div := Int_32 (Temp - LF (Div) * LF (N));
  228.  
  229.       if Div < 0 then
  230.          return Div + N;
  231.       else
  232.          return Div;
  233.       end if;
  234.    end Square_Mod_N;
  235.  
  236.    ---------------
  237.    -- Test_Math --
  238.    ---------------
  239.  
  240.    procedure Test_Math (N : in Int_32) is
  241.    begin
  242.       if Square_Mod_N (N - 1, N) /= 1 then
  243.          raise Program_Error;
  244.       end if;
  245.    end Test_Math;
  246.  
  247.    -----------
  248.    -- Value --
  249.    -----------
  250.  
  251.    function Value (Coded_State : String) return State is
  252.       Start, Stop : Positive := Coded_State'First;
  253.       Out_State : State;
  254.  
  255.    begin
  256.       while Coded_State (Stop) /= ',' loop
  257.          Stop := Stop + 1;
  258.       end loop;
  259.  
  260.       Out_State.X1 := Int_32'Value (Coded_State (Start .. Stop - 1));
  261.       Start := Stop + 1;
  262.  
  263.       loop
  264.          Stop := Stop + 1;
  265.          exit when Coded_State (Stop) = ',';
  266.       end loop;
  267.  
  268.       Out_State.X2 := Int_32'Value (Coded_State (Start .. Stop - 1));
  269.       Start := Stop + 1;
  270.  
  271.       loop
  272.          Stop := Stop + 1;
  273.          exit when Coded_State (Stop) = ',';
  274.       end loop;
  275.  
  276.       Out_State.P := Int_32'Value (Coded_State (Start .. Stop - 1));
  277.       Out_State.Q := Int_32'Value (Coded_State (Stop + 1 .. Coded_State'Last));
  278.       Out_State.X := Euclid (Out_State.P, Out_State.Q);
  279.       Out_State.Scale := 1.0
  280.         / (Longer_Float (Out_State.P) * Longer_Float (Out_State.Q));
  281.  
  282.       --  now do SOME sanity checks.
  283.  
  284.       if Out_State.Q < 31 or else Out_State.P < 31
  285.         or else Out_State.X1 not in 2 .. Out_State.P - 1
  286.         or else Out_State.X2 not in 2 .. Out_State.Q - 1
  287.       then
  288.          raise Constraint_Error;
  289.       end if;
  290.  
  291.       Test_Math (Out_State.P);
  292.       Test_Math (Out_State.Q);
  293.       return Out_State;
  294.    end Value;
  295.  
  296. begin
  297.    Test_Math (94_833_359);
  298. end Ada.Numerics.Float_Random;
  299.