File : balls.adb

with Reps.Ops, Reps.IO, Rounding;
use Reps.IO, Rounding;

pragma Optimize (Off);
pragma Elaborate_All (Rounding);

package body Balls is

--- basic

  function IsNumeric(S: Ball) return Boolean is
  begin
    return False;
  end IsNumeric;

  function IsSharp(S: Ball) return Boolean is
  begin
    return (S.R=Zero) and then (S.B=Zero);
  end IsSharp;

  function IsNumbers(S: Ball) return Boolean is
  begin
    return (S.B=Zero);
  end IsNumbers;

  procedure ModNumbers(S: in out Ball) is
  begin
    S.C := Zero;
    S.R := Zero;
  end ModNumbers;

  procedure ToBalls(S: in out Ball) is
  begin
    S.B := (abs(S.C)+S.R)+S.B;
    S.R := Zero;
    S.C := Zero;
  end ToBalls;

  function Approx(S: Ball) return Rep is
  begin
    return S.C;
  end Approx;

--- sets

  function Center0(S: Ball) return Boolean is
  begin
    return (S.C=Zero);
  end Center0;

  function Contains0(S: Ball) return Boolean is
    RB: constant Rep := S.R+S.B;
  begin
    if abs(S.C)>RB then return False; end if;
    if (abs(S.C)-RB)>Zero then raise Undefined; end if;
    return True;
  end Contains0;

  function Contains(S1,S2: Ball) return Boolean is
  begin
    if S2.B>S1.B then return False; end if;
    if S2.R>S1.R then return False; end if;
    if (((S1.B-S2.B)+(S1.R-S2.R))-abs(S2.C-S1.C))<Zero then return False; end if;
    if (((S2.B-S1.B)+(S2.R-S1.R))+abs(S2.C-S1.C))>Zero then raise Undefined; end if;
    return True;
  end Contains;

  function BallAt0(R: Rep) return Ball is
  begin
    return (Zero,Zero,abs(R));
  end BallAt0;

  function BallAt0(S: Ball) return Ball is
  begin
    return (Zero,Zero,(abs(S.C)+S.R)+S.B);
  end BallAt0;

  function DiskAt0(R: Rep) return Ball is
  begin
    return (Zero,abs(R),Zero);
  end DiskAt0;

  function DiskAt0(S: Ball) return Ball is
  begin
    return (Zero,(abs(S.C)+S.R)+S.B,Zero);
  end DiskAt0;

  function Disk(R1,R2: Rep) return Ball is
    Sum: constant Rep := R1+R2;
    Eps: constant Rep := Sum+(-R1-R2);
  begin
    return (Half*Sum,Half*(abs(R2-R1)+Eps),Zero);
  end Disk;

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

  function Center(S: Ball) return Ball is
  begin
    return (S.C,Zero,Zero);
  end Center;

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

  function ModCenter(S: Ball) return Ball is
  begin
    return (Zero,S.R,S.B);
  end ModCenter;

  function Union(S1,S2: Ball) return Ball is
    R1:  constant Rep := RMin(-(S1.R-S1.C),-(S2.R-S2.C));
    R2:  constant Rep := RMax(S1.C+S1.R,S2.C+S2.R);
    Sum: constant Rep := R1+R2;
    Eps: constant Rep := Sum+(-R1-R2);
  begin
    return (Half*Sum,Half*((R2-R1)+Eps),RMax(S1.B,S2.B));
  end Union;

  function Intersection(S1,S2: Ball) return Ball is
    R1:  constant Rep := RMax(-(S1.R-S1.C),-(S2.R-S2.C));
    R2:  constant Rep := RMin(S1.C+S1.R,S2.C+S2.R);
    Sum: constant Rep := R1+R2;
    Eps: constant Rep := Sum+(-R1-R2);
  begin
    if (S1.B /= Zero) or (S2.B /= Zero) then raise Undefined; end if;
    if R2<R1 then raise Undefined; end if;
    return (Half*Sum,Half*((R2-R1)+Eps),Zero);
  end Intersection;

