with Standard_Mathematical_Functions;    use Standard_Mathematical_Functions;

package body Multprec_Floating_Numbers is

-- NOTES ON THE IMPLEMENTATION :
--   0) See the notes in the bodies of Multprec_{Natural,Integer}_Numbers.
--   1) 0.0 must have zero fraction and zero exponent.
--   2) The truncate and round operations all depend on the decimal radix.
--   3) The function "Decimal_Places" occupies a central place.
--   4) The normalized representation has a nonzero least significant digit.

  tol : constant double_float := 10.0**(-8);

-- AUXILIARIES :

  function Truncate ( f : double_float ) return integer is

    i : integer := integer(f);

  begin
    if i >= 0
     then if double_float(i) > f + tol
           then i := i-1;
          end if;
     else if double_float(i) < f - tol
           then i := i+1;
          end if;
    end if;
    return i;
  end Truncate;

  function Max_Size ( i1,i2 : Integer_Number ) return natural is

  -- DESCRIPTION :
  --   Return the maximal size of i1 and i2.

    sz1 : constant natural := Size(i1);
    sz2 : constant natural := Size(i2);

  begin
    if sz1 > sz2
     then return sz1;
     else return sz2;
    end if;
  end Max_Size;

  procedure Shift_Left ( n : in out Natural_Number; ns : out natural ) is

  -- DESCRIPTION :
  --   Shifts n until the leading coefficient is larger than Base/10.
  --   The number ns on return equals the number of shifts.

  begin
    if (Empty(n) or else Equal(n,0))
     then ns := 0;
     else ns := 0;
          while (Coefficient(n,Size(n)) < Base/10) loop
            Mul(n,10);
            ns := ns+1;
          end loop;
    end if;
  end Shift_Left;

  procedure Shift_Right ( n : in out Natural_Number; ns : out natural ) is

  -- DESCRIPTION :
  --   Shifts n until the last digit in the least significant coefficient
  --   of n is different from zero.
  --   The number ns on return equals the number of shifts.

    cff : natural;

  begin
    if (Empty(n) or else Equal(n,0))
     then ns := 0;
     else ns := 0;
          cff := Coefficient(n,0) mod 10;
          while (cff = 0) loop
            Div(n,10);
            ns := ns+1;
            cff := Coefficient(n,0) mod 10;
          end loop;
    end if;
  end Shift_Right;

  procedure Shift_Right ( i : in out Integer_Number; ns : out natural ) is

  -- DESCRIPTION :
  --   Applies the Shift_Right to the number representation of i.

    ni : Natural_Number;

  begin
    if Empty(i)
     then ns := 0;
     else --ni := Unsigned(i);
          Copy(Unsigned(i),ni);  -- avoid sharing
          Shift_Right(ni,ns);
          if Sign(i) < 0
           then Clear(i); i := Create(ni); Min(i);
           else Clear(i); i := Create(ni);
          end if;
          Clear(ni);
    end if;
  end Shift_Right;

  procedure Normalize ( f : in out Floating_Number ) is

  -- DESCRIPTION :
  --   Normalizes the representation of f so that its fraction has as
  --   rightmost digit a nonzero number.

    ns : natural;

  begin
    Shift_Right(f.fraction,ns);
    Add(f.exponent,ns);
  end Normalize;

  procedure Adjust_Expo ( f : in out double_float;
                          expo : in Integer_Number ) is

  -- DESCRIPTION :
  --   Multiplies f with 10**expo.

    cnt : Integer_Number;

  begin
    Copy(expo,cnt);
    while cnt > 0 loop
      f := f*10.0;
      Sub(cnt,1);
    end loop;
    while cnt < 0 loop
      f := f/10.0;
      Add(cnt,1);
    end loop;
    Clear(cnt);
  end Adjust_Expo;

