with Globals;
use Globals;

pragma Elaborate_All(Globals);

package body Flts is

  function RMin(R1,R2,R3: Flt) return Flt is
  begin
    return RMin(RMin(R1,R2),R3);
  end RMin;

  function RMax(R1,R2,R3: Flt) return Flt is
  begin
    return RMax(RMax(R1,R2),R3);
  end RMax;

  function Sqr(R: Flt) return Flt is
  begin
    return R*R;
  end Sqr;

  function NPower(R: Flt; N: Natural) return Flt is
    I: Integer := N;
    X: Flt := R;
    Y: Flt := One;
  begin
    if N /= 0 then
      loop
        if (I rem 2)=1 then
          Y := X*Y;
        end if;
        I := I/2;
        exit when I=0;
        X := X*X;
      end loop;
    end if;
    return Y;
  end NPower;

  function Convert(K: Integer) return Flt is
  begin
    return Flt(K);
  end Convert;

  function Max(P: RadVec) return Radius is
    R: Radius := P(P'First);
  begin
    for I in P'First+1 .. P'Last loop
      R := RMax(P(I),R);
    end loop;
    return R;
  end Max;

  procedure Swap(P1,P2: in out RadPair) is
    P: constant RadPair := P1;
  begin
    P1 := P2;
    P2 := P;
  end Swap;

  function Frac(N: Integer; D: Positive) return Flt is
  begin
    return Flt(N)/Flt(D);
  end Frac;

  function Between(A,B: Flt; P: Natural; Q: Positive) return Flt is
  begin
    if P=0 then return A; end if;
    if P=Q then return B; end if;
    if P>Q then raise Domain_Error; end if;
    return (Flt(Q-P)*A+Flt(P)*B)/Flt(Q);
  end Between;

  function Scaled(Eps: in Flt; Binary_Digits: in Positive) return Flt is
    use Flt_EF;
  begin
    if Binary_Digits=64 then
      return Eps;
    else
      return Exp(Log(Eps)*Flt(Binary_Digits)/Flt(64));
    end if;
  end Scaled;

  function Epsilon(Eps: in Flt; Binary_Digits: in Positive) return Flt is
    Log2:   constant Flt := 0.65625; -- less
    use Flt_EF;
  begin
    if Binary_Digits <= 64 then
      return Eps;
    else
      return Eps*Exp(-Log2*Flt(Binary_Digits-64));
    end if;
  end Epsilon;

  function Simple_Random return Flt is
  begin
    return Two*Flt(Random(Gen))-One;
  end Simple_Random;

  procedure Reset_Random(I: Integer) is
  begin
    Reset(Gen,I);
  end Reset_Random;

  ------------------------------

  function Min(P1,P2: RadPair) return RadPair is
  begin
    return (RMin(P1(1),P2(1)),RMin(P1(2),P2(2)));
  end Min;

  function Max(P1,P2: RadPair) return RadPair is
  begin
    return (RMax(P1(1),P2(1)),RMax(P1(2),P2(2)));
  end Max;

  function Between(P1,P2: RadPair; R: Flt := Half) return RadPair is
    --- no claim of accuracy
    P: RadPair;
  begin
    P(1) := (One-R)*P1(1)+R*P2(1);
    P(2) := (One-R)*P1(2)+R*P2(2);
    return P;
  end Between;

  function PowersOf(R: Radius; MaxPower: Positive) return RadVec is
    V: RadVec(0 .. MaxPower);
  begin
    V(0) := One;
    V(1) := R;
    for N in 2 .. MaxPower loop
      V(N) := V(N/2)*V(N-N/2);
    end loop;
    return V;
  end PowersOf;

end Flts;