--- order

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

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

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

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

  function Min(S1,S2: Ball) return Ball is
    U1: constant Rep := Sup(S1);
    L2: constant Rep := Inf(S2);
  begin
    if U1 <= L2 then return S1; end if;
    declare
      L1: constant Rep := Inf(S1);
      U2: constant Rep := Sup(S2);
    begin
      if U2 <= L1 then
        return S2;
      else
        return Disk(RMin(L1,L2),RMin(U1,U2));
      end if;
    end;
  end Min;

  function Max(S1,S2: Ball) return Ball is
    U1: constant Rep := Sup(S1);
    L2: constant Rep := Inf(S2);
  begin
    if U1 <= L2 then return S2; end if;
    declare
      L1: constant Rep := Inf(S1);
      U2: constant Rep := Sup(S2);
    begin
      if U2 <= L1 then
        return S1;
      else
        return Disk(RMax(L1,L2),RMax(U1,U2));
      end if;
    end;
  end Max;

  function Inf(S: Ball) return Rep is
  begin
    if S.B /= Zero then raise Undefined; end if;
    return -(S.R-S.C);
  end Inf;

  function Sup(S: Ball) return Rep is
  begin
    if S.B /= Zero then raise Undefined; end if;
    return S.C+S.R;
  end Sup;