-- CONSTRUCTORS :

  function Create ( i : integer ) return Floating_Number is

    res : Floating_Number;

  begin
    res.fraction := Create(i);
    res.exponent := Create(0);
    Normalize(res);
    return res;
  end Create;

  function Create ( i : Integer_Number ) return Floating_Number is

    res : Floating_Number;

  begin
    res.fraction := +i;
    res.exponent := Create(0);
    Normalize(res);
    return res;
  end Create;

  function Create ( fraction,exponent : integer ) return Floating_Number is

    res : Floating_Number;

  begin
    res.fraction := Create(fraction);
    res.exponent := Create(exponent);
    Normalize(res);
    return res;
  end Create;

  function Create ( fraction : Integer_Number;
                    exponent : integer ) return Floating_Number is

    res : Floating_Number;

  begin
    res.fraction := +fraction;
    res.exponent := Create(exponent);
    Normalize(res);
    return res;
  end Create;

  function Create ( fraction : integer;
                    exponent : Integer_Number ) return Floating_Number is

    res : Floating_Number;

  begin
    res.fraction := Create(fraction);
    res.exponent := +exponent;
    Normalize(res);
    return res;
  end Create;

  function Create ( fraction,exponent : Integer_Number )
                  return Floating_Number is

    res : Floating_Number;

  begin
    res.fraction := +fraction;
    res.exponent := +exponent;
    Normalize(res);
    return res;
  end Create;

  function Create ( f : double_float ) return Floating_Number is

  -- NOTE :
  --   To avoid rounding errors, comparison operations are prefered
  --   above relying on intermediate results of calculations, although
  --   rounding errors still occur!

    res : Floating_Number;
    ten : Integer_Number;
    wrkflt,fltexp,lowflt,uppflt,incflt : double_float;
    intexp : integer;
    fraction,exponent,intinc : Integer_Number;

  begin
    if f = 0.0
     then res := Create(0);
          return res;
    end if;
    if f > 0.0
     then wrkflt := f;
     else wrkflt := -f;
    end if;
    fltexp := LOG10(wrkflt);
    intexp := Truncate(fltexp);
    if intexp < 15
     then wrkflt := wrkflt*(10.0**(15-intexp));
          exponent := Create(intexp-15);
          intexp := 15;
     else exponent := Create(0);
    end if;
    lowflt := 10.0**intexp;
    uppflt := 10.0**(intexp+1);
    while lowflt > wrkflt loop
      intexp := intexp-1;
      lowflt := 10.0**intexp;
      uppflt := 10.0**(intexp+1);
    end loop;
    ten := Create(10);
    fraction := ten**intexp;
    for k in 1..16 loop
      incflt := 10.0**intexp;
      intinc := ten**intexp;
      uppflt := lowflt;
      for i in 1..10 loop
        uppflt := uppflt + incflt;
        exit when (uppflt > wrkflt);
        lowflt := uppflt;
        Add(fraction,intinc);
      end loop;
      intexp := intexp - 1;
      Clear(intinc);
    end loop;
    Clear(ten);
    if f < 0.0
     then Min(fraction);
    end if;
    res := Create(fraction,exponent);
    Normalize(res);
    return res;
  end Create;

-- SELECTORS :

  function Fraction ( f : Floating_Number ) return Integer_Number is
  begin
    return f.fraction;
  end Fraction;

  function Exponent ( f : Floating_Number ) return Integer_Number is
  begin
    return f.exponent;
  end Exponent;

  function Size_Fraction ( f : Floating_Number ) return natural is
  begin
    return Size(Fraction(f));
  end Size_Fraction;

  function Size_Exponent ( f : Floating_Number ) return natural is
  begin
    return Size(Exponent(f));
  end Size_Exponent;

  function Decimal_Places_Fraction ( f : Floating_Number ) return natural is
  begin
    return Decimal_Places(Fraction(f));
  end Decimal_Places_Fraction;

  function Decimal_Places_Exponent ( f : Floating_Number ) return natural is
  begin
    return Decimal_Places(Exponent(f));
  end Decimal_Places_Exponent;

  function Decimal_to_Size ( deci : natural ) return natural is

    res : natural;
    exp : constant natural := Multprec_Natural_Numbers.Exponent;

  begin
    res := deci/exp;
    if res*exp < deci
     then res := res+1;
    end if;
    return res;
  end Decimal_to_Size;

  function AbsVal ( f : Floating_Number ) return Floating_Number is

    res : Floating_Number;

  begin
    if Negative(f.fraction)
     then res.fraction := -f.fraction;
     else res.fraction := +f.fraction;
    end if;
    res.exponent := +f.exponent;
    return res;
  end AbsVal;

