File : intervals-ops.adb

with Reps.Ops, Round_Up;

pragma Elaborate_All (Reps.Ops);

package body Intervals.Ops is

  NegOne:  constant Rep := -One;
  NegTwo:  constant Rep := -Two;
  NegHalf: constant Rep := -Half;
  Three:   constant Rep := Rep(3);
  Quarter: constant Rep := Half/Two;
  Epsilon: constant Rep := Half**99;

  Interval_Zero: constant Interval := (Zero,Zero);
  Interval_One:  constant Interval := (One,Zero);
  PlusMinus:     constant Interval := (Zero,One);

  function Scal(R1,R2: Rep) return Interval is
    Sum: constant Rep := R1+R2;
  begin
    if R1 >= R2 then
      return (Half*Sum,Half*((R1-R2)+(Sum+(-R1-R2))));
    end if;
    return (Half*Sum,Half*((R2-R1)+(Sum+(-R1-R2))));
  end Scal;

  function Ball(S: Interval) return Interval is
  begin
    return (Zero,Abs(S.C)+S.R);
  end Ball;

  function Inf(S: Interval) return Rep is
  begin
    return -(S.R-S.C);
  end Inf;

  function Sup(S: Interval) return Rep is
  begin
    return S.C+S.R;
  end Sup;

  pragma Inline (Inf,Sup);

  function SupAbs(S: Interval) return Rep is
  begin
    return Abs(S.C)+S.R;
  end SupAbs;

  function "<"(S1,S2: Interval) return Boolean is
  begin
    if Sup(S1)<Inf(S2) then return True; end if;
    if Sup(S2)<=Inf(S1) then return False; end if;
    raise Undefined;
    return False;
  end "<";

  function IntFloor(S: Interval) return Integer is
    R: constant Rep := Rep'Floor(Inf(S));
  begin
    if Rep'Floor(Sup(S)) /= R then
      raise Undefined;
    end if;
    return Integer(R);
  end IntFloor;

  function IntCeiling(S: Interval) return Integer is
    R: constant Rep := Rep'Ceiling(Sup(S));
  begin
    if Rep'Ceiling(Inf(S)) /= R then
      raise Undefined;
    end if;
    return Integer(R);
  end IntCeiling;

  function "abs"(S: Interval) return Interval is
    R: Rep;
  begin
    if S.C >= S.R then
      return S;
    elsif -S.C >= S.R then
      return (-S.C,S.R);
    end if;
    R := Half*(Abs(S.C)+S.R);
    return (R,R);
  end "abs";

  function Min(S1,S2: Interval) return Interval is
    U1: constant Rep := Sup(S1);
    L2: constant Rep := Inf(S2);
    L1,U2: Rep;
  begin
    if U1 <= L2 then
      return S1;
    end if;
    L1 := Inf(S1);
    U2 := Sup(S2);
    if U2 <= L1 then
      return S2;
    end if;
    return Scal(Rep'Min(L1,L2),Rep'Min(U1,U2));
  end Min;

  function Max(S1,S2: Interval) return Interval is
    U1: constant Rep := Sup(S1);
    L2: constant Rep := Inf(S2);
    L1,U2: Rep;
  begin
    if U1 <= L2 then
      return S2;
    end if;
    L1 := Inf(S1);
    U2 := Sup(S2);
    if U2 <= L1 then
      return S1;
    end if;
    return Scal(Rep'Max(L1,L2),Rep'Max(U1,U2));
   end Max;

  function Cap(R: Radius; S: Interval) return Interval is
-- bound on characteristic_function for closed ball of radius R
  begin
    if Abs(S.C)+S.R <= R then
      return S;
    end if;
    return Scal(Rep'Max(-R,Rep'Min(Zero,Inf(S))),Rep'Min(R,Rep'Max(Zero,Sup(S))));
  end Cap;

  function Up(R: Rep; Dummy: Interval) return Rep is
-- only used in high level procedures to increase domains
  begin
    return R+Epsilon;         -- avoid Gnat's Succ in Round_Up mode
  end Up;

  function Widen(S: Interval) return Interval is
