------------------------------------------------------------------------------
--                                                                          --
--                         GNAT RUNTIME COMPONENTS                          --
--                                                                          --
--         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          --
--                                                                          --
--                                 B o d y                                  --
--                                                                          --
--                            $Revision: 1.9 $                              --
--                                                                          --
--           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.Discrete_Random is

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

   function Square_Mod_N (X, N : Int_32) return Int_32;
   pragma Inline (Square_Mod_N);

   procedure Test_Math (N : Int_32);

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

   function Create return Pointer is
   begin
      return new State' (X1 => 2999 ** 2, X2 => 1439 ** 2,
        P => 94_833_359, Q => 47_416_679, FP => 94_833_359.0,
        Offset =>  Floating (Result_Subtype'Pos (Result_Subtype'First)) + 0.5,
        Scale  => (Floating (Result_Subtype'Pos (Result_Subtype'Last)) -
                   Floating (Result_Subtype'Pos (Result_Subtype'First)) + 1.0)
                   / (94_833_359.0 * 47_416_679.0));
   end Create;

   -----------
   -- 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);
   end Image;

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

   function Random (Gen : Generator) return Result_Subtype is
      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);
      Temp := Gen.Point.X2 - Gen.Point.X1;

      if Temp < 0 then
         Temp := Temp + Gen.Point.Q;
      end if;

      if Temp < 0 then
         Temp := Temp + Gen.Point.Q;
      end if;

      --  duplication is not an error, it is a loop unwinding.

      Temp := Int_32 (Gen.Point.Offset + (Floating (Temp) *
         Floating (Gen.Point.P) + Floating (Gen.Point.X1)) * Gen.Point.Scale);

      if Temp > Int_32 (Result_Subtype'Pos (Result_Subtype'Last)) then
         return Result_Subtype'First;
      else
         return Result_Subtype'Val (Temp);
      end if;

      --  Pathological, but there do exist cases where the rounding implicit
      --  in calculating the scale factor will cause rounding to 'Last + 1.
      --  In those cases, returning 'First results in the least bias.
   end Random;

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

   procedure Reset (Gen : Generator; Initiator : Integer) is
      X1, X2, P, Q : Int_32;
   begin
      P := 94_833_359;
      Q := 47_416_679;
      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, FP => Floating (P),
        Offset =>  Floating (Result_Subtype'Pos (Result_Subtype'First)) + 0.5,
        Scale  => (Floating (Result_Subtype'Pos (Result_Subtype'Last)) -
                   Floating (Result_Subtype'Pos (Result_Subtype'First)) + 1.0)
                  / (Floating (P) * Floating (Q)));
   end Reset;

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

   procedure Reset (Gen : Generator) is
      use Ada.Calendar;
      Now : Time := Clock;
      X1, X2, P, Q : Int_32;
   begin
      P := 94_833_359;
      Q := 47_416_679;
      X1 := Int_32 (Year (Now)) * 12 * 31 +
            Int_32 (Month (Now) * 31) + Int_32 (Day (Now));
      X2 := Int_32  (Seconds (Now) * Duration (1000.0));
      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, FP => Floating (P),
         Offset =>  Floating (Result_Subtype'Pos (Result_Subtype'First)) + 0.5,
         Scale  => (Floating (Result_Subtype'Pos (Result_Subtype'Last)) -
                    Floating (Result_Subtype'Pos (Result_Subtype'First)) + 1.0)
                  / (Floating (P) * Floating (Q)));

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

   end Reset;

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

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

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

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

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

   function Square_Mod_N (X, N : Int_32) return Int_32 is
      Temp : Floating := Floating (X) * Floating (X);
      Div  : Int_32   := Int_32 (Temp / Floating (N));
   begin
      Div := Int_32 (Temp - Floating (Div) * Floating (N));
      if Div < 0 then
         return Div + N;
      else
         return Div;
      end if;
   end Square_Mod_N;


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

   procedure Test_Math (N : 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));
      Out_State.Q  := Int_32'Value (Coded_State
                                       (Stop + 1 .. Coded_State'Last));
      Out_State.P  := Out_State.Q * 2 + 1;
      Out_State.FP := Floating (Out_State.P);
      Out_State.Scale :=
          (Floating (Result_Subtype'Pos (Result_Subtype'Last)) -
           Floating (Result_Subtype'Pos (Result_Subtype'First)) + 1.0)
       /  (Floating (Out_State.P) * Floating (Out_State.Q));

      --  now do SOME sanity checks.

      if Out_State.Q < 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);
      return Out_State;
   end Value;

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