-- COMPARISON AND COPYING :

  function Equal ( f1 : Floating_Number; f2 : double_float ) return boolean is

    ff2 : Floating_Number := Create(f2);
    res : boolean := Equal(f1,ff2);

  begin
    Clear(ff2);
    return res;
  end Equal;

  function Equal ( f1,f2 : Floating_Number ) return boolean is

    res : boolean;
    f1deci : constant natural := Decimal_Places(f1.fraction);
    f2deci : constant natural := Decimal_Places(f2.fraction);
    f1expo : Integer_Number := f1.exponent + f1deci;
    f2expo : Integer_Number := f2.exponent + f2deci;
    mulfra : Integer_Number;

  begin
    if Equal(f1expo,f2expo)
     then if f1deci = f2deci
           then res := Equal(f1.fraction,f2.fraction);
           else if f1deci < f2deci
                 then mulfra := f1.fraction*10;
                      for i in 1..(f2deci-f1deci-1) loop
                        Mul(mulfra,10);
                      end loop;
                      res := Equal(mulfra,f2.fraction);
                 else mulfra := f2.fraction*10;
                      for i in 1..(f1deci-f2deci-1) loop
                        Mul(mulfra,10);
                      end loop;
                      res := Equal(f1.fraction,mulfra);
                end if;
                Clear(mulfra);
          end if;
     else res := false;
    end if;
    Clear(f1expo); Clear(f2expo);
    return res;
  end Equal;

  function "<" ( f1 : double_float; f2 : Floating_Number ) return boolean is

    ff1 : Floating_Number := Create(f1);
    res : boolean := (ff1 < f2);

  begin
    Clear(ff1);
    return res;
  end "<";

  function "<" ( f1 : Floating_Number; f2 : double_float ) return boolean is

    ff2 : Floating_Number := Create(f2);
    res : boolean := (f1 < ff2);

  begin
    Clear(ff2);
    return res;
  end "<";

  function "<" ( f1,f2 : Floating_Number ) return boolean is

    res : boolean;
    f1deci : constant natural := Decimal_Places(f1.fraction);
    f2deci : constant natural := Decimal_Places(f2.fraction);
    f1expo : Integer_Number := f1.exponent + f1deci;
    f2expo : Integer_Number := f2.exponent + f2deci;
    mulfra : Integer_Number;

  begin
    if (f1expo < f2expo)
     then if Sign(f2.fraction) > 0
           then res := true;
           elsif Sign(f2.fraction) < 0
               then res := false;
               else res := (Sign(f1.fraction) < 0);
          end if;
     elsif (f1expo > f2expo)
         then if Sign(f1.fraction) > 0
               then res := false;
               elsif Sign(f1.fraction) < 0
                   then res := true;
                   else res := (Sign(f2.fraction) > 0);
              end if;
         else if f1deci = f2deci
               then res := (f1.fraction < f2.fraction);
               else if f1deci < f2deci
                     then mulfra := f1.fraction*10;
                          for i in 1..(f2deci-f1deci-1) loop
                            Mul(mulfra,10);
                          end loop;
                          res := (mulfra < f2.fraction);
                     else mulfra := f2.fraction*10;
                          for i in 1..(f1deci-f2deci-1) loop
                            Mul(mulfra,10);
                          end loop;
                          res := (f1.fraction < mulfra);
                    end if;
                    Clear(mulfra);
              end if;
    end if;
    Clear(f1expo); Clear(f2expo);
    return res;
  end "<";

  function ">" ( f1 : double_float; f2 : Floating_Number ) return boolean is

    ff1 : Floating_Number := Create(f1);
    res : boolean := (ff1 > f2);

  begin
    Clear(ff1);
    return res;
  end ">";

  function ">" ( f1 : Floating_Number; f2 : double_float ) return boolean is

    ff2 : Floating_Number := Create(f2);
    res : boolean := (f1 > ff2);

  begin
    Clear(ff2);
    return res;
  end ">";

  function ">" ( f1,f2 : Floating_Number ) return boolean is

    res : boolean;
    f1deci : constant natural := Decimal_Places(f1.fraction);
    f2deci : constant natural := Decimal_Places(f2.fraction);
    f1expo : Integer_Number := f1.exponent + f1deci;
    f2expo : Integer_Number := f2.exponent + f2deci;
    mulfra : Integer_Number;

  begin
    if (f1expo > f2expo)
     then if Sign(f1.fraction) > 0
           then res := true;
           elsif Sign(f1.fraction) < 0
               then res := false;
               else res := (Sign(f2.fraction) < 0);
          end if;
     elsif (f1expo < f2expo)
         then if Sign(f2.fraction) > 0
               then res := false;
               elsif Sign(f2.fraction) < 0
                   then res := true;
                   else res := (Sign(f1.fraction) > 0);
              end if;
         else if f1deci = f2deci
               then res := (f1.fraction > f2.fraction);
               else if f1deci < f2deci
                     then mulfra := f1.fraction*10;
                          for i in 1..(f2deci-f1deci-1) loop
                            Mul(mulfra,10);
                          end loop;
                          res := (mulfra > f2.fraction);
                     else mulfra := f2.fraction*10;
                          for i in 1..(f1deci-f2deci-1) loop
                            Mul(mulfra,10);
                          end loop;
                          res := (f1.fraction > mulfra);
                    end if;
                    Clear(mulfra);
              end if;
    end if;
    Clear(f1expo); Clear(f2expo);
    return res;
  end ">";

  procedure Copy ( f1 : in Floating_Number; f2 : in out Floating_Number ) is
  begin
    Clear(f2);
    Copy(f1.fraction,f2.fraction);
    Copy(f1.exponent,f2.exponent);
  end Copy;

