with Rounding.FPU, Globals, Strings, IPowers, QPowers, Roots, Flts.Std.EFuns;
use Globals, Strings, Flts.Std.EFuns;

pragma Elaborate_All (Rounding.FPU,Globals,Strings,IPowers,QPowers,Roots,Flts.Std.EFuns);
pragma Optimize (Off);

package body Flts.Std.Balls is

  package Txt_IO renames Ada.Text_IO;

  --- auxiliary

  function ReturnTrue(Dummy: Ball) return Boolean is
    pragma Unreferenced(Dummy);
  begin
    return True;
  end ReturnTrue;

  procedure NOP(Dummy: in out Ball) is
    pragma Unreferenced(Dummy);
  begin
    null;
  end NOP;

  --- basic

  function Info(Dummy: Ball) return Scalar_Info is
    pragma Unreferenced(Dummy);
  begin
    return BInfo;
  end Info;

  function IsSharp(S: Ball) return Boolean is
  begin
    return S.R=RZero;
  end IsSharp;

  function IsZero(S: Ball) return Boolean is
  begin
    return S.C=RZero and then S.R=RZero;
  end IsZero;

  function "="(S1,S2: Ball) return Boolean is
  begin
    return (S1.C=S2.C) and then (S1.R=S2.R);
  end "=";

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

  function ReturnZero(Dummy: in Ball) return Ball is
    pragma Unreferenced(Dummy);
  begin
    return BZero;
  end  ReturnZero;

  function Copy(S: Ball) return Ball is
  begin
    return S;
  end Copy;

  procedure Copy(S1: in Ball; S2: in out Ball) is
  begin
    S2 := S1;
  end Copy;

  procedure Swap(S1,S2: in out Ball) is
    S: constant Ball := S1;
  begin
    S1 := S2;
    S2 := S;
  end Swap;

  --- sets

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

  function Contains0(S: Ball) return Logical is
  begin
    if S.R<abs(S.C) then
      return False;
    else
      return True;
    end if;
  end Contains0;

  function Contains(S1,S2: Ball) return Logical is
    Y: Rep;
  begin
    if S1.R<S2.R then return False; end if;
    if S2.C=S1.C then return True; end if;
    Y := abs(S2.C-S1.C);
    if ((S1.R-S2.R)-Y)<RZero then return False; end if;
    if ((S2.R-S1.R)+Y)>RZero then
      return Uncertain;
    end if;
    return True;
  end Contains;

  procedure BallAt0(F: in Flt; S: in out Ball) is
  begin
    S.R := abs(Rad(F));
    S.C := RZero;
  end BallAt0;

  function BallAt0(F: Flt) return Ball is
  begin
    return (RZero,abs(Rad(F)));
  end BallAt0;

  procedure ToErr(S: in out Ball) is
  begin
    S.R := abs(S.C)+S.R;
    S.C := RZero;
  end ToErr;

  procedure ToErr(S1: in Ball; S2: in out Ball) is
  begin
    S2.R := abs(S1.C)+S1.R;
    S2.C := RZero;
  end ToErr;

  function ToErr(S: Ball) return Ball is
  begin
    return (RZero,abs(S.C)+S.R);
  end ToErr;

  procedure Center(S: in out Ball) is
  begin
    S.R := RZero;
  end Center;

  procedure Center(S1: in Ball; S2: in out Ball) is
  begin
    S2.C := S1.C;
    S2.R := RZero;
  end Center;

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

  procedure ModCenter(S: in out Ball) is
  begin
    S.C := RZero;
  end ModCenter;

  procedure ModCenter(S1: in Ball; S2: in out Ball) is
  begin
    S2.C := RZero;
    S2.R := S1.R;
  end ModCenter;

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

  procedure ErrMult(B: in Radius; S: in out Ball) is
  begin
    S.R := abs(Rad(B)*S.R);
  end ErrMult;

  procedure Union(S1: in Ball; S2: in out Ball) is
  begin
    Enclose(RMin(Inf(S1),Inf(S2)),RMax(Sup(S1),Sup(S2)),S2);
  end Union;

  function Union(S1,S2: Ball) return Ball is
    S3: Ball;
  begin
    Enclose(RMin(Inf(S1),Inf(S2)),RMax(Sup(S1),Sup(S2)),S3);
    return S3;
  end Union;

  procedure Intersection(S1: in Ball; S2: in out Ball; Empty: out Logical) is
    R1: Rad renames S1.R;
    R2: Rad renames S2.R;
    C1: Rep renames S1.C;
    C2: Rep renames S2.C;
    YL,YU: Rep;
  begin
    Empty := Uncertain;
    if C1=C2 then
      Empty := False;
      if R1<R2 then R2 := R1; end if;
    else
      YL := Max(C1-R1,C2-R2);
      YU := -Max(-C1-R1,-C2-R2);
      if YL <= YU then
        Empty := False;
      end if;
      YL := Max(-(R1-C1),-(R2-C2));
      YU := Min(C1+R1,C2+R2);
      if YU < YL then
        Empty := True;
      end if;
      C2 := RHalf*(YL+YU);
      R2 := (YU-YL)+(C2+RHalf*(-YL-YU));
    end if;
  end Intersection;

  --- order

  function Sign(S: Ball) return Integer is
  begin
    if Inf(S)>Zero then return 1; end if;
    if Sup(S)<Zero then return -1; end if;
    if not IsZero(S) then
      raise Not_Certain;
    end if;
    return 0;
  end Sign;

  function Compare(S1,S2: Ball) return Integer is
  begin
    if Inf(S1)>Sup(S2) then return 1; end if;
    if Sup(S1)<Inf(S2) then return -1; end if;
    if S1 /= S2 then
      raise Not_Certain;
    end if;
    return 0;
  end Compare;

  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 Not_Certain;
    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 Not_Certain;
    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 Not_Certain;
    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 Not_Certain;
    end if;
    return False;
  end ">=";

  procedure Min(S1: in Ball; S2: in out Ball) is
    L1: constant Flt := Inf(S1);
    U2: constant Flt := Sup(S2);
  begin
    if U2>L1 then
      declare
        U1: constant Flt := Sup(S1);
        L2: constant Flt := Inf(S2);
      begin
        if U1>L2 then
          Enclose(RMin(L1,L2),RMin(U1,U2),S2);
        else
          S2 := S1;
        end if;
      end;
    end if;
  end Min;

  procedure Min(S1,S2: in Ball; S3: in out Ball) is
    L1: constant Flt := Inf(S1);
    U2: constant Flt := Sup(S2);
  begin
    if U2>L1 then
      declare
        U1: constant Flt := Sup(S1);
        L2: constant Flt := Inf(S2);
      begin
        if U1>L2 then
          Enclose(RMin(L1,L2),RMin(U1,U2),S3);
        else
          S3 := S1;
        end if;
      end;
    else
      S3 := S2;
    end if;
  end Min;

  function Min(S1,S2: Ball) return Ball is
    S3: Ball;
  begin
    Min(S1,S2,S3);
    return S3;
  end Min;

  procedure Max(S1: in Ball; S2: in out Ball) is
    U1: constant Flt := Sup(S1);
    L2: constant Flt := Inf(S2);
  begin
    if U1>L2 then
      declare
        L1: constant Flt := Inf(S1);
        U2: constant Flt := Sup(S2);
      begin
        if U2>L1 then
          Enclose(RMax(L1,L2),RMax(U1,U2),S2);
        else
          S2 := S1;
        end if;
      end;
    end if;
  end Max;

  procedure Max(S1,S2: in Ball; S3: in out Ball) is
    U1: constant Flt := Sup(S1);
    L2: constant Flt := Inf(S2);
  begin
    if U1>L2 then
      declare
        L1: constant Flt := Inf(S1);
        U2: constant Flt := Sup(S2);
      begin
        if U2>L1 then
          Enclose(RMax(L1,L2),RMax(U1,U2),S3);
        else
          S3 := S1;
        end if;
      end;
    else
      S3 := S2;
    end if;
  end Max;

  function Max(S1,S2: Ball) return Ball is
    S3: Ball;
  begin
    Max(S1,S2,S3);
    return S3;
  end Max;

  function Inf(S: Ball) return Flt is
  begin
    return -Flt(S.R-S.C);
  end Inf;

  function Sup(S: Ball) return Flt is
  begin
    return Flt(S.C+S.R);
  end Sup;

  --- addition and multiplication

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

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

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

  procedure Add(I: in Integer; S: in out Ball) is
    C: Rep renames S.C;
    NegC: constant Rep := Rep(-I)-C;
  begin
    C := Rep(I)+C;
    S.R := S.R+(C+NegC);
  end Add;

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

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

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

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

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

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

  procedure Mult(R: in Flt; S,Dummy: in out Ball) is
    pragma Unreferenced(Dummy);
    C: Rep renames S.C;
    NegC: constant Rep := Rep(-R)*C;
  begin
    C := Rep(R)*C;
    S.R := abs(Rep(R))*S.R+(C+NegC);
  end Mult;

  procedure Mult(R: in Flt; S: in out Ball) is
    C: Rep renames S.C;
    NegC: constant Rep := Rep(-R)*C;
  begin
    C := Rep(R)*C;
    S.R := abs(Rep(R))*S.R+(C+NegC);
  end Mult;

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

  function "*"(R: Flt; S: Ball) return Ball is
    NegC: constant Rep := Rep(-R)*S.C;
    T: Ball;
  begin
    T.C := Rep(R)*S.C;
    T.R := abs(Rep(R))*S.R+(T.C+NegC);
    return T;
  end "*";

  procedure AddProd(R: in Flt; S1: in Ball; S2,Dummy: in out Ball) is
    pragma Unreferenced(Dummy);
    C: Rep renames S2.C;
    NegC: constant Rep := Rep(-R)*S1.C-C;
  begin
    C := Rep(R)*S1.C+C;
    S2.R := abs(Rep(R))*S1.R+S2.R+(C+NegC);
  end AddProd;

  procedure AddProd(R: in Flt; S1: in Ball; S2: in out Ball) is
    C: Rep renames S2.C;
    NegC: constant Rep := Rep(-R)*S1.C-C;
  begin
    C := Rep(R)*S1.C+C;
    S2.R := abs(Rep(R))*S1.R+S2.R+(C+NegC);
  end AddProd;

  procedure Mult(Q: in Rational; S: in out Ball) is
    N: constant LInt := LNum(Q);
    D: constant LInt := LDen(Q);
  begin
    if N /= 1 then Mult(Flt(N),S); end if;
    if D /= 1 then Div(Flt(D),S); end if;
  end Mult;

  procedure Mult(S1: in Ball; S2: in out Ball) is
    C2: Rep renames S2.C;
    R2: Rad renames S2.R;
    NegC: constant Rep := (-S1.C)*C2;
  begin
    R2 := abs(S1.C)*R2+S1.R*(abs(C2)+R2);
    C2 := S1.C*C2;
    R2 := R2+(C2+NegC);
  end Mult;

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

  function "*"(S1,S2: Ball) return Ball is
    NegC: constant Rep := (-S1.C)*S2.C;
    S3: Ball;
  begin
    S3.C := S1.C*S2.C;
    S3.R := abs(S1.C)*S2.R+S1.R*(abs(S2.C)+S2.R)+(S3.C+NegC);
    return S3;
  end "*";

  procedure AddProd(S1,S2: in Ball; S3,Dummy: in out Ball) is
    pragma Unreferenced(Dummy);
    C: Rep renames S3.C;
    NegC: constant Rep := (-S1.C)*S2.C-C;
  begin
    C := S1.C*S2.C+C;
    S3.R := abs(S1.C)*S2.R+S1.R*(abs(S2.C)+S2.R)+S3.R+(C+NegC);
  end AddProd;

  procedure AddProd(S1,S2: in Ball; S3: in out Ball) is
    C: Rep renames S3.C;
    NegC: constant Rep := (-S1.C)*S2.C-C;
  begin
    C := S1.C*S2.C+C;
    S3.R := abs(S1.C)*S2.R+S1.R*(abs(S2.C)+S2.R)+S3.R+(C+NegC);
  end AddProd;

  procedure SumProd(S1,S2,S3: in Ball; S4: in out Ball) is
    NegC: constant Rep := (-S1.C)*S2.C-S3.C;
  begin
    S4.C := S1.C*S2.C+S3.C;
    S4.R := abs(S1.C)*S2.R+S1.R*(abs(S2.C)+S2.R)+S3.R+(S4.C+NegC);
  end SumProd;

  procedure SubProd(S1,S2: in Ball; S3,Dummy: in out Ball) is
    pragma Unreferenced(Dummy);
    C: Rep renames S3.C;
    NegC: constant Rep := S1.C*S2.C-C;
  begin
    C := (-S1.C)*S2.C+C;
    S3.R := abs(S1.C)*S2.R+S1.R*(abs(S2.C)+S2.R)+S3.R+(C+NegC);
  end SubProd;

  procedure Div(R: in Flt; S,Dummy: in out Ball) is
    pragma Unreferenced(Dummy);
    C: Rep renames S.C;
    NegC: constant Rep := C/Rep(-R);
  begin
    C := C/Rep(R);
    S.R := S.R/abs(Rep(R))+(C+NegC);
  end Div;

  procedure Div(R: in Flt; S: in out Ball) is
    C: Rep renames S.C;
    NegC: constant Rep := C/Rep(-R);
  begin
    C := C/Rep(R);
    S.R := S.R/abs(Rep(R))+(C+NegC);
  end Div;

  procedure Quot(S1: in Ball; R: in Flt; S2: in out Ball) is
    NegC: constant Rep := S1.C/Rep(-R);
  begin
    S2.C := S1.C/Rep(R);
    S2.R := S1.R/abs(Rep(R))+(S2.C+NegC);
  end Quot;

  function "/"(S: Ball; R: Flt) return Ball is
    NegC: constant Rep := S.C/Rep(-R);
    T: Ball;
  begin
    T.C := S.C/Rep(R);
    T.R := S.R/abs(Rep(R))+(T.C+NegC);
    return T;
  end "/";

  procedure Div(S1: in Ball; S2: in out Ball) is
    S: Ball;
  begin
    Inv(S1,S);
    Mult(S,S2);
  end Div;

  procedure Quot(S1,S2: in Ball; S3: in out Ball) is
  begin
    Inv(S2,S3);
    Mult(S1,S3);
  end Quot;

  function "/"(S1,S2: Ball) return Ball is
    S3: Ball;
  begin
    Inv(S2,S3);
    Mult(S1,S3);
    return S3;
  end "/";

  procedure Inv(S1: in Ball; S2: in out Ball) is
    R:     constant Rep := S1.R-abs(S1.C);
    NegC:  constant Rep := Rep(-1)/S1.C;
  begin
    if R >= RZero then
      raise Undefined with "Flts.Std.Balls.Inv error: possible division by zero";
    end if;
    S2.C := Rep(1)/S1.C;
    S2.R := (S1.R/abs(S1.C))/(-R)+(S2.C+NegC);
  end Inv;

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

  --- functions

  procedure Norm(S1: in Ball; S2: in out Ball) is
  begin
    S2.C := abs(S1.C);
    if S2.C >= S1.R then
      S2.R := S1.R;
    else
      S2.C := RHalf*(S2.C+S1.R);
      S2.R := S2.C;
    end if;
  end Norm;

  function Norm(S: Ball) return Ball is
    R: Rep := abs(S.C);
  begin
    if R >= S.R then
      return (R,S.R);
    else
      R := RHalf*(R+S.R);
      return (R,R);
    end if;
  end Norm;

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

  procedure Sqr(S: in out Ball) is
    C: Rep renames S.C;
    R: Rad renames S.R;
    NegC: constant Rep := (-C)*C;
  begin
    R := (abs(C+C)+R)*R;
    C := C*C;
    R := R+(C+NegC);
  end Sqr;

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

  procedure Sqr(S1: in Ball; S2: in out Ball) is
    C1: Rep renames S1.C;
    NegC: constant Rep := (-C1)*C1;
  begin
    S2.C := C1*C1;
    S2.R := (abs(C1+C1)+S1.R)*S1.R+(S2.C+NegC);
  end Sqr;

  procedure Sqrt(S1: in Ball; S2: in out Ball) is
    --- IEEE 754 Sqrt has error <= ULP/2
    C1: Rep renames S1.C;
    R1: Rad renames S1.R;
  begin
    if R1 /= C1 then
      Enclose(Rep'Pred(Sqrt(-(R1-C1))),Sqrt(C1+R1),S2); -- raises exception if C1<R1
    elsif C1=RZero then
      S2 := BZero;
    else
      S2.C := Sqrt(RHalf*C1);
      S2.R := S2.C;
    end if;
  end Sqrt;

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

  procedure Root(K: Positive; S1: in Ball; S2: in out Ball) is
    YL,YU: Rep;
  begin
    Root(K,S1.C,S1.R,YL,YU);
    S2.C := RHalf*(YL+YU);
    S2.R := (YU-YL)+(S2.C+RHalf*(-YL-YU));
  end Root;

  function Root(K: Positive; S: Ball) return Ball is
    T: Ball;
  begin
    Root(K,S,T);
    return T;
  end Root;

  procedure Exp(S1: in Ball; S2: in out Ball) is
    YL,YU: Rep;
  begin
    Exp(S1.C,S1.R,YL,YU);
    S2.C := RHalf*(YL+YU);
    S2.R := (YU-YL)+(S2.C+RHalf*(-YL-YU));
  end Exp;

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

  procedure Log(S1: in Ball; S2: in out Ball) is
    YL,YU: Rep;
  begin
    Log(S1.C,S1.R,YL,YU);
    S2.C := RHalf*(YL+YU);
    S2.R := (YU-YL)+(S2.C+RHalf*(-YL-YU));
  end Log;

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

  procedure ArcCos(S1: in Ball; S2: in out Ball) is
    YL,YU: Rep;
  begin
    ArcCos(S1.C,S1.R,YL,YU);
    S2.C := RHalf*(YL+YU);
    S2.R := (YU-YL)+(S2.C+RHalf*(-YL-YU));
  end ArcCos;

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

  procedure ArcSin(S1: in Ball; S2: in out Ball) is
    YL,YU: Rep;
  begin
    ArcSin(S1.C,S1.R,YL,YU);
    S2.C := RHalf*(YL+YU);
    S2.R := (YU-YL)+(S2.C+RHalf*(-YL-YU));
  end ArcSin;

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

  procedure Cos(S1: in Ball; S2: in out Ball) is
    YL,YU: Rep;
  begin
    Cos(S1.C,YL,YU);
    S2.C := RHalf*(YL+YU);
    S2.R := S1.R+(YU-YL)+(S2.C+RHalf*(-YL-YU));
  end Cos;

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

  procedure Sin(S1: in Ball; S2: in out Ball) is
    YL,YU: Rep;
  begin
    Sin(S1.C,YL,YU);
    S2.C := RHalf*(YL+YU);
    S2.R := S1.R+(YU-YL)+(S2.C+RHalf*(-YL-YU));
  end Sin;

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

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

  function Pi return Ball is
    --- 2*ArcCos(0)
    YL,YU: Rep;
    S: Ball;
  begin
    ArcCos(RZero,RZero,YL,YU);
    S.R := -YU-YL;
    S.C := YU+YL;
    S.R := (YU-YL)+(S.C+S.R);
    return S;
  end Pi;

  --- conversion and i/o

  procedure Assign(I: in Integer; S: in out Ball) is
  begin
    S.C := Rep(I);
    S.R := RZero;
  end Assign;

  procedure Assign(Q: in Rational; S: in out Ball) is
  begin
    S.C := Rep(LNum(Q));
    S.R := RZero;
    Div(Flt(LDen(Q)),S);
  end Assign;

  procedure Assign(R: in Flt; S: in out Ball) is
  begin
    S.C := Rep(R);
    S.R := RZero;
  end Assign;

  function Scal(I: Integer) return Ball is
  begin
    return (Rep(I),RZero);
  end Scal;

  function Scal(Q: Rational) return Ball is
    S: Ball;
  begin
    S.C := Rep(LNum(Q));
    S.R := RZero;
    Div(Flt(LDen(Q)),S);
    return S;
  end Scal;

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

  procedure Enclose(R1,R2: in Rep; S: in out Ball) is
    H1:   constant Rep := RepHalf*R1;
    H2:   constant Rep := RepHalf*R2;
    NegC: constant Rep := -H1-H2;
  begin
    S.C := H1+H2;
    if H1 <= H2 then
      S.R := (H2-H1)+(S.C+NegC);
    else
      S.R := (H1-H2)+(S.C+NegC);
    end if;
  end Enclose;

  procedure Enclose(R1,R2: in Flt; S: in out Ball) is
    H1:   constant Rep := Rep(Half*R1);
    H2:   constant Rep := Rep(Half*R2);
    NegC: constant Rep := -H1-H2;
  begin
    S.C := H1+H2;
    if H1 <= H2 then
      S.R := (H2-H1)+(S.C+NegC);
    else
      S.R := (H1-H2)+(S.C+NegC);
    end if;
  end Enclose;

  function Enclose(R1,R2: Flt) return Ball is
    S: Ball;
  begin
    Enclose(R1,R2,S);
    return S;
  end Enclose;

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

  procedure Show1(N: in String; S: in Ball; NewLine: in Boolean := True) is
    Digs: constant Natural := 0; -- use default
  begin
    Show0(N & DecStr(S.C,Digs),False);
    if S.R /= RZero then
      Show0(" [" & DecStr(S.R,Digs) & "]",False);
    end if;
    Show0(NewLine);
  end Show1;

  procedure Show2(N: in String; S1,S2: in Ball; NewLine: in Boolean := True) is
  begin
    Show1(N,S1,False);
    Show1("   ",S2,NewLine);
  end Show2;

  procedure Put(F: in File_Type; S: in Ball; Decimal: in Boolean := False) is
  begin
    if Decimal then
      Txt_IO.Put_Line(F,DecStr(S.C));
    elsif S.R=RZero then
      Txt_IO.Put_Line(F,HexStr(S.C));
    else
      Txt_IO.Put_Line(F,HexStr(S.C) & "," & HexStr(S.R));
    end if;
  end Put;

  procedure Get(F: in File_Type; S: in out Ball; Decimal: in Boolean := False) is
    V: constant String := Strings.Get_Next_Line(F);
    Rounded: Boolean;
    N: Integer := V'First;
    E: Rad := RZero;
  begin
    SetZero(S);
    while N <= V'Last loop
      exit when V(N)=',';
      N := N+1;
    end loop;
    if Decimal then
      ValDec(V(V'First .. N-1),S.C);
    else
      ValHex(V(V'First .. N-1),S.C,Rounded);
      if Rounded then
        E := S.C+(Tiny-S.C); --- cheap
      end if;
      if N<V'Last then
        ValHex(V(N+1 .. V'Last),S.R,Rounded);
        if Rounded then
          S.R := S.R+Tiny;
        end if;
      end if;
      S.R := S.R+E;
    end if;
  end Get;

  procedure Write(FileName: in String; S: in Ball; Decimal: in Boolean := False) is
    F: File_Type;
  begin
    if Verbosity>0 then Show0("Writing " & FileName); end if;
    Create(F,Out_File,FileName);
    Put(F,S,Decimal);
    Close(F);
  end Write;

  procedure Read(FileName: in String; S: in out Ball; Decimal: in Boolean := False) is
    F: File_Type;
  begin
    if Verbosity>0 then Show0("Reading " & FileName); end if;
    Open(F,In_File,FileName);
    Get(F,S,Decimal);
    Close(F);
  end Read;

  procedure PutStd(S: in Ball; C,R: out LLFloat) is
  begin
    C := S.C;
    R := S.R;
  end PutStd;

  procedure PutStd(S: in Ball; C,R,B: out LLFloat) is
  begin
    C := S.C;
    R := S.R;
    B := RZero;
  end PutStd;

  procedure GetStd(C,R: in LLFloat; S: in out Ball) is
  begin
    S.C := C;
    S.R := R;
  end GetStd;

  procedure GetStd(C,R,B: in LLFloat; S: in out Ball) is
  begin
    S.C := C;
    S.R := R+B;
  end GetStd;

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

  procedure GetNum(R: in Rep; S: in out Ball) is
  begin
    S.C := R;
    S.R := RZero;
  end GetNum;

  --- misc

  function Get_Precision(S: Ball) return Positive is
  begin
    return Get_Precision(S.C);
  end Get_Precision;

  procedure Proper_Rounding(Dummy: in Ball) is
    pragma Unreferenced(Dummy);
  begin
    Rounding.FPU.Set_Rounding_Mode((Up,Strict));
  end Proper_Rounding;

  procedure Proper_Rounding is
  begin
     Rounding.FPU.Set_Rounding_Mode((Up,Strict));
  end Proper_Rounding;

  procedure Free_Cache(Dummy: in Ball) is
    pragma Unreferenced(Dummy);
  begin
    null;
  end Free_Cache;

  procedure Put_Numeric(N: in Positive; S: in Ball) is
  begin
    Put_Numeric(N,S.C);
  end Put_Numeric;

  procedure Get_Numeric(N: in Positive; S: in out Ball) is
  begin
    Get_Numeric(N,S.C);
    S.R := RZero;
  end Get_Numeric;

  --------------------------------------------------------------------
  --------------------------------------------------------------------
  package BIP is new IPowers (Scalar => Ball);
  package BQP is new QPowers (Scalar => Ball, IPower => BIP.IPower);
  --------------------------------------------------------------------
  --------------------------------------------------------------------

  procedure IPower(I: in Integer; S1: in Ball; S2,Tmp: in out Ball) renames BIP.Ipower;

  function "**"(S: Ball; I: Integer) return Ball is
    T: Ball;
  begin
    BIP.IPower(I,S,T);
    return T;
  end "**";

  procedure QPower(Q: in Rational; S1: in Ball; S2: in out Ball) renames BQP.QPower;

  -------------------------------------------
  -------------------------------------------
  package SR is new Roots (Scalar => Ball);
  -------------------------------------------
  -------------------------------------------

  procedure Roots2(B,C: in Ball; U1,U2,V: in out Ball) is
  begin
    SR.NewtonRoots2(B,C,U1,U2,V,64);
  end Roots2;

  procedure Roots3(B,C,D: in Ball; U0,U1,U2,V: in out Ball) is
  begin
    SR.NewtonRoots3(B,C,D,U0,U1,U2,V,64);
  end Roots3;

  procedure Rounding_Check(P: in Positive) is
    --- assumes Rounding_Mode=(Up,Strict) and P=3
    U: constant Rep := RepOne/Rep(P);
    V: constant Rep := -(RepOne/Rep(-P));
  begin
    if U=V then
      Show0("Flts.Std.Balls.Rounding_Check error:");
      Show0("Your compiler interferes with rounding by re-arranging expressions.");
      Show0("With gcc use the flag -frounding-math for IEEE 754 compliance.");
      raise Sorry;
    end if;
  end Rounding_Check;

begin

  Flts.A_Numeric_Mode := Rounding.FPU.Default_Rounding_Mode'Access;
  Rounding.FPU.Set_Rounding_Mode((Up,Strict));
  Rounding_Check(3);

end Flts.Std.Balls;
