--  Copyright (2008) Cédric Coussinet (cedric.coussinet@nomoseed.net)
--
--  This program is free software: you can redistribute it and/or modify
--  it under the terms of the GNU Affero General Public License as published
--  by the Free Software Foundation, either version 3 of the License, or
--  (at your option) any later version.
--
--  This program 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 Affero General Public License for more details.

--  You should have received a copy of the GNU Affero General Public License
--  along with this program. If not, see <http://www.gnu.org/licenses/>

with Ada.Numerics.Elementary_Functions;

package body Standard_Random is

   type Unsigned_32 is mod 2 ** 32;
   for Unsigned_32'Size use 32;

   A : Unsigned_32 := 1;
   B : constant Unsigned_32 := 16807;
   C : constant Unsigned_32 := Unsigned_32 (Rand_Max + 1);

   function Random return Positive is
   begin
      A := (A * B) mod C;
      return Positive (A);
   end Random;

   procedure Reset is
   begin
      A := 1;
   end Reset;

   procedure Reset (Initiator : in Positive) is
   begin
      A := Unsigned_32 (Initiator);
   end Reset;

   V1Fac : Float;
   Have_V1Fac : Boolean := False;

   function Random_Gaussian (Average : in Float;
                             Sigma   : in Float) return Float is
      use Ada.Numerics.Elementary_Functions;
      Fac, Rsq, V1, V2 : Float;
   begin
      if Have_V1Fac then
         Have_V1Fac := False;
         return V1Fac * Sigma  + Average;
      else
         loop
            V1 := 2.0 * Float (Random) / Float (Rand_Max) - 1.0;
            V2 := 2.0 * Float (Random) / Float (Rand_Max) - 1.0;
            Rsq := V1 ** 2 + V2 ** 2;
            exit when Rsq < 1.0 and then Rsq /= 0.0;
         end loop;
         Fac := Sqrt (-2.0 * Log (Rsq) / Rsq);
         --  Transformation de Box-Muller
         --  un des deux nombres est gardé pour la prochaine fois
         V1Fac := V1 * Fac;
         Have_V1Fac := True;
         return (V2 * Fac * Sigma + Average);
      end if;
   end Random_Gaussian;

end Standard_Random;