-- EXPANDING AND SHORTENING THE MANTISSA :

  function Expand ( i : Integer_Number; k : natural )
                  return Array_of_Naturals is

  -- DESCRIPTION :
  --   Returns the expanded coefficient vector of i.

    cff : constant Array_of_Naturals := Coefficients(Unsigned(i));
    res : Array_of_Naturals(cff'first..cff'last+k);

  begin
    res(cff'range) := cff(cff'range);
    res(cff'last+1..res'last) := (cff'last+1..res'last => 0);
    return res;
  end Expand;

  procedure Truncate ( i : in Integer_Number; k : in natural;
                       trc : out Array_of_Naturals; exp : out integer ) is

  -- DESCRIPTION :
  --   On return is the truncated coefficient vector of i and the number
  --   of decimal places that are lost, or gained when exp is positive.

    n : Natural_Number;
    ns : natural;          -- #shifts = #decimal places gained

  begin
    Copy(Unsigned(i),n);
    Shift_Left(n,ns);
    exp := k*Multprec_Natural_Numbers.Exponent - ns;
    declare
      cff : constant Array_of_Naturals := Coefficients(n);
    begin
      for i in trc'range loop
        trc(i) := cff(i+k);
      end loop;
      Clear(n);
    end;
  end Truncate;

  procedure Truncate ( i : in out Integer_Number; k : in natural;
                       exp : out integer ) is

    n : Natural_Number := Unsigned(i);
    neg : constant boolean := Negative(i);
    ns : natural;
    newint : Integer_Number;

  begin
    Shift_Left(n,ns);
    exp := k*Multprec_Natural_Numbers.Exponent - ns;
    declare
      cff : constant Array_of_Naturals := Coefficients(n);
      trc : Array_of_Naturals(cff'first..cff'last-k);
    begin
      for i in trc'range loop
        trc(i) := cff(i+k);
      end loop;
      newint := Create(trc);
      i := newint; --Copy(newint,i); --Clear(newint);
      if neg
       then Min(i);
      end if;
    end;
    Clear(n);
  end Truncate;

  procedure Round ( i : in Integer_Number; k : in natural;
                    trc : out Array_of_Naturals; exp : out integer ) is

  -- DESCRIPTION :
  --   On return is the rounded coefficient vector of i and the number
  --   of decimal places that are lost, or gained when exp is positive.

    n : Natural_Number := Unsigned(i);
    ns : natural;          -- #shifts = #decimal places gained

  begin
    Shift_Left(n,ns);
    exp := k*Multprec_Natural_Numbers.Exponent - ns;
    declare
      cff : constant Array_of_Naturals := Coefficients(n);
    begin
      for i in trc'range loop
        trc(i) := cff(i+k);
      end loop;
      if k > 0
       then if cff(k-1) >= Multprec_Natural_Numbers.Base/2
             then trc(0) := trc(0)+1;
            end if;
      end if;
    end;
    Clear(n);
  end Round;

  function Expand ( f : Floating_Number; k : natural ) return Floating_Number is

    res : Floating_Number;
    cff : constant Array_of_Naturals := Expand(f.fraction,k);

  begin
    res.fraction := Create(Create(cff));
    if Negative(f.fraction)
     then Min(res.fraction);
    end if;
    Copy(f.exponent,res.exponent);
    Normalize(res);
    return res;
  end Expand;

  procedure Expand ( f : in out Floating_Number; k : in natural ) is

    cff : constant Array_of_Naturals := Expand(f.fraction,k);
    neg : constant boolean := Negative(f.fraction);

  begin
    Clear(f.fraction);
    f.fraction := Create(Create(cff));
    if neg
     then Min(f.fraction);
    end if;
    Normalize(f);
  end Expand;

  function Round ( f : Floating_Number ) return double_float is

    res : double_float;
    szf,ns : natural;
    n : Natural_Number;
    expo : Integer_Number;

  begin
    if (Empty(f.fraction) or else Equal(f.fraction,0))
     then res := 0.0;
     else Copy(Unsigned(f.fraction),n);
          Shift_Left(n,ns);
          expo := f.exponent - ns;
          szf := Size(n);
          while Coefficient(n,szf) = 0 and szf > 0 loop
            szf := szf - 1;
          end loop;
          res := double_float(Coefficient(n,szf));
          if szf >= 1
           then res := res*double_float(Multprec_Natural_Numbers.Base)
                     + double_float(Coefficient(n,szf-1));
                if szf >= 2
                 then if Coefficient(n,szf-2) > Multprec_Natural_Numbers.Base/2
                       then res := res+1.0;
                      end if;
                      ns := (szf - 1)*Multprec_Natural_Numbers.Exponent;
                      Add(expo,ns);
                end if;
          end if;
          Adjust_Expo(res,expo);
          Clear(n); Clear(expo);
    end if;
    return res;
  end Round;

  function Round ( f : Floating_Number; k : natural ) return Floating_Number is

    res : Floating_Number;
    cff : Array_of_Naturals(0..(Size(f.fraction)-k));
    exp : integer;

  begin
    Round(f.fraction,k,cff,exp);
    res.fraction := Create(Create(cff));
    if Negative(f.fraction)
     then Min(res.fraction);
    end if;
    Copy(f.exponent,res.exponent);
    Add(res.exponent,exp);
    Normalize(res);
    return res;
  end Round;

  procedure Round ( f : in out Floating_Number; k : in natural ) is

    exp : integer;
    n : Natural_Number;
    neg : constant boolean := Negative(f.fraction);
    ns : natural;

  begin
    Copy(Unsigned(f.fraction),n);
    Shift_Left(n,ns);
    exp := k*Multprec_Natural_Numbers.Exponent - ns;
    declare
      cff : constant Array_of_Naturals := Coefficients(n);
      trc : Array_of_Naturals(cff'first..cff'last-k);
    begin
      for i in trc'range loop
        trc(i) := cff(i+k);
      end loop;
      if k > 0
       then if cff(k-1) >= Multprec_Natural_Numbers.Base/2
             then trc(0) := trc(0)+1;
            end if;
      end if;
      Clear(f.fraction); 
      f.fraction := Create(trc);
      if neg
       then Min(f.fraction);
      end if;
    end;
    Clear(n);
    Add(f.exponent,exp);
    Normalize(f);
  end Round;

  function Trunc ( f : Floating_Number ) return double_float is

    res : double_float;
    szf,ns : natural;
    n : Natural_Number;
    expo : Integer_Number;

  begin
    if (Empty(f.fraction) or else Equal(f.fraction,0))
     then res := 0.0;
     else n := Unsigned(f.fraction);
          Shift_Left(n,ns);
          expo := f.exponent - ns;
          szf := Size(n);
          while Coefficient(n,szf) = 0 and szf > 0 loop
            szf := szf - 1;
          end loop;
          res := double_float(Coefficient(f.fraction,szf));
          if szf >= 1
           then res := res*double_float(Multprec_Natural_Numbers.Base)
                     + double_float(Coefficient(f.fraction,szf-1));
                if szf >= 2
                 then ns := (szf - 1)*Multprec_Natural_Numbers.Exponent;
                      Add(expo,ns);
                end if;
          end if;
          Adjust_Expo(res,expo);
          Clear(n); Clear(expo);
    end if;
    return res;
  end Trunc;

  function Trunc ( f : Floating_Number; k : natural ) return Floating_Number is
 
    res : Floating_Number;
    cff : Array_of_Naturals(0..(Size(f.fraction)-k));
    exp : integer;

  begin
    Truncate(f.fraction,k,cff,exp);
    res.fraction := Create(Create(cff));
    if Negative(f.fraction)
     then Min(res.fraction);
    end if;
    Copy(f.exponent,res.exponent);
    Add(res.exponent,exp);
    return res;
  end Trunc;

  procedure Trunc ( f : in out Floating_Number; k : in natural ) is

    exp : integer;
    n : Natural_Number;
    neg : constant boolean := Negative(f.fraction);
    ns : natural;

  begin
    Copy(Unsigned(f.fraction),n);
    Shift_Left(n,ns);
    exp := k*Multprec_Natural_Numbers.Exponent - ns;
    declare
      cff : constant Array_of_Naturals := Coefficients(n);
      trc : Array_of_Naturals(cff'first..cff'last-k);
    begin
      for i in trc'range loop
        trc(i) := cff(i+k);
      end loop;
      Clear(f.fraction);
      f.fraction := Create(trc);
      if neg
       then Min(f.fraction);
      end if;
    end;
    Clear(n);
    Add(f.exponent,exp);
  end Trunc;

  procedure Set_Size ( f : in out Floating_Number; k : in natural ) is

    sz : constant natural := Size_Fraction(f);

  begin
    if sz > k
     then Round(f,sz-k);
     elsif sz < k
         then Expand(f,k-sz);
    end if;
  end Set_Size;

