with Rounding.FPU, Strings, IPowers, QPowers, Roots, MPFR.Floats.EFuns;
use Strings, MPFR.Floats.EFuns;

pragma Elaborate_All (Rounding.FPU,Strings,IPowers,QPowers,Roots,MPFR.Floats.EFuns);
pragma Optimize (Off);

package body MPFR.Floats.Balls is

  package Txt_IO renames Ada.Text_IO;

  --- basic rounding

  procedure Re_Format(P: in MP_Prec_T; V: in out MPFR_T) is
    --- increase precision if necessary (erases value)
  begin
    if V(1).MP_Prec < P then
      MPFR_Set_Prec(V,P);
    end if;
  end Re_Format;

  procedure Re_Format(V: in out MPFR_T) is
    --- reformat with precision MP_Precision if necessary (erases value)
  begin
    if V(1).MP_Prec /= MP_Precision then
      MPFR_Set_Prec(V,MP_Precision);
    end if;
  end Re_Format;

  procedure Fix_Precision(V: in out MPFR_T; Rounded: out Int) is
  begin
    if V(1).MP_Prec=MP_Precision then
      Rounded := 0;
    else
      declare
        W: MPFR_T := V;
      begin
        MPFR_Init2(V,MP_Precision);
        MPFR_Set(Rounded,V,W,GMP_RNDN);
        MPFR_Clear(W);
      end;
    end if;
  end Fix_Precision;

  function ULP_Weight(V: MPFR_T) return LLFloat is
    --- weight of the ulp (unit in the last place)
    --- both IEEE and MPFR represent 1 as binary(0.1E+1)
    V1: MPFR_Struct renames V(1);
  begin
    if MPFR_Zero_P(V) /= 0 then
      raise MPFR_Exception with "underflow - should not happen";
    end if;
    return LLFloat'Compose(ROne,1-Int(V1.MP_Prec)+Int(V1.MPFR_Exp));
  end ULP_Weight;

  procedure SumErr(S1,S2: in MPFloat; S3: in out MPFloat; Err: out LLFloat) is
    Rounded: Int;
  begin
    Re_Format(S3.Val);
    MPFR_Add(Rounded,S3.Val,S1.Val,S2.Val,GMP_RNDN);
    if Rounded=0 then Err := RZero; else Err := ULP_Weight(S3.Val); end if;
    if Problem then Report("Balls.SumErr"); end if;
  end SumErr;

  procedure DiffErr(S1,S2: in MPFloat; S3: in out MPFloat; Err: out LLFloat) is
    Rounded: Int;
  begin
    Re_Format(S3.Val);
    MPFR_Sub(Rounded,S3.Val,S1.Val,S2.Val,GMP_RNDN);
    if Rounded=0 then Err := RZero; else Err := ULP_Weight(S3.Val); end if;
    if Problem then Report("Balls.DiffErr"); end if;
  end DiffErr;

  --- 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 (MPFR_Zero_P(S.C.Val) /= 0) and then S.R=RZero;
  end IsZero;

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

  procedure SetZero(S: in out Ball) is
  begin
    MPFR_Set_Zero(S.C.Val,0);
    S.R := RZero;
  end SetZero;

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

  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
    V1: MPFR_T renames S1.C.Val;
    V2: MPFR_T renames S2.C.Val;
    Rounded: Int;
  begin
    Re_Format(V1(1).MP_Prec,V2);
    MPFR_Set(Rounded,V2,V1,GMP_RNDN);
    S2.R := S1.R;
  end Copy;

  procedure Swap(S1,S2: in out Ball) is
    R: constant Rad := S1.R;
  begin
    MPFR_Swap(S1.C.Val,S2.C.Val);
    S1.R := S2.R; S2.R := R;
  end Swap;

  --- sets

  function Center0(S: Ball) return Boolean is
  begin
    return (MPFR_Zero_P(S.C.Val) /= 0);
  end Center0;

  function Contains0(S: Ball) return Logical is
    SR: constant Rep := Scal(Flt(S.R));
    SC: constant Rep := Norm(S.C);
  begin
    if Compare(SR,SC)<0 then
      return False;
    else
      return True;
    end if;
  exception
    when Not_Certain => return Uncertain;
  end Contains0;

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

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

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

  procedure ToErr(S: in out Ball) is
  begin
    S.R := Rad(MaxNorm(S.C))+S.R;
    MPFR_Set_Zero(S.C.Val,0);
  end ToErr;

  procedure ToErr(S1: in Ball; S2: in out Ball) is
  begin
    MPFR_Set_Zero(S2.C.Val,0);
    S2.R := Rad(MaxNorm(S1.C))+S1.R;
  end ToErr;

  function ToErr(S: Ball) return Ball is
  begin
    return (CZero,Rad(MaxNorm(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
    Copy(S1.C,S2.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
    MPFR_Set_Zero(S.C.Val,0);
  end ModCenter;

  procedure ModCenter(S1: in Ball; S2: in out Ball) is
  begin
    MPFR_Set_Zero(S2.C.Val,0);
    S2.R := S1.R;
  end ModCenter;

  function ModCenter(S: Ball) return Ball is
  begin
    return (CZero,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
    Eps: Rad;
    RU1,RU2,YL,YU,Y: Rep;
  begin
    if S1.C=S2.C then
      S2.R := Rad'Max(S1.R,S2.R);
    else
      Val(S1.R,RU1);
      Val(S2.R,RU2);
      Min(-(RU1-S1.C),-(RU2-S2.C),YL);
      Max(S1.C+RU1,S2.C+RU2,YU);
      SumErr(YL,YU,Y,Eps);
      Prod(Half,Y,S2.C);
      S2.R := RHalf*(Upper(YU-YL)+Eps);
    end if;
  end Union;

  function Union(S1,S2: Ball) return Ball is
    S: Ball := S2;
  begin
    Union(S1,S);
    return S;
  end Union;

  procedure Intersection(S1: in Ball; S2: in out Ball; Empty: out Logical) is
  begin
    if S1.C=S2.C then
      Empty := False;
      S2.R := Rad'Min(S1.R,S2.R);
    else
      Empty := Uncertain;
      declare
        Eps: Rad;
        RU1,RU2,YL,YU,Y: Rep;
      begin
        Val(S1.R,RU1);
        Val(S2.R,RU2);
        Max(S1.C-RU1,S2.C-RU2,YL);
        Max(-S1.C-RU1,-S2.C-RU2,YU);
        Neg(YU);
        if YL <= YU then
          Empty := False;
        end if;
        Max(-(RU1-S1.C),-(RU2-S2.C),YL);
        Min(S1.C+RU1,S2.C+RU2,YU);
        if YU < YL then
          Empty := True;
        end if;
        SumErr(YL,YU,Y,Eps);
        Prod(Half,Y,S2.C);
        S2.R := RHalf*(Upper(YU-YL)+Eps);
      end;
    end if;
  end Intersection;

  --- order

  procedure RBall(R1,R2: in Rep; S: in out Ball) is
    Eps: Rad;
    Tmp: Rep;
  begin
    if Compare(R2,R1)>0 then
      Diff(R2,R1,Tmp);
    else
      Diff(R1,R2,Tmp);
    end if;
    SumErr(R1,R2,S.C,Eps);
    S.R := RHalf*(Upper(Tmp)+Eps);
    Mult(Half,S.C,Tmp);
  end RBall;

  procedure RSup(S: in Ball; R: in out Rep) is
    Rounded: Int;
    V: MPFR_T renames R.Val;
  begin
    Re_Format(V);
    MPFR_Set_Ld(Rounded,V,Long_Double(S.R),GMP_RNDU);
    MPFR_Add(Rounded,V,S.C.Val,V,GMP_RNDU);
    if Problem then Report("Balls.RSup"); end if;
  end RSup;

  function RSup(S: Ball) return Rep is
    R: Rep;
  begin
    RSup(S,R);
    return R;
  end RSup;

  procedure RInf(S: in Ball; R: in out Rep) is
    Rounded: Int;
    V: MPFR_T renames R.Val;
  begin
    Re_Format(V);
    MPFR_Set_Ld(Rounded,V,Long_Double(S.R),GMP_RNDU);
    MPFR_Sub(Rounded,V,S.C.Val,V,GMP_RNDD);
    if Problem then Report("Balls.RInf"); end if;
  end RInf;

  function RInf(S: Ball) return Rep is
    R: Rep;
  begin
    RInf(S,R);
    return R;
  end RInf;

  function Sign(S: Ball) return Integer is
  begin
    if Sign(RInf(S))>0 then return 1; end if;
    if Sign(RSup(S))<0 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 Compare(RInf(S1),RSup(S2))>0 then return  1; end if;
    if Compare(RSup(S1),RInf(S2))<0 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 RSup(S1) < RInf(S2) then return True; end if;
    if RInf(S1) < RSup(S2) then
      raise Not_Certain;
    end if;
    return False;
  end "<";

  function "<="(S1,S2: Ball) return Boolean is
  begin
    if RSup(S1) <= RInf(S2) then return True; end if;
    if RInf(S1) <= RSup(S2) then
      raise Not_Certain;
    end if;
    return False;
  end "<=";

  function ">"(S1,S2: Ball) return Boolean is
  begin
    if RInf(S1) > RSup(S2) then return True; end if;
    if RSup(S1) > RInf(S2) then
      raise Not_Certain;
    end if;
    return False;
  end ">";

  function ">="(S1,S2: Ball) return Boolean is
  begin
    if RInf(S1) >= RSup(S2) then return True; end if;
    if RSup(S1) >= RInf(S2) then
      raise Not_Certain;
    end if;
    return False;
  end ">=";

  procedure Min(S1: in Ball; S2: in out Ball) is
    RU1,RU2,RL1,RL2: Rep;
  begin
    RInf(S1,RL1);
    RSup(S2,RU2);
    if RL1<RU2 then
      RSup(S1,RU1);
      RInf(S2,RL2);
      if RL2<RU1 then
        RBall(Min(RL1,RL2),Min(RU1,RU2),S2);
      else
        Copy(S1,S2);
      end if;
    end if;
  end Min;

  procedure Min(S1,S2: in Ball; S3: in out Ball) is
  begin
    Copy(S2,S3);
    Min(S1,S3);
  end Min;

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

  procedure Max(S1: in Ball; S2: in out Ball) is
    RU1,RU2,RL1,RL2: Rep;
  begin
    RSup(S1,RU1);
    RInf(S2,RL2);
    if RL2<RU1 then
      RInf(S1,RL1);
      RSup(S2,RU2);
      if RL1<RU2 then
        RBall(Max(RL1,RL2),Max(RU1,RU2),S2);
      else
        Copy(S1,S2);
      end if;
    end if;
  end Max;

  procedure Max(S1,S2: in Ball; S3: in out Ball) is
  begin
    Copy(S2,S3);
    Max(S1,S3);
  end Max;

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

  function Inf(S: Ball) return Flt is
  begin
    return Flt(Lower(RInf(S)));
  end Inf;

  function Sup(S: Ball) return Flt is
  begin
    return Flt(Upper(RSup(S)));
  end Sup;

  --- addition and multiplication

  procedure Neg(S: in out Ball) is
    Rounded: Int;
  begin
    MPFR_Neg(Rounded,S.C.Val,S.C.Val,GMP_RNDN);
  end Neg;

  procedure Neg(S1: in Ball; S2: in out Ball) is
    Rounded: Int;
  begin
    Re_Format(S1.C.Val(1).MP_Prec,S2.C.Val);
    MPFR_Neg(Rounded,S2.C.Val,S1.C.Val,GMP_RNDN);
    S2.R := S1.R;
  end Neg;

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

  procedure Add(I: in Integer; S: in out Ball) is
    Rounded: Int;
    V: MPFR_T renames S.C.Val;
  begin
    Fix_Precision(V,Rounded);
    if Rounded /= 0 then
      S.R := S.R+ULP_Weight(V);
    end if;
    MPFR_Add_Si(Rounded,V,V,Long(I),GMP_RNDN);
    if Rounded /= 0 then
      S.R := S.R+ULP_Weight(V);
    end if;
    if Problem then Report("Balls.Add"); end if;
  end Add;

  procedure Add(S1: in Ball; S2: in out Ball) is
    Rounded: Int;
    V2: MPFR_T renames S2.C.Val;
  begin
    Fix_Precision(V2,Rounded);
    if Rounded /= 0 then
      S2.R := S2.R+ULP_Weight(V2);
    end if;
    MPFR_Add(Rounded,V2,V2,S1.C.Val,GMP_RNDN);
    if Rounded=0 then
      S2.R := S2.R+S1.R;
    else
      S2.R := S2.R+S1.R+ULP_Weight(V2);
    end if;
    if Problem then Report("Balls.Add"); end if;
  end Add;

  procedure Sum(S1,S2: in Ball; S3: in out Ball) is
    Rounded: Int;
    V3: MPFR_T renames S3.C.Val;
  begin
    Re_Format(V3);
    MPFR_Add(Rounded,V3,S1.C.Val,S2.C.Val,GMP_RNDN);
    if Rounded=0 then
      S3.R := S1.R+S2.R;
    else
      S3.R := S1.R+S2.R+ULP_Weight(V3);
    end if;
    if Problem then Report("Sum"); end if;
  end Sum;

  function "+"(S1,S2: Ball) return Ball is
    S3: Ball;
  begin
    Sum(S1,S2,S3);
    return S3;
  end "+";

  procedure Sub(S1: in Ball; S2: in out Ball) is
    Rounded: Int;
    V2: MPFR_T renames S2.C.Val;
  begin
    Fix_Precision(V2,Rounded);
    if Rounded /= 0 then
      S2.R := S2.R+ULP_Weight(V2);
    end if;
    MPFR_Sub(Rounded,V2,V2,S1.C.Val,GMP_RNDN);
    if Rounded=0 then
      S2.R := S2.R+S1.R;
    else
      S2.R := S2.R+S1.R+ULP_Weight(V2);
    end if;
    if Problem then Report("Balls.Sub"); end if;
  end Sub;

  procedure Diff(S1,S2: in Ball; S3: in out Ball) is
    Rounded: Int;
    V3: MPFR_T renames S3.C.Val;
  begin
    Re_Format(V3);
    MPFR_Sub(Rounded,V3,S1.C.Val,S2.C.Val,GMP_RNDN);
    if Rounded=0 then
      S3.R := S1.R+S2.R;
    else
      S3.R := S1.R+S2.R+ULP_Weight(V3);
    end if;
    if Problem then Report("Balls.Diff"); end if;
  end Diff;

  function "-"(S1,S2: Ball) return Ball is
    S3: Ball;
  begin
    Diff(S1,S2,S3);
    return S3;
  end "-";

  procedure Mult(R: in Flt; S,Tmp: in out Ball) is
    Rounded: Int;
    VS: MPFR_T renames S.C.Val;
    VT: MPFR_T renames Tmp.C.Val;
  begin
    Fix_Precision(VS,Rounded);
    if Rounded /= 0 then
      S.R := S.R+ULP_Weight(VS);
    end if;
    Re_Format(VT);
    MPFR_Set_Ld(Rounded,VT,Long_Double(R),GMP_RNDN);
    MPFR_Mul(Rounded,VS,VT,VS,GMP_RNDN);
    if Rounded=0 then
      S.R := abs(Rad(R))*S.R;
    else
      S.R := abs(Rad(R))*S.R+ULP_Weight(VS);
    end if;
    if Problem then Report("Balls.Mult"); end if;
  end Mult;

  procedure Mult(R: in Flt; S: in out Ball) is
    Rounded: Int;
    VS: MPFR_T renames S.C.Val;
    VT: MPFR_T;
  begin
    Fix_Precision(VS,Rounded);
    if Rounded /= 0 then
      S.R := S.R+ULP_Weight(VS);
    end if;
    MPFR_Init2(VT,MP_Precision);
    MPFR_Set_Ld(Rounded,VT,Long_Double(R),GMP_RNDN);
    MPFR_Mul(Rounded,VS,VT,VS,GMP_RNDN);
    MPFR_Clear(VT);
    if Rounded=0 then
      S.R := abs(Rad(R))*S.R;
    else
      S.R := abs(Rad(R))*S.R+ULP_Weight(VS);
    end if;
    if Problem then Report("Balls.Mult"); end if;
  end Mult;

  procedure Prod(R: in Flt; S1: in Ball; S2: in out Ball) is
    Rounded: Int;
    V2: MPFR_T renames S2.C.Val;
  begin
    Re_Format(V2);
    MPFR_Set_Ld(Rounded,V2,Long_Double(R),GMP_RNDN);
    MPFR_Mul(Rounded,V2,S1.C.Val,V2,GMP_RNDN);
    if Rounded=0 then
      S2.R := abs(Rad(R))*S1.R;
    else
      S2.R := abs(Rad(R))*S1.R+ULP_Weight(V2);
    end if;
    if Problem then Report("Balls.Prod"); end if;
  end Prod;

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

  procedure AddProd(R: in Flt; S1: in Ball; S2,Tmp: in out Ball) is
    Rounded: Int;
    V2: MPFR_T renames S2.C.Val;
    VT: MPFR_T renames Tmp.C.Val;
  begin
    Fix_Precision(V2,Rounded);
    if Rounded /= 0 then
      S2.R := S2.R+ULP_Weight(V2);
    end if;
    Reformat(VT);
    MPFR_Set_Ld(Rounded,VT,Long_Double(R),GMP_RNDN);
    MPFR_FMA(Rounded,V2,VT,S1.C.Val,V2,GMP_RNDN);
    if Rounded=0 then
      S2.R := abs(Rad(R))*S1.R+S2.R;
    else
      S2.R := abs(Rad(R))*S1.R+S2.R+ULP_Weight(V2);
    end if;
    if Problem then Report("Balls.AddProd"); end if;
  end AddProd;

  procedure AddProd(R: in Flt; S1: in Ball; S2: in out Ball) is
    Tmp: Ball;
  begin
    AddProd(R,S1,S2,Tmp);
  end AddProd;

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

  procedure Mult(S1: in Ball; S2: in out Ball) is
    Rounded: Int;
    A1: constant Rad := Rad(MaxNorm(S1.C));
    A2: constant Rad := Rad(MaxNorm(S2.C))+S2.R;
    V2: MPFR_T renames S2.C.Val;
  begin
    Fix_Precision(V2,Rounded);
    if Rounded /= 0 then
      S2.R := S2.R+ULP_Weight(V2);
    end if;
    MPFR_Mul(Rounded,V2,V2,S1.C.Val,GMP_RNDN);
    if Rounded=0 then
      S2.R := A1*S2.R+S1.R*A2;
    else
      S2.R := A1*S2.R+S1.R*A2+ULP_Weight(V2);
    end if;
    if Problem then Report("Balls.Mult"); end if;
  end Mult;

  procedure Prod(S1,S2: in Ball; S3: in out Ball) is
    Rounded: Int;
    A1: constant Rad := Rad(MaxNorm(S1.C));
    A2: constant Rad := Rad(MaxNorm(S2.C))+S2.R;
    V3: MPFR_T renames S3.C.Val;
  begin
    Re_Format(V3);
    MPFR_Mul(Rounded,V3,S1.C.Val,S2.C.Val,GMP_RNDN);
    if Rounded=0 then
      S3.R := A1*S2.R+S1.R*A2;
    else
      S3.R := A1*S2.R+S1.R*A2+ULP_Weight(V3);
    end if;
    if Problem then Report("Balls.Prod"); end if;
  end Prod;

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

  procedure AddProd(S1,S2: in Ball; S3,Tmp: in out Ball) is
    --- not using slow MPFR_FMA
    Rounded: Int;
    A1: constant Rad := Rad(MaxNorm(S1.C));
    A2: constant Rad := Rad(MaxNorm(S2.C))+S2.R;
    V3: MPFR_T renames S3.C.Val;
    VT: MPFR_T renames Tmp.C.Val;
  begin
    Re_Format(VT);
    MPFR_Mul(Rounded,VT,S1.C.Val,S2.C.Val,GMP_RNDN);
    if Rounded=0 then
      Tmp.R := A1*S2.R+S1.R*A2;
    else
      Tmp.R := A1*S2.R+S1.R*A2+ULP_Weight(VT);
    end if;
    Fix_Precision(V3,Rounded);
    if Rounded /= 0 then
      S3.R := S3.R+ULP_Weight(V3);
    end if;
    MPFR_Add(Rounded,V3,V3,VT,GMP_RNDN);
    if Rounded=0 then
      S3.R := S3.R+Tmp.R;
    else
      S3.R := S3.R+Tmp.R+ULP_Weight(V3);
    end if;
    if Problem then Report("Balls.AddProd"); end if;
  end AddProd;

  procedure AddProd(S1,S2: in Ball; S3: in out Ball) is
    Rounded: Int;
    A1: constant Rad := Rad(MaxNorm(S1.C));
    A2: constant Rad := Rad(MaxNorm(S2.C))+S2.R;
    V3: MPFR_T renames S3.C.Val;
  begin
    Fix_Precision(V3,Rounded);
    if Rounded /= 0 then
      S3.R := S3.R+ULP_Weight(V3);
    end if;
    MPFR_FMA(Rounded,V3,S1.C.Val,S2.C.Val,V3,GMP_RNDN);  --- slow!
    if Rounded=0 then
      S3.R := A1*S2.R+S1.R*A2+S3.R;
    else
      S3.R := A1*S2.R+S1.R*A2+S3.R+ULP_Weight(V3);
    end if;
    if Problem then Report("Balls.AddProd"); end if;
  end AddProd;

  procedure SumProd(S1,S2,S3: in Ball; S4: in out Ball) is
    Rounded: Int;
    A1: constant Rad := Rad(MaxNorm(S1.C));
    A2: constant Rad := Rad(MaxNorm(S2.C))+S2.R;
    V4: MPFR_T renames S4.C.Val;
  begin
    Re_Format(V4);
    MPFR_FMA(Rounded,V4,S1.C.Val,S2.C.Val,S3.C.Val,GMP_RNDN);
    if Rounded=0 then
      S4.R := A1*S2.R+S1.R*A2+S3.R;
    else
      S4.R := A1*S2.R+S1.R*A2+S3.R+ULP_Weight(V4);
    end if;
    if Problem then Report("Balls.SumProd"); end if;
  end SumProd;

  procedure SubProd(S1,S2: in Ball; S3,Tmp: in out Ball) is
    --- not using slow MPFR_FMA
    Rounded: Int;
    A1: constant Rad := Rad(MaxNorm(S1.C));
    A2: constant Rad := Rad(MaxNorm(S2.C))+S2.R;
    V3: MPFR_T renames S3.C.Val;
    VT: MPFR_T renames Tmp.C.Val;
  begin
    Re_Format(VT);
    MPFR_Mul(Rounded,VT,S1.C.Val,S2.C.Val,GMP_RNDN);
    if Rounded=0 then
      Tmp.R := A1*S2.R+S1.R*A2;
    else
      Tmp.R := A1*S2.R+S1.R*A2+ULP_Weight(VT);
    end if;
    Fix_Precision(V3,Rounded);
    if Rounded /= 0 then
      S3.R := S3.R+ULP_Weight(V3);
    end if;
    MPFR_Sub(Rounded,V3,V3,VT,GMP_RNDN);
    if Problem then Report("Balls.SubProd"); end if;
    if Rounded=0 then
      S3.R := S3.R+Tmp.R;
    else
      S3.R := S3.R+Tmp.R+ULP_Weight(V3);
    end if;
  end SubProd;

  procedure Div(R: in Flt; S,Tmp: in out Ball) is
    Rounded: Int;
    VS: MPFR_T renames S.C.Val;
    VT: MPFR_T renames Tmp.C.Val;
  begin
    Fix_Precision(VS,Rounded);
    if Rounded /= 0 then
      S.R := S.R+ULP_Weight(VS);
    end if;
    Re_Format(VT);
    MPFR_Set_Ld(Rounded,VT,Long_Double(R),GMP_RNDN);
    MPFR_Div(Rounded,VS,VS,VT,GMP_RNDN);
    if Rounded=0 then
      S.R := S.R/abs(Rad(R));
    else
      S.R := S.R/abs(Rad(R))+ULP_Weight(VS);
    end if;
    if Problem then Report("Balls.Div"); end if;
  end Div;

  procedure Div(R: in Flt; S: in out Ball) is
    Rounded: Int;
    VS: MPFR_T renames S.C.Val;
    VT: MPFR_T;
  begin
    Fix_Precision(VS,Rounded);
    if Rounded /= 0 then
      S.R := S.R+ULP_Weight(VS);
    end if;
    MPFR_Init2(VT,MP_Precision);
    MPFR_Set_Ld(Rounded,VT,Long_Double(R),GMP_RNDN);
    MPFR_Div(Rounded,VS,VS,VT,GMP_RNDN);
    MPFR_Clear(VT);
    if Rounded=0 then
      S.R := S.R/abs(Rad(R));
    else
      S.R := S.R/abs(Rad(R))+ULP_Weight(VS);
    end if;
    if Problem then Report("Balls.Div"); end if;
  end Div;

  procedure Quot(S1: in Ball; R: in Flt; S2: in out Ball) is
    Rounded: Int;
    V2: MPFR_T renames S2.C.Val;
  begin
    Re_Format(V2);
    MPFR_Set_Ld(Rounded,V2,Long_Double(R),GMP_RNDN);
    MPFR_Div(Rounded,V2,S1.C.Val,V2,GMP_RNDN);
    if Rounded=0 then
      S2.R := S1.R/abs(Rad(R));
    else
      S2.R := S1.R/abs(Rad(R))+ULP_Weight(V2);
    end if;
    if Problem then Report("Balls.Quot"); end if;
  end Quot;

  function "/"(S: Ball; R: Flt) return Ball is
    T: Ball;
  begin
    Quot(S,R,T);
    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
    Rounded: Int;
    A,X: Rad;
    C1: Rep renames S1.C;
    V2: MPFR_T renames S2.C.Val;
  begin
    if Sign(C1)<0 then
      A := -Upper(C1);
    else
      A := Lower(C1);
    end if;
    X := S1.R-A;
    if X >= RZero then
      raise Undefined with "Balls.Inv error: possible division by zero";
    end if;
    Re_Format(V2);
    MPFR_Set_Si(Rounded,V2,1,GMP_RNDN);
    MPFR_Div(Rounded,V2,V2,C1.Val,GMP_RNDN);
    if Rounded=0 then
      S2.R := (S1.R/A)/(-X);
    else
      S2.R := (S1.R/A)/(-X)+ULP_Weight(V2);
    end if;
    if Problem then Report("Balls.Inv"); end if;
  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
    A: constant Rad := Rad(MaxNorm(S1.C));
  begin
    if S1.R>A then
      S2.R := RHalf*(A+S1.R);
      Val(S2.R,S2.C);
    elsif Sign(S1.C)<0 then
      Neg(S1.C,S2.C);
      S2.R := S1.R;
    else
      Copy(S1.C,S2.C);
      S2.R := S1.R;
    end if;
  end Norm;

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

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

  procedure Sqr(S: in out Ball) is
    Rounded: Int;
    A: constant Rad := Rad(MaxNorm(S.C));
    V: MPFR_T renames S.C.Val;
  begin
    Fix_Precision(V,Rounded);
    if Rounded /= 0 then
      S.R := S.R+ULP_Weight(V);
    end if;
    MPFR_Mul(Rounded,V,V,V,GMP_RNDN);
    if Rounded=0 then
      S.R := (A+A+S.R)*S.R;
    else
      S.R := (A+A+S.R)*S.R+ULP_Weight(V);
    end if;
    if Problem then Report("Balls.Sqr"); end if;
  end Sqr;

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

  procedure Sqrt(S1: in Ball; S2: in out Ball) is
    --- assumig correct rounding for Sqrt
    use LLFloat_EF;
  begin
    if IsZero(S1) then
      SetZero(S2);
    else
      declare
        C1: Rep renames S1.C;
        C2: Rep renames S2.C;
      begin
        Val(S1.R,C2);
        if C1=C2 then
          S2.R := Sqrt(RHalf*S1.R);
          C2 := Val(S2.R);
        else
          RBall(-((C2-C1)/Sqrt(C1-C2)),Sqrt(C1+C2),S2);
        end if;
      end;
    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
    DX,YL,YU: Rep;
    DY: Rep renames DX;
    Tmp: Rep renames S2.C;
    Eps: Rad;
  begin
    Val(S1.R,DX);
    Root(K,S1.C,DX,YL,YU);
    Mult(Half,YU,Tmp);
    Mult(NegHalf,YL,Tmp);
    Sum(YU,YL,DY);
    DiffErr(YU,YL,S2.C,Eps);
    S2.R := Upper(DY)+Eps;
  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
    DX,YL,YU: Rep;
    DY: Rep renames DX;
    Tmp: Rep renames S2.C;
    Eps: Rad;
  begin
    Val(S1.R,DX);
    Exp(S1.C,DX,YL,YU);
    Mult(Half,YU,Tmp);
    Mult(NegHalf,YL,Tmp);
    Sum(YU,YL,DY);
    DiffErr(YU,YL,S2.C,Eps);
    S2.R := Upper(DY)+Eps;
  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
    DX,YL,YU: Rep;
    DY: Rep renames DX;
    Tmp: Rep renames S2.C;
    Eps: Rad;
  begin
    Val(S1.R,DX);
    Log(S1.C,DX,YL,YU);
    Mult(Half,YU,Tmp);
    Mult(NegHalf,YL,Tmp);
    Sum(YU,YL,DY);
    DiffErr(YU,YL,S2.C,Eps);
    S2.R := Upper(DY)+Eps;
  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
    DX,YL,YU: Rep;
    DY: Rep renames DX;
    Tmp: Rep renames S2.C;
    Eps: Rad;
  begin
    Val(S1.R,DX);
    ArcCos(S1.C,DX,YL,YU);
    Mult(Half,YU,Tmp);
    Mult(NegHalf,YL,Tmp);
    Sum(YU,YL,DY);
    DiffErr(YU,YL,S2.C,Eps);
    S2.R := Upper(DY)+Eps;
  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
    DX,YL,YU: Rep;
    DY: Rep renames DX;
    Tmp: Rep renames S2.C;
    Eps: Rad;
  begin
    Val(S1.R,DX);
    ArcSin(S1.C,DX,YL,YU);
    Mult(Half,YU,Tmp);
    Mult(NegHalf,YL,Tmp);
    Sum(YU,YL,DY);
    DiffErr(YU,YL,S2.C,Eps);
    S2.R := Upper(DY)+Eps;
  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,DY: Rep;
    Tmp: Rep renames S2.C;
    Eps: Rad;
  begin
    Cos(S1.C,YL,YU);
    Mult(Half,YU,Tmp);
    Mult(NegHalf,YL,Tmp);
    Sum(YU,YL,DY);
    DiffErr(YU,YL,S2.C,Eps);
    S2.R := Upper(DY)+Eps+S1.R;
  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,DY: Rep;
    Tmp: Rep renames S2.C;
    Eps: Rad;
  begin
    Sin(S1.C,YL,YU);
    Mult(Half,YU,Tmp);
    Mult(NegHalf,YL,Tmp);
    Sum(YU,YL,DY);
    DiffErr(YU,YL,S2.C,Eps);
    S2.R := Upper(DY)+Eps+S1.R;
  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)
    Eps: Rad;
    YL,YU,DY: Rep;
    S: Ball;
  begin
    ArcCos(CZero,CZero,YL,YU);
    Diff(YU,YL,DY);
    SumErr(YL,YU,S.C,Eps);
    S.R := Upper(DY)+Eps;
    return S;
  end Pi;

  --- conversion and i/o

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

  procedure Assign(Q: in Rational; S: in out Ball) is
    D: constant LInt := LDen(Q);
  begin
    Assign(Flt(LNum(Q)),S.C);
    S.R := RZero;
    if D /= 1 then Div(Flt(D),S); end if;
  end Assign;

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

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

  function Scal(Q: Rational) return Ball is
    D: constant LInt := LDen(Q);
    S: Ball;
  begin
    Assign(Flt(LNum(Q)),S.C);
    S.R := RZero;
    if D /= 1 then Div(Flt(D),S); end if;
    return S;
  end Scal;

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

  procedure Enclose(R1,R2: in Flt; S: in out Ball) is
    H1: constant Rad := Rad(Half*R1);
    H2: constant Rad := Rad(Half*R2);
    C: constant Rad := H1+H2;
    NegC: constant Rad := -H1-H2;
  begin
    Val(C,S.C);
    S.R := abs(H2-H1)+(C+NegC);
  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 Approx(S.C);
  end Approx;

  procedure Show1(N: in String; S: in Ball; NewLine: in Boolean := True) is
    Digs: constant Natural := 0; -- use default
    use LLFloat_FS;
  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
    use LLFloat_FS;
  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
    use LLFLoat_FS;
    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 := ULP_Weight(S.C.Val);
      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
    PutStd(S.C,C,R);
    R := S.R+R;
  end PutStd;

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

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

  procedure GetNum(R: in Rep; S: in out Ball) is
  begin
    Copy(R,S.C);
    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
    if not Check_Flags then
      MPFR_Clear_Flags;
      Check_Flags := True;
    end if;
    Rounding.FPU.Set_Rounding_Mode((Up,Strict));
  end Proper_Rounding;

  procedure Proper_Rounding is
  begin
    if not Check_Flags then
      MPFR_Clear_Flags;
      Check_Flags := True;
    end if;
    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;

begin

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

end MPFR.Floats.Balls;
