with Strings, Ints, QPowers, Roots;
use Strings, Ints;

pragma Elaborate_All (Strings,Ints,QPowers,Roots);

package body MPFR.Floats is

  package Txt_IO renames Ada.Text_IO;

  procedure Initialize(S: in out MPFloat) is
    --- used during declarations that do not assign values
  begin
    MPFR_Init2(S.Val,D_Precision);
  end Initialize;

  --- An assignment dest := source triggers
  --- (1) finalize(dest)
  --- (2) bitwise copy of source into dest
  --- (3) adjust(dest)

  procedure Adjust(S: in out MPFloat) is
    --- also adjusts precision
    Rounded: Int;
    V1: MPFR_T renames S.Val;
    V0: constant MPFR_T := V1;
  begin
    if MPFR_NaN_P(V0)=0 then
      MPFR_Init2(V1,V0(1).MP_Prec);
      MPFR_Set(Rounded,V1,V0,GMP_RNDN);
    else
      MPFR_Init2(V1,D_Precision); --- copy NaN
    end if;
  end Adjust;

  procedure Finalize(S: in out MPFloat) is
    --- also used during assignments, to finalize previous content
  begin
    MPFR_Clear(S.Val);
  end Finalize;

  procedure Reformat(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 Reformat;

  procedure Reformat(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 Reformat;

  procedure Fix_Precision(V: in out MPFR_T) is
    --- adjust precision if necessary, but keep (rounded) value
    --- use GMP_RNDN
  begin
    if V(1).MP_Prec /= MP_Precision then
      declare
        Rounded: Int;
        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 Num_Problem return Boolean is
    --- ignore underflow
  begin
    if MPFR_Overflow_P>0 or MPFR_NaNFlag_P>0 or MPFR_ERangeFlag_P>0 then
      if MPFR_Overflow_P>0   then Show0("Overflow ",False);  end if;
      if MPFR_NaNFlag_P>0    then Show0("NaN ",False);       end if;
      if MPFR_ERangeFlag_P>0 then Show0("ERange ",False);    end if;
      return True;
    else
      return False;
    end if;
  end Num_Problem;

  function Problem return Boolean is
  begin
    if MPFR_Underflow_P>0 or MPFR_Overflow_P>0 or MPFR_NaNFlag_P>0 or MPFR_ERangeFlag_P>0 then
      if MPFR_Underflow_P>0  then Show0("Underflow ",False); end if;
      if MPFR_Overflow_P>0   then Show0("Overflow ",False);  end if;
      if MPFR_NaNFlag_P>0    then Show0("NaN ",False);       end if;
      if MPFR_ERangeFlag_P>0 then Show0("ERange ",False);    end if;
      return True;
    else
      return False;
    end if;
  end Problem;

  procedure Report(Where: in String) is
  begin
    Show0("error occured in " & Where,False);
    raise MPFR_Exception;
  end Report;

  procedure Clear_Flags is
  begin
    if Check_Flags and then Problem then Report("Clear_Flags"); end if;
    MPFR_Clear_Flags;
  end Clear_Flags;

  function Raised_Inexact return Boolean is
    I: constant Int := MPFR_InexFlag_P;
  begin
    Clear_Flags;
    return (I>0);
  end Raised_Inexact;

  ---  auxiliary

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

  procedure NOP(Dummy1: in Flt; Dummy2: in out MPFloat) is
    pragma Unreferenced(Dummy1,Dummy2);
  begin
    null;
  end NOP;

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

  function IsZero(S: MPFloat) return Boolean is
  begin
    return (MPFR_Zero_P(S.Val) /= 0);
  end IsZero;

  procedure SetZero(S: in out MPFloat) is
  begin
    MPFR_Set_Zero(S.Val,0);
  end SetZero;

  procedure SetZero(Dummy: in Flt; S: in out MPFloat) is
    pragma Unreferenced(Dummy);
  begin
    MPFR_Set_Zero(S.Val,0);
  end SetZero;

  procedure SetZero(Dummy: in MPFloat; S: in out MPFloat) is
    pragma Unreferenced(Dummy);
  begin
    MPFR_Set_Zero(S.Val,0);
  end SetZero;

  function ReturnZero(Dummy: MPFloat) return Flt is
    pragma Unreferenced(Dummy);
  begin
    return Zero;
  end ReturnZero;

  function ReturnZero(Dummy: Flt) return MPFloat is
    pragma Unreferenced(Dummy);
    V: MPFR_T;
  begin
    MPFR_Init2(V,Limb_Size);
    MPFR_Set_Zero(V,0);
    return (Controlled with V);
  end ReturnZero;

  function ReturnZero(Dummy: MPFloat) return MPFloat is
    pragma Unreferenced(Dummy);
    V: MPFR_T;
  begin
    MPFR_Init2(V,Limb_Size);
    MPFR_Set_Zero(V,0);
    return (Controlled with V);
  end ReturnZero;

  procedure Copy(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S1.Val(1).MP_Prec,S2.Val);
    MPFR_Set(Rounded,S2.Val,S1.Val,GMP_RNDN);
  end Copy;

  function Copy(S: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,S.Val(1).MP_Prec);
    MPFR_Set(Rounded,V,S.Val,GMP_RNDN);
    return (Controlled with V);
  end Copy;

  function Approx(S: MPFloat) return Flt is
    A: constant Long_Double := MPFR_Get_Ld(S.Val,GMP_RNDN);
  begin
    if Check_Flags and then Num_Problem then Report("Approx"); end if;
    return Flt(A);
  end Approx;

  function FirstArg(Q: Rational; Dummy: Flt) return MPFloat is
    pragma Unreferenced(Dummy);
    D: constant LInt := LDen(Q);
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Set_Si(Rounded,V,Long(LNum(Q)),GMP_RNDN);
    if D /= 1 then
      declare
        U: MPFR_T;
      begin
        MPFR_Init2(U,MP_Precision);
        MPFR_Set_Si(Rounded,U,Long(D),GMP_RNDN);
        MPFR_Div(Rounded,V,V,U,GMP_RNDN);
        MPFR_Clear(U);
      end;
    end if;
    if Check_Flags and then Num_Problem then Report("FirstArg"); end if;
    return (Controlled with V);
  end FirstArg;

  function FirstArg(S: MPFloat; Dummy: Flt) return MPFloat is
    pragma Unreferenced(Dummy);
    Rounded: Int;
    VS: MPFR_T renames S.Val;
    VT: MPFR_T;
  begin
    MPFR_Init2(VT,VS(1).MP_Prec);
    MPFR_Set(Rounded,VT,VS,GMP_RNDN);
    return (Controlled with VT);
  end FirstArg;

  procedure Between(S1,S2: in Flt; S3: in out MPFloat) is
    --- approx Half*(S1+S2)
    Rounded: Int;
  begin
    Reformat(LD_Precision,S3.Val);
    MPFR_Set_Ld(Rounded,S3.Val,Long_Double(Half*(S1+S2)),GMP_RNDN);
  end Between;

  function Between(S1,S2: Flt) return MPFloat is
    --- approx Half*(S1+S2)
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,LD_Precision);
    MPFR_Set_Ld(Rounded,V,Long_Double(Half*(S1+S2)),GMP_RNDN);
    return (Controlled with V);
  end Between;

  procedure Between(S1: in MPFloat; S2: in out MPFloat) is
    --- approx Half*(S1+S2)
    Rounded: Int;
    V2: MPFR_T renames S2.Val;
  begin
    Fix_Precision(V2);
    MPFR_Add(Rounded,V2,S1.Val,S2.Val,GMP_RNDN);
    MPFR_Mul_D(Rounded,V2,V2,Double(Half),GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Between"); end if;
  end Between;

  function Between(S1,S2: MPFloat) return MPFloat is
    --- approx Half*(S1+S2)
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Add(Rounded,V,S1.Val,S2.Val,GMP_RNDN);
    MPFR_Mul_D(Rounded,V,V,Double(Half),GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Between"); end if;
    return (Controlled with V);
  end Between;

  --- basic

  function Info(Dummy: MPFloat) return Scalar_Info is
    pragma Unreferenced(Dummy);
    I: Scalar_Info;
  begin
    I.RepType    := "MPFR Float";
    I.IsStandard := False;          --- rest is default
    return I;
  end Info;

  procedure Details(S: MPFloat) is
    V: MPFR_T renames S.Val;
    V1: MPFR_Struct renames V(1);
    ULP: constant Flt := Flt'Compose(1.0,1-Int(V1.MP_Prec)+Int(V1.MPFR_Exp));
  begin
    Show1("MP_Prec   : ",Integer(V1.MP_Prec));
    Show1("MPFR_Sign : ",Integer(V1.MPFR_Sign));
    Show1("MP_Exp    : ",Integer(V1.MPFR_Exp));
    Show0("Dec Value : " & Image(S,10));
    Show0("Hex Value : " & Image(S,16));
    Show0("Bin Value : " & Image(S,2));
    Show0("Value-ULP : " & Image(S-Scal(ULP),2));
    Show0("Value+ULP : " & Image(S+Scal(ULP),2));
    Show0;
  end Details;

  procedure SetSmallZero(S: in out MPFloat) is
    V: MPFR_T renames S.Val;
  begin
    if V(1).MP_Prec > Limb_Size then
      MPFR_Set_Prec(V,Limb_Size);
    end if;
    MPFR_Set_Zero(V,0);
  end SetSmallZero;

  procedure Swap(S1,S2: in out MPFloat) is
  begin
    MPFR_Swap(S1.Val,S2.Val);
  end Swap;

  --- sets

  function Contains0(S: MPFloat) return Logical is
  begin
    Clear_Flags;
    declare
      I: constant Int := MPFR_Zero_P(S.Val);
    begin
      if Num_Problem then
        MPFR_Clear_Flags;
        return Uncertain;
      elsif I=0 then
        return False;
      else
        return True;
      end if;
    end;
  end Contains0;

  function Contains(S1,S2: MPFloat) return Logical is
  begin
    Clear_Flags;
    declare
      I: constant Int := MPFR_Cmp(S1.Val,S2.Val);
    begin
      if Num_Problem then
        MPFR_Clear_Flags;
        return Uncertain;
      elsif I=0 then
        return True;
      else
        return False;
      end if;
    end;
  end Contains;

  procedure Intersection(S1: in MPFloat; S2: in out MPFloat; Empty: out Logical) is
  begin
    Clear_Flags;
    declare
      I: constant Int := MPFR_Cmp(S1.Val,S2.Val);
    begin
      if Num_Problem then
        MPFR_Clear_Flags;
        Empty := Uncertain;
      elsif I=0 then
        Empty := False;
      else
        Between(S1,S2);
        Empty := True;
      end if;
    end;
  end Intersection;

  --- order

  function Sign(S: MPFloat) return Integer is
  begin
    return Integer(MPFR_Sgn(S.Val));
  end Sign;

  function Compare(R: Flt; S: MPFloat) return Integer is
  begin
    Clear_Flags;
    declare
      I: constant Int := MPFR_Cmp_Ld(S.Val,Long_Double(R));
    begin
      if Problem then
        MPFR_Clear_Flags;
        raise Not_Certain;
      end if;
      return Integer(I);
    end;
  end Compare;

  function Compare(S1,S2: MPFloat) return Integer is
  begin
    Clear_Flags;
    declare
      I: constant Int := MPFR_Cmp(S1.Val,S2.Val);
    begin
      if Problem then
        MPFR_Clear_Flags;
        raise Not_Certain;
      end if;
      return Integer(I);
    end;
  end Compare;

  function "<"(S1,S2: MPFloat) return Boolean is
  begin
    return (Compare(S1,S2) < 0);
  end "<";

  function "<="(S1,S2: MPFloat) return Boolean is
  begin
    return (Compare(S1,S2) <= 0);
  end "<=";

  function "="(S1,S2: MPFloat) return Boolean is
  begin
    return (Compare(S1,S2) = 0);
  end "=";

  function ">="(S1,S2: MPFloat) return Boolean is
  begin
    return (Compare(S1,S2) >= 0);
  end ">=";

  function ">"(S1,S2: MPFloat) return Boolean is
  begin
    return (Compare(S1,S2) > 0);
  end ">";

  procedure Min(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    if MPFR_Cmp(S1.Val,S2.Val) < 0 then
      Reformat(S1.Val(1).MP_Prec,S2.Val);
      MPFR_Set(Rounded,S2.Val,S1.Val,GMP_RNDN);
      if Check_Flags and then Problem then Report("Min"); end if;
    end if;
  end Min;

  procedure Min(S1,S2: in MPFloat; S3: in out MPFloat) is
    Rounded: Int;
  begin
    if MPFR_Cmp(S1.Val,S2.Val) < 0 then
      Reformat(S1.Val(1).MP_Prec,S3.Val);
      MPFR_Set(Rounded,S3.Val,S1.Val,GMP_RNDN);
    else
      Reformat(S2.Val(1).MP_Prec,S3.Val);
      MPFR_Set(Rounded,S3.Val,S2.Val,GMP_RNDN);
    end if;
    if Check_Flags and then Problem then Report("Min"); end if;
  end Min;

  function Min(S1,S2: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    if MPFR_Cmp(S1.Val,S2.Val) < 0 then
      MPFR_Init2(V,S1.Val(1).MP_Prec);
      MPFR_Set(Rounded,V,S1.Val,GMP_RNDN);
    else
      MPFR_Init2(V,S2.Val(1).MP_Prec);
      MPFR_Set(Rounded,V,S2.Val,GMP_RNDN);
    end if;
    if Check_Flags and then Problem then Report("Min"); end if;
    return (Controlled with V);
  end Min;

  procedure Max(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    if MPFR_Cmp(S1.Val,S2.Val) > 0 then
      Reformat(S1.Val(1).MP_Prec,S2.Val);
      MPFR_Set(Rounded,S2.Val,S1.Val,GMP_RNDN);
      if Check_Flags and then Problem then Report("Max"); end if;
    end if;
  end Max;

  procedure Max(S1,S2: in MPFloat; S3: in out MPFloat) is
    Rounded: Int;
  begin
    if MPFR_Cmp(S1.Val,S2.Val) > 0 then
      Reformat(S1.Val(1).MP_Prec,S3.Val);
      MPFR_Set(Rounded,S3.Val,S1.Val,GMP_RNDN);
    else
      Reformat(S2.Val(1).MP_Prec,S3.Val);
      MPFR_Set(Rounded,S3.Val,S2.Val,GMP_RNDN);
    end if;
    if Check_Flags and then Problem then Report("Min"); end if;
  end Max;

  function Max(S1,S2: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    if MPFR_Cmp(S1.Val,S2.Val) > 0 then
      MPFR_Init2(V,S1.Val(1).MP_Prec);
      MPFR_Set(Rounded,V,S1.Val,GMP_RNDN);
    else
      MPFR_Init2(V,S2.Val(1).MP_Prec);
      MPFR_Set(Rounded,V,S2.Val,GMP_RNDN);
    end if;
    if Check_Flags and then Problem then Report("Max"); end if;
    return (Controlled with V);
  end Max;

  function Sup(S: MPFloat) return Flt is
    A: constant Long_Double := MPFR_Get_Ld(S.Val,GMP_RNDU);
  begin
    if Check_Flags and then Problem then Report("Sup"); end if;
    return Flt(A);
  end Sup;

  function Inf(S: MPFloat) return Flt is
    A: constant Long_Double := MPFR_Get_Ld(S.Val,GMP_RNDD);
  begin
    if Check_Flags and then Problem then Report("Inf"); end if;
    return Flt(A);
  end Inf;

  --- arithmetic

  function "-"(S: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,S.Val(1).MP_Prec);
    MPFR_Neg(Rounded,V,S.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report(" -"); end if;
    return (Controlled with V);
  end "-";

  procedure Neg(S: in out MPFloat) is
    Rounded: Int;
  begin
    MPFR_Neg(Rounded,S.Val,S.Val,GMP_RNDN); --- just adjusts sign
  end Neg;

  procedure Neg(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S1.Val(1).MP_Prec,S2.Val);
    MPFR_Neg(Rounded,S2.Val,S1.Val,GMP_RNDN);
  end Neg;

  procedure Add(I: in Integer; S: in out MPFloat) is
    Rounded: Int;
  begin
    Fix_Precision(S.Val);
    MPFR_Add_Si(Rounded,S.Val,S.Val,Long(I),GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Add"); end if;
  end Add;

  function "+"(S1,S2: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Add(Rounded,V,S1.Val,S2.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("+"); end if;
    return (Controlled with V);
  end "+";

  procedure Add(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Fix_Precision(S2.Val);
    MPFR_Add(Rounded,S2.Val,S2.Val,S1.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Add"); end if;
  end Add;

  procedure Sum(S1,S2: in MPFloat; S3: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S3.Val);
    MPFR_Add(Rounded,S3.Val,S1.Val,S2.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Sum"); end if;
  end Sum;

  function "-"(S1,S2: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Sub(Rounded,V,S1.Val,S2.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("-"); end if;
    return (Controlled with V);
  end "-";

  procedure Sub(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Fix_Precision(S2.Val);
    MPFR_Sub(Rounded,S2.Val,S2.Val,S1.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Sub"); end if;
  end Sub;

  procedure Diff(S1,S2: in MPFloat; S3: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S3.Val);
    MPFR_Sub(Rounded,S3.Val,S1.Val,S2.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Diff"); end if;
  end Diff;

  function "*"(R: Flt; S: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Set_Ld(Rounded,V,Long_Double(R),GMP_RNDN);
    MPFR_Mul(Rounded,V,S.Val,V,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("*"); end if;
    return (Controlled with V);
  end "*";

  procedure Mult(R: in Flt; S,Tmp: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(Tmp.Val);
    MPFR_Set_Ld(Rounded,Tmp.Val,Long_Double(R),GMP_RNDN);
    Fix_Precision(S.Val);
    MPFR_Mul(Rounded,S.Val,Tmp.Val,S.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Mult"); end if;
  end Mult;

  procedure Mult(R: in Flt; S: in out MPFloat) is
    Rounded: Int;
  begin
    if R=Zero then
      MPFR_Set_Zero(S.Val,0);
    elsif R=One then
      null;
    elsif R=NegOne then
      MPFR_Neg(Rounded,S.Val,S.Val,GMP_RNDN); --- just adjusts sign
    else
      declare
        Tmp: MPFloat;
      begin
        Mult(R,S,Tmp);
      end;
    end if;
  end Mult;

  procedure Prod(R: in Flt; S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S2.Val);
    MPFR_Set_Ld(Rounded,S2.Val,Long_Double(R),GMP_RNDN);
    MPFR_Mul(Rounded,S2.Val,S1.Val,S2.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Prod"); end if;
  end Prod;

  procedure AddProd(R: in Flt; S1: in MPFloat; S2,Tmp: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(Tmp.Val);
    MPFR_Set_Ld(Rounded,Tmp.Val,Long_Double(R),GMP_RNDN);
    Fix_Precision(S2.Val);
    MPFR_FMA(Rounded,S2.Val,Tmp.Val,S1.Val,S2.Val,GMP_RNDN); --- slow!
    if Check_Flags and then Num_Problem then Report("AddProd"); end if;
  end AddProd;

  procedure AddProd(R: in Flt; S1: in MPFloat; S2: in out MPFloat) is
  begin
    if R=Zero then
      null;
    elsif R=One then
      Add(S1,S2);
    elsif R=NegOne then
      Sub(S1,S2);
    else
      declare
        Tmp: MPFloat;
      begin
        AddProd(R,S1,S2,Tmp);
      end;
    end if;
  end AddProd;

  procedure Mult(Q: in Rational; S: in out MPFloat) is
    Tmp: MPFloat;
  begin
    Mult(Flt(LNum(Q)),S,Tmp);
    Div(Flt(LDen(Q)),S,Tmp);
  end Mult;

  function "*"(S1,S2: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Mul(Rounded,V,S1.Val,S2.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("*"); end if;
    return (Controlled with V);
  end "*";

  procedure Mult(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Fix_Precision(S2.Val);
    MPFR_Mul(Rounded,S2.Val,S2.Val,S1.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Mult"); end if;
  end Mult;

  procedure Prod(S1,S2: in MPFloat; S3: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S3.Val);
    MPFR_Mul(Rounded,S3.Val,S1.Val,S2.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Prod"); end if;
  end Prod;

  procedure AddProd(S1,S2: in MPFloat; S3,Tmp: in out MPFloat) is
    --- not using slow MPFR_FMA
    Rounded: Int;
  begin
    Reformat(Tmp.Val);
    MPFR_Mul(Rounded,Tmp.Val,S1.Val,S2.Val,GMP_RNDN);
    Fix_Precision(S3.Val);
    MPFR_Add(Rounded,S3.Val,S3.Val,Tmp.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("AddProd"); end if;
  end AddProd;

  procedure AddProd(S1,S2: in MPFloat; S3: in out MPFloat) is
    Rounded: Int;
  begin
    Fix_Precision(S3.Val);
    MPFR_FMA(Rounded,S3.Val,S1.Val,S2.Val,S3.Val,GMP_RNDN);  --- slow!
    if Check_Flags and then Num_Problem then Report("AddProd"); end if;
  end AddProd;

  procedure SumProd(S1,S2,S3: in MPFloat; S4: in out MPFloat) is
    --- not using slow MPFR_FMA
    Rounded: Int;
  begin
    Reformat(S4.Val);
    MPFR_Mul(Rounded,S4.Val,S1.Val,S2.Val,GMP_RNDN);
    MPFR_Add(Rounded,S4.Val,S4.Val,S3.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("SumProd"); end if;
  end SumProd;

  procedure SubProd(S1,S2: in MPFloat; S3,Tmp: in out MPFloat) is
    --- not using slow MPFR_FMA
    Rounded: Int;
  begin
    Reformat(Tmp.Val);
    if GMP_RNDN=GMP_RNDU then
      MPFR_Mul(Rounded,Tmp.Val,S1.Val,S2.Val,GMP_RNDD);
    elsif GMP_RNDN=GMP_RNDD then
      MPFR_Mul(Rounded,Tmp.Val,S1.Val,S2.Val,GMP_RNDU);
    else
      MPFR_Mul(Rounded,Tmp.Val,S1.Val,S2.Val,GMP_RNDN);
    end if;
    Fix_Precision(S3.Val);
    MPFR_Sub(Rounded,S3.Val,S3.Val,Tmp.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("SubProd"); end if;
  end SubProd;

  procedure Div(R: in Flt; S,Tmp: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(Tmp.Val);
    MPFR_Set_Ld(Rounded,Tmp.Val,Long_Double(R),GMP_RNDN);
    Fix_Precision(S.Val);
    MPFR_Div(Rounded,S.Val,S.Val,Tmp.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Div"); end if;
  end Div;

  procedure Div(R: in Flt; S: in out MPFloat) is
    Tmp: MPFloat;
  begin
    Div(R,S,Tmp);
  end Div;

  procedure Quot(S1: in MPFloat; R: in Flt; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S2.Val);
    MPFR_Set_Ld(Rounded,S2.Val,Long_Double(R),GMP_RNDN);
    MPFR_Div(Rounded,S2.Val,S1.Val,S2.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Quot"); end if;
  end Quot;

  function "/"(S: MPFloat; R: Flt) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Set_Ld(Rounded,V,Long_Double(R),GMP_RNDN);
    MPFR_Div(Rounded,V,S.Val,V,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("/"); end if;
    return (Controlled with V);
  end "/";

  procedure Div(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Fix_Precision(S2.Val);
    MPFR_Div(Rounded,S2.Val,S2.Val,S1.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Div"); end if;
  end Div;

  procedure Quot(S1,S2: in MPFloat; S3: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S3.Val);
    MPFR_Div(Rounded,S3.Val,S1.Val,S2.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Quot"); end if;
  end Quot;

  function "/"(S1,S2: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Div(Rounded,V,S1.Val,S2.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("/"); end if;
    return (Controlled with V);
  end "/";

  procedure Inv(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S2.Val);
    MPFR_Set_Si(Rounded,S2.Val,1,GMP_RNDN);
    MPFR_Div(Rounded,S2.Val,S2.Val,S1.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Inv"); end if;
  end Inv;

  function Inv(S: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Set_Si(Rounded,V,1,GMP_RNDN);
    MPFR_Div(Rounded,V,V,S.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Inv"); end if;
    return (Controlled with V);
  end Inv;

  procedure IPower(I: Integer; S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S2.Val);
    MPFR_Pow_Si(Rounded,S2.Val,S1.Val,Long(I),GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("IPower"); end if;
  end IPower;

  procedure IPower(I: Integer; S1: in MPFloat; S2,Dummy: in out MPFloat) is
    pragma Unreferenced(Dummy);
    Rounded: Int;
  begin
    Reformat(S2.Val);
    MPFR_Pow_Si(Rounded,S2.Val,S1.Val,Long(I),GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("IPower"); end if;
  end IPower;

  function "**"(S: MPFloat; I: Integer) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Pow_Si(Rounded,V,S.Val,Long(I),GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("**"); end if;
    return (Controlled with V);
  end "**";

  --- fun

  function "abs"(S: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,S.Val(1).MP_Prec);
    MPFR_Abs(Rounded,V,S.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("abs"); end if;
    return (Controlled with V);
  end "abs";

  procedure Norm(S: in out MPFloat) is
    Rounded: Int;
  begin
    if MPFR_Sgn(S.Val)<0 then
      MPFR_Neg(Rounded,S.Val,S.Val,GMP_RNDN); --- just adjusts sign
    end if;
  end Norm;

  procedure Norm(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S1.Val(1).MP_Prec,S2.Val);
    MPFR_Abs(Rounded,S2.Val,S1.Val,GMP_RNDN);
    if Check_Flags and then Problem then Report("Norm"); end if;
  end Norm;

  function MaxNorm(S: MPFloat) return Radius is
    I: constant Int := MPFR_Sgn(S.Val);
    N: Long_Double;
  begin
    if I<0 then
      N := -MPFR_Get_Ld(S.Val,GMP_RNDD);
    elsif I>0 then
      N := MPFR_Get_Ld(S.Val,GMP_RNDU);
    else
      N := 0.0;
    end if;
    if Check_Flags and then Problem then Report("MaxNorm"); end if;
    return Flt(N);
  end MaxNorm;

  procedure Sqr(S: in out MPFloat) is
    Rounded: Int;
  begin
    Fix_Precision(S.Val);
    MPFR_Sqr(Rounded,S.Val,S.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Sqr"); end if;
  end Sqr;

  procedure Sqr(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S2.Val);
    MPFR_Sqr(Rounded,S2.Val,S1.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Sqr"); end if;
  end Sqr;

  function Sqr(S: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Sqr(Rounded,V,S.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Sqr"); end if;
    return (Controlled with V);
  end Sqr;

  procedure Sqrt(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S2.Val);
    MPFR_Sqrt(Rounded,S2.Val,S1.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Sqrt"); end if;
  end Sqrt;

  function Sqrt(S: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Sqrt(Rounded,V,S.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Sqrt"); end if;
    return (Controlled with V);
  end Sqrt;

  procedure Root(K: Positive; S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S2.Val);
    MPFR_Root(Rounded,S2.Val,S1.Val,Unsigned_Long(K),GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Root"); end if;
  end Root;

  function Root(K: Positive; S: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Root(Rounded,V,S.Val,Unsigned_Long(K),GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Root"); end if;
    return (Controlled with V);
  end Root;

  procedure Exp(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S2.Val);
    MPFR_Exp(Rounded,S2.Val,S1.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Exp"); end if;
  end Exp;

  function Exp(S: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Exp(Rounded,V,S.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Exp"); end if;
    return (Controlled with V);
  end Exp;

  procedure Log(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S2.Val);
    MPFR_Log(Rounded,S2.Val,S1.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Log"); end if;
  end Log;

  function Log(S: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Log(Rounded,V,S.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Log"); end if;
    return (Controlled with V);
  end Log;

  procedure ArcCos(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S2.Val);
    MPFR_ACos(Rounded,S2.Val,S1.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("ArcCos"); end if;
  end ArcCos;

  function ArcCos(S: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_ACos(Rounded,V,S.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("ArcCos"); end if;
    return (Controlled with V);
  end ArcCos;

  procedure ArcSin(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S2.Val);
    MPFR_ASin(Rounded,S2.Val,S1.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("ArcSin"); end if;
  end ArcSin;

  function ArcSin(S: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_ASin(Rounded,V,S.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("ArcSin"); end if;
    return (Controlled with V);
  end ArcSin;

  procedure Cos(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S2.Val);
    MPFR_Cos(Rounded,S2.Val,S1.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Cos"); end if;
  end Cos;

  function Cos(S: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Cos(Rounded,V,S.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Cos"); end if;
    return (Controlled with V);
  end Cos;

  procedure Sin(S1: in MPFloat; S2: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(S2.Val);
    MPFR_Sin(Rounded,S2.Val,S1.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Sin"); end if;
  end Sin;

  function Sin(S: MPFloat) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Sin(Rounded,V,S.Val,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Sin"); end if;
    return (Controlled with V);
  end Sin;

  procedure Simple_Random(S: in out MPFloat) is
  begin
    Assign(Simple_Random,S);
  end Simple_Random;

  function Pi return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Const_Pi(Rounded,V,GMP_RNDN);
    if Check_Flags and then Num_Problem then Report("Pi"); end if;
    return (Controlled with V);
  end Pi;

  --- conversion, io

  procedure Assign(I: in Integer; S: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(Integer'Size,S.Val);
    MPFR_Set_Si(Rounded,S.Val,Long(I),GMP_RNDN);
    if Check_Flags and then Problem then Report("Assign"); end if;
  end Assign;

  procedure Assign(Q: in Rational; S: in out MPFloat) is
    D: constant LInt := LDen(Q);
    Rounded: Int;
    V: MPFR_T renames S.Val;
  begin
    Reformat(V);
    MPFR_Set_Si(Rounded,V,Long(LNum(Q)),GMP_RNDN);
    if D /= 1 then
      declare
        U: MPFR_T;
      begin
        MPFR_Init2(U,MP_Precision);
        MPFR_Set_Si(Rounded,U,Long(D),GMP_RNDN);
        MPFR_Div(Rounded,V,V,U,GMP_RNDN);
        MPFR_Clear(U);
      end;
    end if;
    if Check_Flags and then Num_Problem then Report("Scal"); end if;
  end Assign;

  procedure Assign(R: in Flt; S: in out MPFloat) is
    Rounded: Int;
  begin
    Reformat(LD_Precision,S.Val);
    MPFR_Set_Ld(Rounded,S.Val,Long_Double(R),GMP_RNDN);
    if Check_Flags and then Problem then Report("Assign"); end if;
  end Assign;

  function Scal(I: Integer) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,Integer'Size);
    MPFR_Set_Si(Rounded,V,Long(I),GMP_RNDN);
    if Check_Flags and then Problem then Report("Scal"); end if;
    return (Controlled with V);
  end Scal;

  function Scal(Q: Rational) return MPFloat is
    Rounded: Int;
    U,V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Set_Si(Rounded,V,Long(LNum(Q)),GMP_RNDN);
    MPFR_Init2(U,MP_Precision);
    MPFR_Set_Si(Rounded,U,Long(LDen(Q)),GMP_RNDN);
    MPFR_Div(Rounded,V,V,U,GMP_RNDN);
    MPFR_Clear(U);
    if Check_Flags and then Num_Problem then Report("Scal"); end if;
    return (Controlled with V);
  end Scal;

  function Scal(R: Flt) return MPFloat is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,LD_Precision);
    MPFR_Set_Ld(Rounded,V,Long_Double(R),GMP_RNDN);
    if Check_Flags and then Problem then Report("Scal"); end if;
    return (Controlled with V);
  end Scal;

  --- more conversion, io

  function NumDigits(P: MP_Prec_T; Base: Positive) return Size_T is
    N: constant Integer := 1073741823;      -- 2**30-1
    M: constant Integer := Width(N,Base)/30;
  begin
    return Size_T(M*Integer(P));
  end NumDigits;

  function EffLast(F: Positive; N: String) return Positive is
  begin
    for L in reverse F+1 .. N'Last loop
      if N(L) /= '0' then return L; end if;
    end loop;
    return F;
  end EffLast;

  function Image(S: MPFloat; Base: Positive; Digs: Natural := 0) return String is
    B: constant MPFR_Base := MPFR_Base(Base);
    EP: aliased constant MP_Exp_T_Ptr := new MP_Exp_T'(0);
    C: Chars_Ptr;
    D: Size_T := Size_T(Digs);
  begin
    if IsZero(S) then Just_Rounded := False; return " 0.0+0"; end if;
    if Digs=0 and then Strip0(Base,2) /= 1 then
      D := NumDigits(S.Val(1).MP_Prec,Base); -- full precision impossible
    end if;
    Clear_Flags;
    C := MPFR_Get_Str(Null_Ptr,EP,B,D,S.Val,GMP_RNDN);
    Just_Rounded := Raised_Inexact;
    declare
      E: constant Integer := Integer(EP.all-1);
      M: constant String := Value(C);
    begin
      MPFR_Free_Str(C);
      if MPFR_Sgn(S.Val)<0 then
        return '-' & M(2) & '.' & M(3 .. EffLast(3,M)) & Image(E,Base,True);
      else
        return '+' & M(1) & '.' & M(2 .. EffLast(2,M)) & Image(E,Base,True);
      end if;
    end;
  end Image;

  procedure Parse(N: in String; DigCount,ExpPos,NLast: out Integer) is
    --- assumes reasonable input
    DotPos: Integer := 0;
    C: Character;
  begin
    NLast := N'Last;
    for K in N'Range loop
      if N(K) = ',' then
        NLast := K-1;
        exit;
      end if;
    end loop;
    DigCount := 0;
    ExpPos := NLast+1;
    for K in N'First .. NLast loop
      C := N(K);
      if C = '.' then
        DotPos := K;
        exit;
      end if;
      if (C /= ' ') and then (C /= '+') and then (C /= '-') then
        DigCount := DigCount+1;
      end if;
    end loop;
    if DotPos=0 then raise Data_Error with Show0("MPFR.Floats.Parse: missing period"); end if;
    for K in DotPos+1 .. NLast loop
      C := N(K);
      if C = '+' or else C = '-' then
        ExpPos := K;
        exit;
      end if;
      DigCount := DigCount+1;
    end loop;
  end Parse;

  procedure Value(N: in String; Base: in Positive; R: in out MPFloat; Rounded: out Boolean) is
    Representable: Int; -- 0 means yes
    DigCount,ExpPos,NLast,ExpVal: Integer;
    P: MP_Prec_T;
    C: Chars_Ptr;
  begin
    Parse(N,DigCount,ExpPos,NLast);
    if ExpPos>NLast then
      C := New_String(N(N'First .. NLast) & "@+0");
    else
      ExpVal := Value(N(ExpPos .. NLast),Base);
      C := New_String(N(N'First .. ExpPos-1) & "@" & Image(ExpVal,10,True));
    end if;
    P := MP_Prec_T(Width(Base-1,2)*DigCount);
    P := Limb_Size*((P+Limb_Size-1)/Limb_Size);
    if not Check_Flags and then P>MP_Precision then P := MP_Precision; end if;
    Reformat(P,R.Val);
    Clear_Flags;
    MPFR_Set_Str(Representable,R.Val,C,MPFR_Base(Base),GMP_RNDN); -- exact if Base=2**K
    Rounded := Raised_Inexact or (Representable /= 0);
    Free(C);
  end Value;

  function IsInteger(S: MPFloat) return Boolean is
    I: Int; -- nonzero iff S is an integer
  begin
    MPFR_Integer_P(I,S.Val);
    return (I /= 0);
  end IsInteger;

  procedure Ceiling(S1: in MPFloat; S2: in out MPFloat) is
  begin
    MPFR_Ceil(S2.Val,S1.Val);
  end Ceiling;

  function Ceiling(S1: MPFloat) return MPFloat is
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Ceil(V,S1.Val);
    return (Controlled with V);
  end Ceiling;

  procedure Floor(S1: in MPFloat; S2: in out MPFloat) is
  begin
    MPFR_Floor(S2.Val,S1.Val);
  end Floor;

  function Floor(S1: MPFloat) return MPFloat is
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Floor(V,S1.Val);
    return (Controlled with V);
  end Floor;

  procedure Round(S1: in MPFloat; S2: in out MPFloat) is
  begin
    MPFR_Round(S2.Val,S1.Val);
  end Round;

  function Round(S1: MPFloat) return MPFloat is
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Round(V,S1.Val);
    return (Controlled with V);
  end Round;

  function Value(N: String; Base: Positive) return MPFloat is
    R: MPFloat;
  begin
    Value(N,Base,R,Just_Rounded);
    return R;
  end Value;

  procedure ValDec(N: in String; R: in out MPFloat) is
  begin
    Value(N,10,R,Just_Rounded);
  end ValDec;

  function ValDec(N: String) return MPFloat is
    R: MPFloat;
  begin
    Value(N,10,R,Just_Rounded);
    return R;
  end ValDec;

  procedure ValHex(N: in String; R: in out MPFloat; Rounded: out Boolean) is
  begin
    Value(N,16,R,Rounded);
    if Rounded and then IsZero(R) then
      Show0(N);
      raise Sorry with "MPFR.Floats.ValHex: should not happen";
    end if;
  end ValHex;

  function ValHex(N: String) return MPFloat is
    R: MPFloat;
  begin
    ValHex(N,R,Just_Rounded);
    return R;
  end ValHex;

  function DecStr(S: MPFloat; Digs: Natural := 0) return String is
  begin
    return Image(S,10,Digs);
  end DecStr;

  function HexStr(S: MPFloat; Digs: Natural := 0) return String is
  begin
    return Image(S,16,Digs);
  end HexStr;

  procedure Show1(N: in String; S: in MPFloat; NewLine: in Boolean := True) is
    Digs: constant Natural := 0; -- use default
  begin
    Show0(N & DecStr(S,Digs),NewLine);
  end Show1;

  procedure Show2(N: in String; S1,S2: in MPFloat; NewLine: in Boolean := True) is
    Digs: constant Natural := 0; -- use default
  begin
    Show0(N & DecStr(S1,Digs) & "   " & DecStr(S2,Digs),NewLine);
  end Show2;

  procedure Put(F: in File_Type; S: in MPFloat; Decimal: in Boolean := False) is
  begin
    if Decimal then
      Txt_IO.Put_Line(F,DecStr(S));
    else
      Txt_IO.Put_Line(F,HexStr(S));
    end if;
  end Put;

  procedure Get(F: in File_Type; S: in out MPFloat; Decimal: in Boolean := False) is
    V: constant String := Get_Next_Line(F);
  begin
    if Decimal then
      S := ValDec(V);
    else
      S := ValHex(V);
    end if;
  end Get;

  procedure Write(FileName: in String; S: in MPFloat; 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 MPFloat; 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;

  ---------- Rep

  procedure Val(R: in LLFloat; S: in out Rep) is
    Rounded: Int;
  begin
    MPFR_Set_Ld(Rounded,S.Val,Long_Double(R),GMP_RNDN);
    if Check_Flags and then Problem then Report("Val"); end if;
  end Val;

  function Val(R: LLFloat) return Rep is
    Rounded: Int;
    V: MPFR_T;
  begin
    MPFR_Init2(V,MP_Precision);
    MPFR_Set_Ld(Rounded,V,Long_Double(R),GMP_RNDN);
    if Check_Flags and then Problem then Report("Val"); end if;
    return (Controlled with V);
  end Val;

  function Upper(S: Rep) return LLFloat is
    U: constant Long_Double := MPFR_Get_Ld(S.Val,GMP_RNDU);
  begin
    if Check_Flags and then Problem then Report("Upper"); end if;
    return LLFloat(U);
  end Upper;

  function Lower(S: Rep) return LLFloat is
    L: constant Long_Double := MPFR_Get_Ld(S.Val,GMP_RNDD);
  begin
    if Check_Flags and then Problem then Report("Lower"); end if;
    return LLFloat(L);
  end Lower;

  procedure PutStd(S: in MPFloat; C,R: out LLFloat) is
    V: MPFR_T renames S.Val;
  begin
    C := LLFloat(MPFR_Get_Ld(V,GMP_RNDN));
    if Check_Flags then
      R := RHalf*(LLFloat(MPFR_Get_Ld(V,GMP_RNDU))-LLFloat(MPFR_Get_Ld(V,GMP_RNDD)));
      if Problem then Report("PutStd"); end if;
    end if;
  end PutStd;

  procedure PutStd(S: in MPFloat; C,R,B: out LLFloat) is
    V: MPFR_T renames S.Val;
  begin
    C := LLFloat(MPFR_Get_Ld(V,GMP_RNDN));
    if Check_Flags then
      R := RHalf*(LLFloat(MPFR_Get_Ld(V,GMP_RNDU))-LLFloat(MPFR_Get_Ld(V,GMP_RNDD)));
      if Problem then Report("PutStd"); end if;
    end if;
    B := RZero;
  end PutStd;

  procedure GetStd(C,Dummy: in LLFloat; S: in out MPFloat) is
    pragma Unreferenced(Dummy);
    Rounded: Int;
  begin
    MPFR_Set_Ld(Rounded,S.Val,Long_Double(C),GMP_RNDN);
    if Check_Flags and then Problem then Report("GetStd"); end if;
  end GetStd;

  procedure GetStd(C,Dummy1,Dummy2: in LLFloat; S: in out MPFloat) is
    pragma Unreferenced(Dummy1,Dummy2);
    Rounded: Int;
  begin
    MPFR_Set_Ld(Rounded,S.Val,Long_Double(C),GMP_RNDN);
    if Check_Flags and then Problem then Report("GetStd"); end if;
  end GetStd;

  --- misc

  function Exponent(S: in MPFloat) return Integer is
  begin
    return Integer(S.Val(1).MPFR_Exp);
  end Exponent;

  function Scaled(Eps: Flt) return Flt is
  begin
    return Scaled(Eps,Integer(MP_Precision));
  end Scaled;

  function Epsilon(Eps: Flt) return Flt is
  begin
    return Epsilon(Eps,Integer(MP_Precision));
  end Epsilon;

  procedure Set_Precision(Binary_Digits: in Positive) is
  begin
    MP_Precision := MP_Prec_T(Binary_Digits);
    if Check_Flags then
      if Binary_Digits<Long_Double'Machine_Mantissa then  -- avoid inexact conversions
        MP_Precision := LD_Precision;
        Show1("decreased precision only to Long_Double'Machine_Mantissa =",Integer(LD_Precision));
      elsif (Binary_Digits>Initial_Precision) and then (Verbosity>0) then
        Show1("Set_Precision warning: some packages may have been elaborated with Globals.Initial_Precision =",Initial_Precision);
      end if;
    end if;
  end Set_Precision;

  function Set_Precision(Binary_Digits: Positive) return Positive is
  begin
    Set_Precision(Binary_Digits);
    return Binary_Digits;
  end Set_Precision;

  procedure Set_Precision(Binary_Digits: in Positive; Dummy: in MPFloat) is
    pragma Unreferenced(Dummy);
  begin
    Set_Precision(Binary_Digits);
  end Set_Precision;

  function Get_Precision(R: MPFloat) return Positive is
  begin
    return Positive(R.Val(1).MP_Prec);
  end Get_Precision;

  procedure Set_Rounding_Mode(M: in Rounding_Mode) is
    --- ignores M.Direction
  begin
    Clear_Flags;
    Check_Flags := (M.Mask=Strict);
  end Set_Rounding_Mode;

  procedure Proper_Rounding is
  begin
    if Check_Flags then
      if Problem then Report("Proper_Rounding"); end if;
      MPFR_Clear_Flags;
      Check_Flags := False;
    end if;
    if A_Numeric_Mode /= null then A_Numeric_Mode.all; end if;
  end Proper_Rounding;

  procedure Proper_Rounding(Dummy: in MPFloat) is
    pragma Unreferenced(Dummy);
  begin
    if Check_Flags then
      if Problem then Report("Proper_Rounding"); end if;
      MPFR_Clear_Flags;
      Check_Flags := False;
    end if;
    if A_Numeric_Mode /= null then A_Numeric_Mode.all; end if;
  end Proper_Rounding;

  procedure Free_Cache(Dummy: in MPFloat) is
    pragma Unreferenced(Dummy);
  begin
    MPFR_Free_Cache;
  end Free_Cache;

  ---------------------------------------------------------------
  Primitive_Buffer: array (1 .. Primitive_Buffer_Size) of MPFloat;
  ---------------------------------------------------------------

  procedure Put_Numeric(N: in Positive; S: in Rep) is
  begin
    if N <= Primitive_Buffer_Size then
      Copy(S,Primitive_Buffer(N));
    else
      raise Sorry with "increase Primitive_Buffer_Size";
    end if;
  end Put_Numeric;

  procedure Get_Numeric(N: in Positive; S: in out Rep) is
  begin
    if N <= Primitive_Buffer_Size then
      Copy(Primitive_Buffer(N),S);
    else
      raise Sorry with "increase Primitive_Buffer_Size";
    end if;
  end Get_Numeric;

  ------------------------------------------------
  ------------------------------------------------
  package QP is new QPowers (Scalar => MPFloat);
  ------------------------------------------------
  ------------------------------------------------

  procedure QPower(Q: in Rational; S1: in MPFloat; S2: in out MPFloat) renames QP.QPower;

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

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

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

begin

  if Long_Double'Machine_Mantissa /= LLFloat'Machine_Mantissa then
    raise MPFR_Exception with "adapt precision of LLFloat";
  end if;

  if Long'Size < Integer'Size then
    raise MPFR_Exception with "size mismatch";
  end if;

end MPFR.Floats;