-- ARITHMETIC OPERATIONS as functions (no data sharing) :

  function "+" ( f1 : Floating_Number; f2 : double_float )
               return Floating_Number is

    ff2 : Floating_Number := Create(f2);
    res : Floating_Number := f1+ff2;

  begin
    Clear(ff2);
    return res;
  end "+";

  function "+" ( f1 : double_float; f2 : Floating_Number )
               return Floating_Number is
  begin
    return f2+f1;
  end "+";

  function "+" ( f1,f2 : Floating_Number ) return Floating_Number is

    res : Floating_Number;
    maxsz,diff_size,max_diff : integer;
    diff_expo : Integer_Number;
    cnt : natural;

  begin
    if (Empty(f1.fraction) or else Equal(f1.fraction,0))
     then Copy(f2,res);
     elsif (Empty(f2.fraction) or else Equal(f2.fraction,0))
         then Copy(f1,res);
         else maxsz := Max_Size(f1.fraction,f2.fraction) + 1;
              if Equal(f1.exponent,f2.exponent)
               then res.fraction := f1.fraction + f2.fraction;
                    Copy(f1.exponent,res.exponent);
               else diff_expo := f1.exponent - f2.exponent;
                    max_diff := maxsz*Multprec_Natural_Numbers.Exponent;
                    if diff_expo > 0
                     then Copy(f1.fraction,res.fraction);
                          if diff_expo < 2*max_diff
                           then cnt := Create(diff_expo);
                                for i in 1..cnt loop
                                  Mul(res.fraction,10);
                                end loop;
                                Add(res.fraction,f2.fraction);
                                Copy(f2.exponent,res.exponent);
                           else Copy(f1.exponent,res.exponent);
                          end if;
                     else Copy(f2.fraction,res.fraction);
                          Min(diff_expo);
                          if diff_expo < 2*max_diff
                           then cnt := Create(diff_expo);
                                for i in 1..cnt loop
                                  Mul(res.fraction,10);
                                end loop;
                                Add(res.fraction,f1.fraction);
                                Copy(f1.exponent,res.exponent);
                           else Copy(f2.exponent,res.exponent);
                          end if;
                    end if;
                    Clear(diff_expo);
              end if;
              diff_size := Size(res.fraction) + 1 - maxsz;
              if diff_size > 0
               then Round(res,diff_size);
               else Normalize(res);
              end if;
    end if;
    return res;
  end "+";

  function "+" ( f : Floating_Number ) return Floating_Number is

    res : Floating_Number;

  begin
    Copy(f,res);
    return res;
  end "+";

  function "-" ( f : Floating_Number ) return Floating_Number is

    res : Floating_Number;

  begin
    Copy(f,res);
    Min(res.fraction);
    return res;
  end "-";

  function "-" ( f1 : Floating_Number; f2 : double_float )
               return Floating_Number is

    minf2 : double_float := -f2;
    res : Floating_Number := f1+minf2;

  begin
    return res;
  end "-";

  function "-" ( f1 : double_float; f2 : Floating_Number )
               return Floating_Number is

    ff1 : Floating_Number := Create(f1);
    res : Floating_Number := ff1 - f2;

  begin
    Clear(ff1);
    return res;
  end "-";

  function "-" ( f1,f2 : Floating_Number ) return Floating_Number is

    res,minf2 : Floating_Number;

  begin
    minf2.exponent := f2.exponent;
    minf2.fraction := -f2.fraction;
    res := f1 + minf2;
    Clear(minf2.fraction);
    return res;
  end "-";

  function "*" ( f1 : Floating_Number; f2 : double_float )
               return Floating_Number is

    ff2 : Floating_Number := Create(f2);
    res : Floating_Number := f1*ff2;

  begin
    Clear(ff2);
    return res;
  end "*";

  function "*" ( f1 : double_float; f2 : Floating_Number )
               return Floating_Number is

    ff1 : Floating_Number := Create(f1);
    res : Floating_Number := ff1*f2;

  begin
    Clear(ff1);
    return res;
  end "*";

  function "*" ( f1,f2 : Floating_Number ) return Floating_Number is

    res : Floating_Number;
    diffsize : integer;

  begin
    if not (Empty(f1.fraction) or else Equal(f1.fraction,0))
     then if not (Empty(f2.fraction) or else Equal(f2.fraction,0))
           then res.fraction := f1.fraction*f2.fraction;
                res.exponent := f1.exponent+f2.exponent;
                diffsize := Size(res.fraction)
                          - Max_Size(f1.fraction,f2.fraction);
                if diffsize > 0
                 then Round(res,diffsize);
                end if;
                Normalize(res);
          end if;
    end if;
    return res;
  end "*";

  function "**" ( f : double_float; n : Natural_Number )
                return Floating_Number is

    ff : Floating_Number := Create(f);
    res : Floating_Number := ff**n;

  begin
    Clear(ff);
    return res;
  end "**";

  function "**" ( f : Floating_Number; n : Natural_Number )
                return Floating_Number is

    res : Floating_Number;
    cnt : Natural_Number;

  begin
    if (Empty(n) or else Equal(n,0))
     then res := Create(1);
     else Copy(f,res);
          cnt := Create(1);
          while cnt < n loop
            Mul(res,f);
            Add(cnt,1);
          end loop;
          Clear(cnt);
          Normalize(res);
    end if;
    return res;
  end "**";

  function "**" ( f : Floating_Number; i : integer ) return Floating_Number is

    res : Floating_Number;

  begin
    if i = 0
     then res := Create(1.0);
     else if i > 0
           then Copy(f,res);
                for j in 1..(i-1) loop
                  Mul(res,f);
                end loop;
           else res := Create(1);
                for j in 1..(-i) loop
                  Div(res,f);
                end loop;
          end if;
          Normalize(res);
    end if;
    return res;
  end "**";

  function "**" ( f : double_float; i : Integer_Number )
                return Floating_Number is

    ff : Floating_Number := Create(f);
    res : Floating_Number := ff**i;

  begin
    Clear(ff);
    return res;
  end "**";

  function "**" ( f : Floating_Number; i : Integer_Number )
                return Floating_Number is

    res : Floating_Number;
    cnt,n : Natural_Number;

  begin
    if (Empty(i) or else Equal(i,0))
     then res := Create(1);
     else n := Unsigned(i);
          if i > 0
           then Copy(f,res);
                cnt := Create(1);
                while cnt < n loop
                  Mul(res,f);
                  Add(cnt,1);
                end loop;
           else res := Create(1);
                cnt := Create(0);
                while cnt < n loop
                  Div(res,f);
                  Add(cnt,1);
                end loop;
          end if;
          Clear(cnt);
          Normalize(res);
    end if;
    return res;
  end "**";

  function "/" ( f1 : Floating_Number; f2 : double_float )
               return Floating_Number is

    ff2 : Floating_Number := Create(f2);
    res : Floating_Number := f1/ff2;

  begin
    Clear(ff2);
    return res;
  end "/";

  function "/" ( f1 : double_float; f2 : Floating_Number )
               return Floating_Number is

    ff1 : Floating_Number := Create(f1);
    res : Floating_Number := ff1/f2;

  begin
    Clear(ff1);
    return res;
  end "/";

  function Pos_Div ( f1,f2 : Floating_Number ) return Floating_Number is

  -- DESCRIPTION :
  --   This division only applies when f1 and f2 are both positive.

    res : Floating_Number;
    rmd_frac,new_rmd_frac : Integer_Number;
    diffsize : integer;
    maxsz : natural := Max_Size(f1.fraction,f2.fraction);

  begin
    Div(f1.fraction,f2.fraction,res.fraction,rmd_frac);
    res.exponent := f1.exponent - f2.exponent;
    diffsize := Size(res.fraction) - maxsz;
    if not Equal(rmd_frac,0)
     then while diffsize <= 0 loop
            while rmd_frac < f2.fraction loop
              Mul(rmd_frac,10);
              Mul(res.fraction,10);
              Sub(res.exponent,1);
            end loop;
            Div(rmd_frac,f2.fraction,new_rmd_frac);
            Add(res.fraction,rmd_frac);
            Clear(rmd_frac); rmd_frac := new_rmd_frac;
            diffsize := Size(res.fraction) - maxsz;
            exit when Equal(rmd_frac,0);
          end loop;
    end if;
    Clear(rmd_frac);
    if diffsize > 0
     then Round(res,diffsize);
    end if;
    Normalize(res);
    return res;
  end Pos_Div;

  function "/" ( f1,f2 : Floating_Number ) return Floating_Number is

    res,minf1,minf2 : Floating_Number;

  begin
    if not (Empty(f1.fraction) or else Equal(f1.fraction,0))
     then if not (Empty(f2.fraction) or else Equal(f2.fraction,0))
           then if Multprec_Integer_Numbers.Positive(f1.fraction)
                 then if Multprec_Integer_Numbers.Positive(f2.fraction)
                       then res := Pos_Div(f1,f2);
                       else minf2.fraction := -f2.fraction;
                            minf2.exponent := f2.exponent;
                            res := Pos_Div(f1,minf2);
                            Clear(minf2.fraction);
                            Min(res);
                      end if;
                 else minf1.fraction := -f1.fraction;
                      minf1.exponent := f1.exponent;
                      if Multprec_Integer_Numbers.Positive(f2.fraction)
                       then res := Pos_Div(minf1,f2);
                            Min(res);
                       else minf2.fraction := -f2.fraction;
                            minf2.exponent := f2.exponent;
                            res := Pos_Div(minf1,minf2);
                            Clear(minf2.fraction);
                      end if;
                      Clear(minf1.fraction);
                end if;
           else raise NUMERIC_ERROR;
          end if;
    end if;
    return res;
  end "/";