-- only used in high level procedures to increase norms
  begin
    return (S.C,S.R+Epsilon); -- avoid Gnat's Succ in Round_Up mode
  end Widen;

  procedure ErrMult(R: in Rep; S: in out Interval) is
  begin
    S.R := Abs(R)*S.R;
  end ErrMult;

  function "+"(S,T: Interval) return Interval is
    Sum: constant Rep := S.C+T.C;
  begin
    return (Sum,S.R+T.R+(Sum+(-S.C-T.C)));
  end "+";

  function "-"(S,T: Interval) return Interval is
    Diff: constant Rep := S.C-T.C;
  begin
    return (Diff,S.R+T.R+(Diff+(T.C-S.C)));
  end "-";

  function "*"(R: Rep; S: Interval) return Interval is
    Prod: constant Rep := R*S.C;
  begin
    return (Prod,Abs(R)*S.R+(Prod+(-R)*S.C));
  end "*";

  function "/"(S: Interval; R: Rep) return Interval is
    Quot: constant Rep := S.C/R;
  begin
    return (Quot,S.R/Abs(R)+(Quot+S.C/(-R)));
  end "/";

  function Sqr(S: Interval) return Interval is
    Prod: constant Rep := S.C*S.C;
  begin
    return (Prod,(Two*Abs(S.C)+S.R)*S.R+(Prod+S.C*(-S.C)));
  end Sqr;

  function "*"(S,T: Interval) return Interval is
    Prod: constant Rep := S.C*T.C;
  begin
    return (Prod,Abs(S.C)*T.R+S.R*(Abs(T.C)+T.R)+(Prod+S.C*(-T.C)));
  end "*";

  function Inv(S: Interval) return Interval is
    Quot: constant Rep := One/S.C;
  begin
    if Abs(S.C) <= S.R then
      raise Constraint_Error;               -- division by zero
    end if;
    return (Quot,S.R/Abs(S.C)/(-(S.R-Abs(S.C)))+(Quot+NegOne/S.C));
  end Inv;

  function "/"(S,T: Interval) return Interval is
  begin
    return S*Inv(T);
  end "/";

  function "**"(S: Interval; I: Integer) return Interval is
  begin
    if I>0 then
      declare
        J: Integer := I;
        X,Y: Interval;
      begin
        X := S;
        if (J rem 2)=0 then
          Y := Interval_One;
        else
          Y := X;
        end if;
        J := J/2;
        while J>0 loop
          X := Sqr(X);
          if (J rem 2)>0 then
            Y := X*Y;
          end if;
          J := J/2;
        end loop;
        return Y;
      end;
    elsif I<0 then
      return Inv(S)**(-I);
    else
      return Interval_One;
    end if;
  end "**";

  procedure CoshSinh(S: in Interval; Sc,Ss: out Interval) is
    Iter: constant Integer := 6;
  begin
    if S = Interval_Zero then
      Sc := Interval_One;
      Ss := S;
    else
      declare
        K: Integer;
        R: Rep;
        X,T: Interval;
      begin
        K := 0;
        R := SupAbs(S);
        while R>Half loop
          K := K+1;
          R := R*Half;
        end loop;
        X := (Half**K)*S;
        T := Sqr(X);
        Sc := (One,Half);
        Ss := Sc;
        for I in reverse 1 .. Iter loop
          Sc := Interval_One+Sc*T/Rep((2*I)*(2*I-1));
          Ss := Interval_One+Ss*T/Rep((2*I)*(2*I+1));
        end loop;
        Ss := X*Ss;
        for I in 1 .. K loop
          T := Sc;
          Sc := Interval_One+Two*Sqr(Ss);
          Ss := Two*(T*Ss);
        end loop;
      end;
    end if;
  end CoshSinh;

  function Cosh(S: Interval) return Interval is
    Sc,Ss: Interval;
  begin
    CoshSinh(S,Sc,Ss);
    return Sc;
  end Cosh;

  function Sinh(S: Interval) return Interval is
    Sc,Ss: Interval;
  begin
    CoshSinh(S,Sc,Ss);
    return Ss;
  end Sinh;

  function Exp(S: Interval) return Interval is
  begin
    if S=Interval_Zero then return Scal(One); end if;   -- S=0
    if Sup(S)<Zero then return Inv(Exp(-S)); end if;    -- S<0
    declare
      X,Y: Interval;
    begin
      if Inf(S)<Zero then                               -- S contains 0
        CoshSinh(S,X,Y);
        return X+Y;
      end if;
      declare                                           -- S>0
        Iter: constant Integer := 12;
        K: Integer := 0;
        R: Rep;
      begin
        R := Sup(S);
        while R>Quarter loop
          K := K+1;
          R := R*Half;
        end loop;
        X := (Half**K)*S;
        Y := (Zero,Three/Two);
        for I in reverse 1 .. Iter loop
          Y := Interval_One+(X*Y)/Rep(I);
        end loop;
        for I in 1 .. K loop
          Y := Sqr(Y);
        end loop;
        return Y;
      end;
    end;
  end Exp;

  function Log(S: Interval) return Interval is
    use Reps.Ops;
  begin
    if S=Interval_One then return Interval_Zero; end if;
    declare
      Iter: constant Integer := 3;
      S0,LogS0,X,Y: Interval;
    begin
      LogS0 := Scal(Log(Approx(S)));   --- guess
      S0 := Exp(LogS0);
      X := (S0-S)/S0;                  --- Log(S) = Log(S0)+Log(1-X)
      Y := Interval_One-Scal(SupAbs(X));
      if Inf(Y) <= Zero then raise Constraint_Error; end if;
      Y := PlusMinus*Inv(Y)/Rep(Iter+1);
      for I in reverse 1 .. Iter loop
        Y := Inv(Scal(I))+X*Y;
      end loop;
      return LogS0-X*Y;
    end;
  end Log;

  function Lambda(M,N: Integer) return Interval is
    S: Interval;
  begin
    S := (Three/Two,Half);
    for I in 1 .. 100 loop
      S := Interval_One+Inv(S);
    end loop;
    return Scal(M)-Rep(N)*Inv(S);
  end Lambda;

begin

  Round_Up;

end Intervals.Ops;