--- addition and multiplication

  function IsZero(S: Ball) return Boolean is
  begin
    return (S=Ball_Zero);
  end IsZero;

  procedure SetZero(S: out Ball) is
  begin
    S := Ball_Zero;
  end SetZero;

  procedure Neg(S: in out Ball) is
  begin
    S.C := -S.C;
  end Neg;

  procedure Neg(S1: in Ball; S2: out Ball) is
  begin
    S2.C := -S1.C;
    S2.R :=  S1.R;
    S2.B :=  S1.B;
  end Neg;

  function "-"(S: Ball) return Ball is
  begin
    return (-S.C,S.R,S.B);
  end "-";

  procedure Add(S1: in Ball; S2: in out Ball) is
    Sum: constant Rep := S1.C+S2.C;
    Eps: constant Rep := Sum+(-S1.C-S2.C);
  begin
    S2.C := Sum;
    S2.R := (S1.R+S2.R)+Eps;
    S2.B := S1.B+S2.B;
  end Add;

  procedure Sum(S1,S2: in Ball; S3: out Ball) is
    Sum: constant Rep := S1.C+S2.C;
    Eps: constant Rep := Sum+(-S1.C-S2.C);
  begin
    S3.C := Sum;
    S3.R := (S1.R+S2.R)+Eps;
    S3.B := S1.B+S2.B;
  end Sum;

  function "+"(S1,S2: Ball) return Ball is
    Sum: constant Rep := S1.C+S2.C;
    Eps: constant Rep := Sum+(-S1.C-S2.C);
  begin
    return (Sum,(S1.R+S2.R)+Eps,S1.B+S2.B);
  end "+";

  procedure Sub(S1: in Ball; S2: in out Ball) is
    Dif: constant Rep := S2.C-S1.C;
    Eps: constant Rep := Dif+(S1.C-S2.C);
  begin
    S2.C := Dif;
    S2.R := (S1.R+S2.R)+Eps;
    S2.B := S1.B+S2.B;
  end Sub;

  procedure Diff(S1,S2: in Ball; S3: out Ball) is
    Dif: constant Rep := S1.C-S2.C;
    Eps: constant Rep := Dif+(S2.C-S1.C);
  begin
    S3.C := Dif;
    S3.R := (S1.R+S2.R)+Eps;
    S3.B := S1.B+S2.B;
  end Diff;

  function "-"(S1,S2: Ball) return Ball is
    Dif: constant Rep := S1.C-S2.C;
    Eps: constant Rep := Dif+(S2.C-S1.C);
  begin
    return (Dif,(S1.R+S2.R)+Eps,S1.B+S2.B);
  end "-";

  procedure Mult(I: in Integer; S: in out Ball) is
  begin
    Mult(Rep(I),S);
  end Mult;

  procedure Prod(I: in Integer; S1: in Ball; S2: out Ball) is
  begin
    Prod(Rep(I),S1,S2);
  end Prod;

  function "*"(I: Integer; S: Ball) return Ball is
    T: Ball;
  begin
    Prod(Rep(I),S,T);
    return T;
  end "*";

  procedure MultAdd(I: Integer; S1: in Ball; S2: in out Ball) is
    T: Ball;
  begin
    Prod(Rep(I),S1,T);
    Add(T,S2);
  end MultAdd;

  procedure Mult(R: in Rep; S: in out Ball) is
    Pro: constant Rep := R*S.C;
    Eps: constant Rep := Pro+R*(-S.C);
  begin
    S.C := Pro;
    S.R := abs(R)*S.R+Eps;
    S.B := abs(R)*S.B;
  end Mult;

  procedure Prod(R: in Rep; S1: in Ball; S2: out Ball) is
    Pro: constant Rep := R*S1.C;
    Eps: constant Rep := Pro+R*(-S1.C);
  begin
    S2.C := Pro;
    S2.R := abs(R)*S1.R+Eps;
    S2.B := abs(R)*S1.B;
  end Prod;

  function "*"(R: Rep; S: Ball) return Ball is
    T: Ball;
  begin
    Prod(R,S,T);
    return T;
  end "*";

  procedure MultAdd(R: in Rep; S1: in Ball; S2: in out Ball) is
    T: Ball;
  begin
    Prod(R,S1,T);
    Add(T,S2);
  end MultAdd;

  procedure Mult(S1: in Ball; S2: in out Ball) is
    T: Ball;
  begin
    Prod(S1,S2,T);
    S2 := T;
  end Mult;

  procedure Prod(S1,S2: in Ball; S3: out Ball) is
    R1:  constant Rep := abs(S1.C);
    R2:  constant Rep := abs(S2.C)+S2.R;
    Pro: constant Rep := S1.C*S2.C;
    Eps: constant Rep := Pro+S1.C*(-S2.C);
  begin
    S3.C := Pro;
    S3.R := (R1*S2.R+S1.R*R2)+Eps;
    S3.B := (R1+S1.R)*S2.B+S1.B*(R2+S2.B);
  end Prod;

  function "*"(S1,S2: Ball) return Ball is
    T: Ball;
  begin
    Prod(S1,S2,T);
    return T;
  end "*";

  procedure MultAdd(S1,S2: in Ball; S3: in out Ball) is
    T: Ball;
  begin
    Prod(S1,S2,T);
    Add(T,S3);
  end MultAdd;

  function "/"(S: Ball; R: Rep) return Ball is
    Quo: constant Rep := S.C/R;
    Eps: constant Rep := Quo+S.C/(-R);
  begin
    return (Quo,S.R/abs(R)+Eps,S.B/abs(R));
  end "/";

  procedure Inv(S1: in Ball; S2: out Ball) is
    R1:  constant Rep := S1.R-Abs(S1.C);
    R2:  constant Rep := S1.B+R1;
    Quo: constant Rep := One/S1.C;
    Eps: constant Rep := Quo+NegOne/S1.C;
  begin
    if R2 >= Zero then raise Constraint_Error; end if;  -- division by zero
    S2.C := Quo;
    S2.R := (S1.R/abs(S1.C))/(-R1)+Eps;
    S2.B := (S1.B/(-R1))/(-R2);
  end Inv;

  function Inv(S: Ball) return Ball is
    T: Ball;
  begin
    Inv(S,T);
    return T;
  end Inv;

  function "/"(S1,S2: Ball) return Ball is
    T1,T2: Ball;
  begin
    Inv(S2,T1);
    Prod(S1,T1,T2);
    return T2;
  end "/";

  function "**"(S: Ball; I: Integer) return Ball is
    J: Integer := I;
    X,Y: Ball;
  begin
    if I>0 then
      X := S;
      if (J rem 2)=0 then
        Y := Ball_One;
      else
        Y := X;
      end if;
      J := J/2;
      while J>0 loop
        Mult(X,X);           -- identical args ok with our Mult
        if (J rem 2)>0 then
          Mult(X,Y);
        end if;
        J := J/2;
      end loop;
      return Y;
    elsif I<0 then
      return Inv(S)**(-I);
    else
      return Ball_One;
    end if;
  end "**";

--- functions

  function Sqr(S: Ball) return Ball is
    T: Ball;
  begin
    Prod(S,S,T);
    return T;
  end Sqr;

  function Norm(S: Ball) return Ball is
    T: Ball;
  begin
    T.C := abs(S.C);
    T.R := S.R+S.B;
    T.B := Zero;
    if T.R>T.C then
      T.C := Half*(T.C+T.R);
      T.R := T.C;
    end if;
    return T;
  end Norm;

  function MaxNorm(S: Ball) return Radius is
  begin
    return (abs(S.C)+S.R)+S.B;
  end MaxNorm;