-- ARITHMETIC OPERATIONS as procedures for memory management :

  procedure Add ( f1 : in out Floating_Number; f2 : in double_float ) is

    ff2 : Floating_Number := Create(f2);

  begin
    Add(f1,ff2);
    Clear(ff2);
  end Add;

  procedure Add ( f1 : in out Floating_Number; f2 : in Floating_Number ) is

    res : Floating_Number;
    maxsz,diff_size,max_diff : integer;
    diff_expo : Integer_Number;
    cnt : natural;

  begin
    if (Empty(f1.fraction) or else Equal(f1.fraction,0))
     then Copy(f2,f1);
     elsif (Empty(f2.fraction) or else Equal(f2.fraction,0))
         then null;
         else maxsz := Max_Size(f1.fraction,f2.fraction) + 1;
              if Equal(f1.exponent,f2.exponent)
               then Add(f1.fraction,f2.fraction);
               else diff_expo := f1.exponent - f2.exponent;
                    max_diff := maxsz*Multprec_Natural_Numbers.Exponent;
                    if diff_expo > 0
                     then if diff_expo < 2*max_diff
                           then cnt := Create(diff_expo);
                                for i in 1..cnt loop
                                  Mul(f1.fraction,10);
                                end loop;
                                Add(f1.fraction,f2.fraction);
                                Copy(f2.exponent,f1.exponent);
                           else null;
                          end if;
                     else Copy(f2.fraction,res.fraction);
                          Min(diff_expo);
                          if diff_expo < 2*max_diff
                           then cnt := Create(diff_expo);
                                for i in 1..cnt loop
                                  Mul(res.fraction,10);
                                end loop;
                                Add(f1.fraction,res.fraction);
                                Clear(res.fraction);
                           else Copy(f2.exponent,f1.exponent);
                                Clear(f1.fraction);
                                f1.fraction := res.fraction;
                          end if;
                    end if;
                    Clear(diff_expo);
              end if;
              diff_size := Size(f1.fraction) + 1 - maxsz;
              if diff_size > 0
               then Round(f1,diff_size);
               else Normalize(f1);
              end if;
    end if;
  end Add;

  procedure Sub ( f1 : in out Floating_Number; f2 : in double_float ) is

    minf2 : double_float := -f2;

  begin
    Add(f1,minf2);
  end Sub;

  procedure Sub ( f1 : in out Floating_Number; f2 : in Floating_Number ) is

    minf2 : Floating_Number;

  begin
    minf2.exponent := f2.exponent;
    minf2.fraction := -f2.fraction;
    Add(f1,minf2);
    Clear(minf2.fraction);
  end Sub;

  procedure Min ( f : in out Floating_Number ) is
  begin
    Min(f.fraction);
  end Min;

  procedure Mul ( f1 : in out Floating_Number; f2 : in double_float ) is

    ff2 : Floating_Number := Create(f2);

  begin
    Mul(f1,ff2);
    Clear(ff2);
  end Mul;

  procedure Mul ( f1 : in out Floating_Number; f2 : in Floating_Number ) is

    diffsize,maxsize : integer;

  begin
    if Empty(f2.fraction)
     then Clear(f1);
     elsif Equal(f2.fraction,0)
         then Clear(f1);
         elsif Empty(f1.fraction)
             then null;
             elsif Equal(f1.fraction,0)
                 then null;
                 else maxsize := Max_Size(f1.fraction,f2.fraction);
                      Mul(f1.fraction,f2.fraction);
                      Add(f1.exponent,f2.exponent);
                      diffsize := Size(f1.fraction) - maxsize;
                      if diffsize > 0
                       then Round(f1,diffsize);
                      end if;
                      Normalize(f1);
    end if;
  end Mul;

  procedure Div ( f1 : in out Floating_Number; f2 : in double_float ) is

    ff2 : Floating_Number := Create(f2);

  begin
    Div(f1,ff2);
    Clear(ff2);
  end Div;

  procedure Pos_Div ( f1 : in out Floating_Number; f2 : in Floating_Number ) is

  -- DESCRIPTION :
  --   Does the division in case f1 and f2 are both positive.

    rmd_frac,new_rmd_frac : Integer_Number;
    diffsize : integer;
    maxsz : natural := Max_Size(f1.fraction,f2.fraction);

  begin
    Div(f1.fraction,f2.fraction,rmd_frac);
    Sub(f1.exponent,f2.exponent);
    diffsize := Size(f1.fraction) - maxsz;
    if not Equal(rmd_frac,0)
     then while diffsize <= 0 loop
            while rmd_frac < f2.fraction loop
              Mul(rmd_frac,10); 
              Mul(f1.fraction,10);
              Sub(f1.exponent,1);
            end loop;
            Div(rmd_frac,f2.fraction,new_rmd_frac);
            Add(f1.fraction,rmd_frac);
            Clear(rmd_frac); rmd_frac := new_rmd_frac;
            diffsize := Size(f1.fraction) - maxsz;
            exit when Equal(rmd_frac,0);
          end loop;
    end if;
    Clear(rmd_frac);
    if diffsize > 0
     then Round(f1,diffsize);
    end if;
    Normalize(f1);
  end Pos_Div;

  procedure Div ( f1 : in out Floating_Number; f2 : in Floating_Number ) is

    minf1,minf2 : Floating_Number;

  begin
    if not (Empty(f1.fraction) or else Equal(f1.fraction,0))
     then if not (Empty(f2.fraction) or else Equal(f2.fraction,0))
           then if Multprec_Integer_Numbers.Positive(f1.fraction)
                 then if Multprec_Integer_Numbers.Positive(f2.fraction)
                       then Pos_Div(f1,f2);
                       else minf2.fraction := -f2.fraction;
                            minf2.exponent := f2.exponent;
                            Pos_Div(f1,minf2);
                            Clear(minf2.fraction); Min(f1);
                      end if;
                 else minf1 := -f1;
                      if Multprec_Integer_Numbers.Positive(f2.fraction)
                       then Pos_Div(minf1,f2);
                            Clear(f1); f1 := minf1; Min(f1);
                       else minf2.fraction := -f2.fraction;
                            minf2.exponent := f2.exponent;
                            Pos_Div(minf1,minf2);
                            Clear(minf2.fraction);
                            Clear(f1); f1 := minf1;
                      end if;
                end if;
           else raise NUMERIC_ERROR;
          end if;
    end if;
  end Div;

-- DESTRUCTOR :

  procedure Clear ( f : in out Floating_Number ) is
  begin
    Clear(f.fraction);
    Clear(f.exponent);
  end Clear;

end Multprec_Floating_Numbers;