with Ada.Numerics;                      use Ada.Numerics;
with Ada.Numerics.Generic_Elementary_Functions;

with Beta_function, Phi_function;

package body Generic_Random_Functions is

  -- Ada 95 Quality and Style Guide, 7.2.7:
  -- Tests for
  --
  -- (1) absolute "equality" to 0 in storage,
  -- (2) absolute "equality" to 0 in computation,
  -- (3) relative "equality" to 0 in storage, and
  -- (4) relative "equality" to 0 in computation:
  --
  --  abs X <= Float_Type'Model_Small                      -- (1)
  --  abs X <= Float_Type'Base'Model_Small                 -- (2)
  --  abs X <= abs X * Float_Type'Model_Epsilon            -- (3)
  --  abs X <= abs X * Float_Type'Base'Model_Epsilon       -- (4)

  function Almost_zero(X: Real) return Boolean is
  begin
    return  abs X <= Real'Base'Model_Small;
  end Almost_zero;

  package RPhF is new Phi_function (Real);
  function Normal_CDF (x: Real) return Real renames RPhF.Phi;
  function Normal_inverse_CDF (y : Real) return Real renames RPhF.Inverse_Phi;

  package RBF is new Beta_function (Real);
  function Beta_CDF (x, a, b: Real) return Real renames RBF.Regularized_Beta;
  function Beta_inverse_CDF (x, a, b: Real) return Real renames RBF.Inverse_Regularized_Beta;

  package GEF is new Ada.Numerics.Generic_Elementary_Functions(Real);
  use GEF;

  procedure Box_Muller(u1,u2: in Real; n1,n2: out Real) is
    phi, r: Real;
    Two_Pi  : constant:= 2.0 * Pi;
    m2 : constant:= -2.0;
  begin
    r:= Sqrt(m2 * Log(u2));
    phi:= Two_Pi * u1;
    n1:= r * Cos(phi);
    n2:= r * Sin(phi);
  end Box_Muller;

  function Poisson(lambda: Real) return Natural is
    lower_limit, product: Real;
    k: Integer;
  begin
    -- Algo found in:
    -- http://en.wikipedia.org/wiki/Poisson_distribution
    --   #Generating_Poisson-distributed_random_variables
    -- referring to:
    -- Knuth (whom else ?!), The Art of Computer Programming,
    -- Volume 2, Seminumerical algorithms, 3.4.1. Numerical Distributions,
    -- F. Important integer-valued distributions.
    lower_limit:= Exp(-lambda);
    k:= -1;
    product:= 1.0;
    loop
      k:= k + 1;
      product:= product * U;
      exit when product <= lower_limit;
    end loop;
    return k;
  end Poisson;

  function Pareto_inverse_CDF(q, threshold, minus_inv_alpha: Real) return Real is
  begin
    if Almost_zero(q) then
      return Real'Last;
    end if;
    return threshold * q ** (minus_inv_alpha);
  end Pareto_inverse_CDF;

  function Pareto_CDF(x, threshold, alpha: Real) return Real is
  begin
    return 1.0 - (threshold / x) ** alpha;
  end Pareto_CDF;

end Generic_Random_Functions;