--- could trust/use hardware for these functions

  function ExpUp(R: Rep) return Rep is
    Iter: constant Integer := 13;
    K: Integer := 0;
    X: Rep := R;
    Y: Rep := Rep(33)/Rep(32);  -- assumes Iter>8
  begin
    if R<Zero then return One/ExpDown(-R); end if;
    if R=Zero then return One; end if;
    while X>Quarter loop
      K := K+1;
      X := Half*X;
    end loop;
    for I in reverse 1 .. Iter loop
      Y := One+(X*Y)/Rep(I);
    end loop;
    for I in 1 .. K loop
      Y := Y*Y;
    end loop;
    return Y;
  end ExpUp;

  function ExpDown(R: Rep) return Rep is
    Iter: constant Integer := 13;
    K: Integer := 0;
    X: Rep := R;
    Y: Rep := NegOne;
  begin
    if R<Zero then return -(NegOne/ExpUp(-R)); end if;
    if R=Zero then return One; end if;
    while X>Quarter loop
      K := K+1;
      X := Half*X;
    end loop;
    for I in reverse 1 .. Iter loop
      Y := NegOne+(X*Y)/Rep(I);
    end loop;
    for I in 1 .. K loop
      Y := (-Y)*Y;
    end loop;
    return -Y;
  end ExpDown;

  function Exp(S: Ball) return Ball is
  begin
    return Disk(ExpDown(Inf(S)),ExpUp(Sup(S)));
  end Exp;

  function Log(S: Ball) return Ball is
    use Reps.Ops;
    Iter: constant Integer := 3;
    S0,LogS0,X,Y: Ball;
  begin
    if S=Ball_One then return Ball_Zero; end if;
    if S.B /= Zero then raise Undefined; end if;
    LogS0 := Scal(Log(Approx(S)));                  -- guess
    S0 := Exp(LogS0);
    X := (S0-S)/S0;                                 -- Log(S) = Log(S0)+Log(1-X)
    Y := Ball_One-Norm(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 Log;

  function Sqrt(S: Ball) return Ball is
    use Reps.Ops;
    S0,SqrtS0,X,Y: Ball;
  begin
    if S=Ball_Zero then return Ball_Zero; end if;
    if Inf(S)<Zero then raise Undefined; end if;    -- also excludes nonzero S.B
    SqrtS0 := Scal(Sqrt(Approx(S)));                -- guess
    S0 := Sqr(SqrtS0);
    X := (S-S0)/S0;                                 -- Sqrt(S) = Sqrt(S0)*Sqrt(1+X)
    Y := ((PlusMinus*X-Scal(2))*X+Scal(8))*X/Rep(16);
    return SqrtS0*(Ball_One+Y);
  end Sqrt;

--- conversion and i/o

  function Scal(K: Integer) return Ball is
  begin
    return (Rep(K),Zero,Zero);
  end Scal;

  function Scal(R: Rep) return Ball is
  begin
    return (R,Zero,Zero);
  end Scal;

  procedure ShowType(S: in Ball) is
  begin
    Show0("Scalar type is Ball");
  end ShowType;

  procedure Show1(N: in String; S: in Ball) is
  begin
    Show3(N,S.C,S.R,S.B);
  end Show1;

  procedure Show2(N: in String; S1,S2: in Ball) is
  begin
    Show3(N,S1.C,S1.R,S1.B,False);
    Show3("   ",S2.C,S2.R,S2.B,True);
  end Show2;

  procedure Get(F: in File_Type; S: out Ball; Decimal: in Boolean := False) is
  begin
    Get(F,S.C,S.R,S.B,Decimal);
  end Get;

  procedure Put(F: in File_Type; S: in Ball; Decimal: in Boolean := False) is
  begin
    Put(F,S.C,S.R,S.B,Decimal);
  end Put;

  procedure UpMultAdd(R1,R2: in Rep; R3: in out Rep) is
  begin
    R3 := R3+R1*R2;
  end UpMultAdd;

begin

  Restore_Checks;

end Balls;