------------------------------------------------------------------------------
--                                                                          --
--                         GNAT RUNTIME COMPONENTS                          --
--                                                                          --
--            A D A . N U M E R I C S . F L O A T _ R A N D O M             --
--                                                                          --
--                                 B o d y                                  --
--                                                                          --
--                            $Revision: 1.11 $                             --
--                                                                          --
--           Copyright (c) 1992,1993,1994 NYU, All Rights Reserved          --
--                                                                          --
-- The GNAT library is free software; you can redistribute it and/or modify --
-- it under terms of the GNU Library General Public License as published by --
-- the Free Software  Foundation; either version 2, or (at your option) any --
-- later version.  The GNAT library is distributed in the hope that it will --
-- be useful, but WITHOUT ANY WARRANTY;  without even  the implied warranty --
-- of MERCHANTABILITY  or  FITNESS FOR  A PARTICULAR PURPOSE.  See the  GNU --
-- Library  General  Public  License for  more  details.  You  should  have --
-- received  a copy of the GNU  Library  General Public License  along with --
-- the GNAT library;  see the file  COPYING.LIB.  If not, write to the Free --
-- Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.        --
--                                                                          --
------------------------------------------------------------------------------

with Ada.Calendar;
package body Ada.Numerics.Float_Random is

   -----------------------
   -- Local Subprograms --
   -----------------------

   procedure Euclid (P, Q : in Int_32; X, Y : out Int_32; GCD : out Int_32);
   function  Euclid (P, Q : Int_32) return Int_32;
   function Square_Mod_N (X, N : Int_32) return Int_32;

   pragma Inline (Square_Mod_N);
   procedure Test_Math (N : in Int_32);

   ------------
   -- Create --
   ------------

   function Create return Pointer is
   begin
      return new State' (Initial_State);
   end Create;

   ------------
   -- Euclid --
   ------------

   procedure Euclid (P, Q : in Int_32; X, Y : out Int_32; GCD : out Int_32) is

      XT : Int_32 := 1;
      YT : Int_32 := 0;

      procedure Recur
        (P,  Q  : in Int_32;                 --  a (i-1), a (i)
         X,  Y  : in Int_32;                 --  x (i),   y (i)
         XP, YP : in out Int_32;             --  x (i-1), y (i-1)
         GCD    : out Int_32);

      procedure Recur
        (P,  Q  : in Int_32;
         X,  Y  : in Int_32;
         XP, YP : in out Int_32;
         GCD    : out Int_32)
      is
         Quo : Int_32  := P / Q;             --  q <-- |_ a (i-1) / a (i) _|
         XT  : Int_32 := X;                  --  x (i)
         YT  : Int_32 := Y;                  --  y (i)

      begin
         if P rem Q = 0 then                 --  while does not divide
            GCD := Q;
            XP  := X;
            YP  := Y;
         else
            Recur (Q, P - Q * Quo, XP - Quo * X, YP - Quo * Y, XT, YT, Quo);
            --  a (i) <== a (i)
            --  a (i+1) <-- a (i-1) - q*a (i)
            --  x (i+1) <-- x (i-1) - q*x (i)
            --  y (i+1) <-- y (i-1) - q*y (i)
            --  x (i) <== x (i)
            --  y (i) <== y (i)
            XP  := XT;
            YP  := YT;
            GCD := Quo;
         end if;
      end Recur;

   begin
      Recur (P, Q, 0, 1, XT, YT, GCD);
      X := XT;
      Y := YT;
   end Euclid;

   ------------
   -- Euclid --
   ------------

   function Euclid (P, Q : Int_32) return Int_32 is
      X, Y, GCD : Int_32;
   begin
      Euclid (P, Q, X, Y, GCD);
      return X;
   end Euclid;

   -----------
   -- Image --
   -----------

   function Image (Of_State : State) return String is
   begin
      return Int_32'Image (Of_State.X1) & ',' & Int_32'Image (Of_State.X2)
             & ',' &
             Int_32'Image (Of_State.P)  & ',' & Int_32'Image (Of_State.Q)
             & ',' &
             Int_32'Image (Of_State.X);
   end Image;

   ------------
   -- Random --
   ------------

   function Random  (Gen : Generator) return  Uniformly_Distributed is
      X1, X2, Temp : Int_32;

   begin
      Gen.Point.X1 := Square_Mod_N (Gen.Point.X1,  Gen.Point.P);
      Gen.Point.X2 := Square_Mod_N (Gen.Point.X2,  Gen.Point.Q);
      return
        Float ((Longer_Float (((Gen.Point.X2 - Gen.Point.X1) * Gen.Point.X)
        mod Gen.Point.Q) * Longer_Float (Gen.Point.P)
        + Longer_Float (Gen.Point.X1)) * Gen.Point.Scale);
   end Random;

   -----------
   -- Reset --
   -----------

   procedure Reset (Gen : in Generator; Initiator : in Integer) is
      X1, X2, P, Q : Int_32;

   begin
      P := Initial_State.P;
      Q := Initial_State.Q;
      X1 := 2 + Int_32 (Initiator) rem  (P - 3);
      X2 := 2 + Int_32 (Initiator) rem  (Q - 3);
      for I in 1 .. 5 loop
         X1 := Square_Mod_N (X1, P);
         X2 := Square_Mod_N (X2, Q);
      end loop;
      --  eliminate effects of small Initiators.
      Gen.Point.all :=
        (X1 => X1, X2 => X2, P => P, Q => Q, X => Euclid (P, Q),
         Scale => 1.0 / (Longer_Float (P) * Longer_Float (Q)));
   end Reset;

   -----------
   -- Reset --
   -----------

   procedure Reset (Gen : Generator; From_State : State) is
   begin
      Gen.Point.all := From_State;
   end Reset;

   -----------
   -- Reset --
   -----------

   procedure Reset (Gen : Generator) is
      use Ada.Calendar;
      Now : Time := Clock;
      X1, X2, P, Q : Int_32;
   begin
      P := Initial_State.P;
      Q := Initial_State.Q;
      X1 := Int_32 (Year (Now)) * 12 * 31 +
            Int_32 (Month (Now)) * 31 +
            Int_32 (Day (Now));
      X2 := Int_32 (Seconds (Now) * Duration (1000.0));
      --  X2 := Int_32 (1000 * Seconds (Now)); raises Constraint_Error
      X1 := 2 + X1 rem  (P - 3);
      X2 := 2 + X2 rem  (Q - 3);
      for I in 1 .. 5 loop
         X1 := Square_Mod_N (X1, P);
         X2 := Square_Mod_N (X2, Q);
      end loop;
      --  less justification here, but eliminate visible effects of same day
      --  starts.
      Gen.Point.all :=
        (X1 => X1, X2 => X2, P => P, Q => Q,
         X => Euclid (P, Q),
         Scale => 1.0 / (Longer_Float (P) * Longer_Float (Q)));

      --  Why save the actual assignments to the end?  To insure to the
      --  greatest extent possible that an exception won't leave the generator
      --  inconsistant.

   end Reset;

   ----------
   -- Save --
   ----------

   procedure Save (Gen : in Generator; To_State : out State) is
   begin
      To_State := Gen.Point.all;
   end Save;

   ------------------
   -- Square_Mod_N --
   ------------------

   --  don't mess with working code, but I want a function that returns X.
   --  Note rem is used below instead of mod where both arguments are known
   --  positive. This is faster on some hardware.

   function Square_Mod_N (X, N : Int_32) return Int_32 is
      subtype LF is Longer_Float;
      Temp : LF  := LF (X) * LF (X);
      Div : Int_32 := Int_32 (Temp / LF (N));
   begin
      Div := Int_32 (Temp - LF (Div) * LF (N));

      if Div < 0 then
         return Div + N;
      else
         return Div;
      end if;
   end Square_Mod_N;

   ---------------
   -- Test_Math --
   ---------------

   procedure Test_Math (N : in Int_32) is
   begin
      if Square_Mod_N (N - 1, N) /= 1 then
         raise Program_Error;
      end if;
   end Test_Math;

   -----------
   -- Value --
   -----------

   function Value (Coded_State : String) return State is
      Start, Stop : Positive := Coded_State'First;
      Out_State : State;

   begin
      while Coded_State (Stop) /= ',' loop
         Stop := Stop + 1;
      end loop;

      Out_State.X1 := Int_32'Value (Coded_State (Start .. Stop - 1));
      Start := Stop + 1;

      loop
         Stop := Stop + 1;
         exit when Coded_State (Stop) = ',';
      end loop;

      Out_State.X2 := Int_32'Value (Coded_State (Start .. Stop - 1));
      Start := Stop + 1;

      loop
         Stop := Stop + 1;
         exit when Coded_State (Stop) = ',';
      end loop;

      Out_State.P := Int_32'Value (Coded_State (Start .. Stop - 1));
      Out_State.Q := Int_32'Value (Coded_State (Stop + 1 .. Coded_State'Last));
      Out_State.X := Euclid (Out_State.P, Out_State.Q);
      Out_State.Scale := 1.0
        / (Longer_Float (Out_State.P) * Longer_Float (Out_State.Q));

      --  now do SOME sanity checks.

      if Out_State.Q < 31 or else Out_State.P < 31
        or else Out_State.X1 not in 2 .. Out_State.P - 1
        or else Out_State.X2 not in 2 .. Out_State.Q - 1
      then
         raise Constraint_Error;
      end if;

      Test_Math (Out_State.P);
      Test_Math (Out_State.Q);
      return Out_State;
   end Value;

begin
   Test_Math (94_833_359);
end Ada.Numerics.Float_Random;
