with unchecked_deallocation;

package body Multprec_Natural_Numbers is

-- NOTES ON THE CHOICE OF REPRESENTATION AND IMPLEMENTATION :
--   1) The base is decimal, for convenience of input/output.
--      Everything should remain correct if another base is chosen,
--      except of course the function "Decimal_Places".
--      The base is such that factor times the base is still a standard
--      natural number, so that standard arithmetic can be used.
--   2) A natural number is implemented by means of a coefficient array.
--      The advantage over a linked-list representation is that the
--      coefficients can be accessed faster.  The disadvantage of some
--      wasted storage is countered by the aim of implementing higher
--      precision floating-point numbers, where the mantissa is fixed.
--   3) The usage of a pointer to access the coefficient vector permits
--      to have uniformity in input/output, since on input it is not
--      known in advance how long a natural number may be.
--   4) A natural number n for which Empty(n) holds is considered as zero.

  fact : constant natural := 10;
  expo : constant natural := 8;
  the_base : constant natural := fact**expo;

-- DATATRUCTURES :

  type Natural_Number_Rep ( size : natural ) is record
    numb : Array_of_Naturals(0..size);
  end record;

-- AUXILIARY :

  function Size ( n : natural ) return natural is

  -- DESCRIPTION :
  --   Determines the size of the natural number to represent n.

    acc : natural := the_base;

  begin
    for i in 0..n loop
      if n < acc
       then return i;
       else acc := acc*the_base;
      end if;
    end loop;
    return n;
  end Size;

-- CREATORS :

  function Create ( n : natural ) return Array_of_Naturals is

    res : Array_of_Naturals(0..Size(n));
    remainder : natural := n;

  begin
    for i in res'range loop
      res(i) := remainder mod the_base;
      remainder := remainder/the_base;
    end loop;
    return res;
  end Create;

  function Create ( n : natural ) return Natural_Number is

    res : Natural_Number := Create(Create(n));

  begin
    return res;
  end Create;

  function Create ( n : Array_of_Naturals ) return Natural_Number is

    res : Natural_Number;
    res_rep : Natural_Number_Rep(n'last);

  begin
    res_rep.numb := n;
    res := new Natural_Number_Rep'(res_rep);
    return res;
  end Create;

  function Create ( n : Natural_Number ) return natural is

    res : natural;

  begin
    if Empty(n)
     then res := 0;
     else res := Coefficient(n,Size(n));
          for i in reverse 0..Size(n)-1 loop
            res := res*the_base + Coefficient(n,i);
          end loop;
    end if;
    return res;
  end Create;

-- SELECTORS :

  function Base return natural is
  begin
    return the_base;
  end Base;

  function Exponent return natural is
  begin
    return expo;
  end Exponent;

  function Empty ( n : Natural_Number ) return boolean is
  begin
    return ( n = null );
  end Empty;

  function Size ( n : Natural_Number ) return natural is
  begin
    if n = null
     then return 0;
     else return n.size;
    end if;
  end Size;

  function Decimal_Places ( n : natural ) return natural is

    acc : natural := 1;

  begin
    for i in 0..n loop
      if n < acc
       then return i;
       else acc := acc*10;
      end if;
    end loop;
    return n;
  end Decimal_Places;

  function Decimal_Places ( n : Natural_Number ) return natural is
  
    res : natural := 0;

  begin
    if not Empty(n)
     then for i in reverse 0..Size(n) loop
            if n.numb(i) /= 0
             then res := Decimal_Places(n.numb(i)) + i*expo;
                  exit;
            end if;
          end loop;
    end if;
    return res;
  end Decimal_Places;

  function Coefficient ( n : Natural_Number; i : natural ) return natural is
  begin
    if (n = null or else i > n.size)
     then return 0;
     else return n.numb(i);
    end if;
  end Coefficient;

  function Coefficients ( n : Natural_Number ) return Array_of_Naturals is

    na : Array_of_Naturals(0..0) := (0..0 => 0);

  begin
    if Empty(n)
     then return na;
     else return n.numb;
    end if;
  end Coefficients;

-- COMPARISON AND COPYING :
--   Note that these implementations are all independent of the
--   actual representation of Natural_Number as they use the selectors.

  function Equal ( n1,n2 : Array_of_Naturals ) return boolean is
  begin
    if n1'first /= n2'first
     then return false;
     elsif n1'last /= n2'last
         then return false;
         else for i in n1'range loop
                if n1(i) /= n2(i)
                 then return false;
                end if;
              end loop;
              return true;
    end if;
  end Equal;

  function Equal ( n1 : Natural_Number; n2 : natural ) return boolean is
  begin
    if Empty(n1)
     then if n2 = 0
           then return true;
           else return false;
          end if;
     else declare
            n2cff : constant Array_of_Naturals := Create(n2);
          begin
            if n2cff'last > Size(n1)
             then return false;
             else for i in n2cff'range loop
                    if Coefficient(n1,i) /= n2cff(i)
                     then return false;
                    end if;
                  end loop;
                  for i in n2cff'last+1..Size(n1) loop
                    if Coefficient(n1,i) /= 0
                     then return false;
                    end if;
                  end loop;
            end if;
          end;
          return true;
    end if;
  end Equal;

  function Equal ( n1,n2 : Natural_Number ) return boolean is

    min_size : natural;

  begin
    if Empty(n1)
     then return Equal(n2,0);
     else if Empty(n2)
           then return Equal(n1,0);
           else if Size(n1) < Size(n2)
                 then for i in Size(n1)+1..Size(n2) loop
                        if Coefficient(n2,i) /= 0
                         then return false;
                        end if;
                      end loop;
                      min_size := Size(n1);
                 elsif Size(n1) > Size(n2)
                     then for i in Size(n2)+1..Size(n1) loop
                            if Coefficient(n1,i) /= 0
                             then return false;
                            end if;
                          end loop;
                          min_size := Size(n2);
                     else min_size := Size(n1);
                end if;
                return Equal(Coefficients(n1)(0..min_size),
                             Coefficients(n2)(0..min_size));
          end if;
    end if;
  end Equal;

  function "<" ( n1 : Natural_Number; n2 : natural ) return boolean is
  begin
    if Empty(n1)
     then if n2 > 0
           then return true;
           else return false;
          end if;
     else declare
            n2cff : constant Array_of_Naturals := Create(n2);
          begin
            if n2cff'last > Size(n1)
             then return true;
             else if n2cff'last < Size(n1) 
                   then for i in reverse n2cff'last+1..Size(n1) loop
                          if Coefficient(n1,i) /= 0
                           then return false;
                          end if;
                        end loop;
                  end if;
                  for i in reverse n2cff'range loop
                    if Coefficient(n1,i) >= n2cff(i)
                     then return false;
                    end if;
                  end loop;
            end if;
          end;
          return true;
    end if;
  end "<";

  function "<" ( n1 : natural; n2 : Natural_Number ) return boolean is
  begin
    if Empty(n2)
     then return false;
     else declare
            n1cff : constant Array_of_Naturals := Create(n1);
          begin
            if n1cff'last > Size(n2)
             then return false;
             else if n1cff'last < Size(n2)
                   then for i in reverse n1cff'last+1..Size(n2) loop
                          if Coefficient(n2,i) /= 0
                           then return true;
                          end if;
                        end loop;
                  end if;
                  for i in reverse n1cff'range loop
                    if n1cff(i) >= Coefficient(n2,i)
                     then return false;
                    end if;
                  end loop;
            end if;
          end;
          return true;
    end if;
  end "<";

  function "<" ( n1,n2 : Natural_Number ) return boolean is

    min_size : natural;

  begin
    if Empty(n1)
     then if Empty(n2)
           then return false;
           else return true;
          end if;
     else if Empty(n2)
           then return false;
           else if Size(n1) < Size(n2)
                 then for i in Size(n1)+1..Size(n2) loop
                        if Coefficient(n2,i) /= 0
                         then return true;
                        end if;
                      end loop;
                      min_size := Size(n1);
                 elsif Size(n1) > Size(n2)
                     then for i in Size(n2)+1..Size(n1) loop
                            if Coefficient(n1,i) /= 0
                             then return false;
                            end if;
                          end loop;
                          min_size := Size(n2);
                     else min_size := Size(n1);
                end if;
                for i in reverse 0..min_size loop
                  if Coefficient(n1,i) < Coefficient(n2,i)
                   then return true;
                   elsif Coefficient(n1,i) > Coefficient(n2,i)
                       then return false;
                  end if;
                end loop;
                return false;      -- n1 = n2
          end if;
    end if;
  end "<";

  function ">" ( n1 : Natural_Number; n2 : natural ) return boolean is
  begin
    if Empty(n1)
     then return false;
     else declare
            n2cff : constant Array_of_Naturals := Create(n2);
          begin
            if n2cff'last > Size(n1)
             then return false;
             else if n2cff'last < Size(n1)
                   then for i in n2cff'last+1..Size(n1) loop
                          if Coefficient(n1,i) /= 0
                           then return true;
                          end if;
                        end loop;
                  end if;
                  for i in reverse n2cff'range loop
                    if Coefficient(n1,i) <= n2cff(i)
                     then return false;
                    end if;
                  end loop;
            end if;
          end;
          return true;
    end if;
  end ">";
 
  function ">" ( n1 : natural; n2 : Natural_Number ) return boolean is
  begin
    if Empty(n2)
     then if n1 > 0
           then return true;
           else return false;
          end if;
     else declare
            n1cff : constant Array_of_Naturals := Create(n1);
          begin
            if n1cff'last > Size(n2)
             then return true;
             else if n1cff'last < Size(n2)
                   then for i in n1cff'last+1..Size(n2) loop
                          if Coefficient(n2,i) /= 0
                           then return false;
                          end if;
                        end loop;
                  end if;
                  for i in reverse n1cff'range loop
                    if n1cff(i) <= Coefficient(n2,i)
                     then return false;
                    end if;
                  end loop;
            end if;
          end;
          return true;
    end if;
  end ">";

  function ">" ( n1,n2 : Natural_Number ) return boolean is

    min_size : natural;

  begin
    if Empty(n2)
     then if Empty(n1)
           then return false;
           else return true;
          end if;
     else if Empty(n1)
           then return false;
           else if Size(n1) < Size(n2)
                 then for i in Size(n1)+1..Size(n2) loop
                        if Coefficient(n2,i) /= 0
                         then return false;
                        end if;
                      end loop;
                      min_size := Size(n1);
                 elsif Size(n1) > Size(n2)
                     then for i in Size(n2)+1..Size(n1) loop
                            if Coefficient(n1,i) /= 0
                             then return true;
                            end if;
                          end loop;
                          min_size := Size(n2);
                     else min_size := Size(n1);
                end if;
                for i in reverse 0..min_size loop
                  if Coefficient(n1,i) > Coefficient(n2,i)
                   then return true;
                   elsif Coefficient(n1,i) < Coefficient(n2,i)
                       then return false;
                  end if;
                end loop;
                return false;  -- n1 = n2
          end if;
    end if;
  end ">";

  procedure Copy ( n1 : in natural; n2 : in out Natural_Number ) is
  begin
    Clear(n2);
    n2 := Create(n1);
  end Copy;

  procedure Copy ( n1 : in Natural_Number; n2 : in out Natural_Number ) is
  begin
    Clear(n2);
    if not Empty(n1)
     then declare
            n2rep : Natural_Number_Rep(Size(n1));
          begin
            for i in n2rep.numb'range loop
              n2rep.numb(i) := n1.numb(i);
            end loop;
            n2 := new Natural_Number_Rep'(n2rep);
          end;
    end if;
  end Copy;

-- AUXILIARIES FOR ARITHMETIC OPERATIONS :

  function Extend ( n : Array_of_Naturals; carry : natural )
                  return Natural_Number_Rep is

  -- DESCRIPTION :
  --   Extends the array of naturals as much as need to store the carry
  --   obtained from adding a number to a natural number.

  begin
    if carry = 0
     then declare
            nrep : Natural_Number_Rep(n'last);
          begin
            nrep.numb := n;
            return nrep;
          end;
     else declare
            np1 : Array_of_Naturals(n'first..n'last+1);
            carry1 : natural;
          begin
            np1(n'range) := n;
            np1(np1'last) := carry mod the_base;
            carry1 := carry/the_base;
            return Extend(np1,carry1);
          end;
    end if;
  end Extend;

  procedure Extend ( n : in out Natural_Number; carry : in natural ) is

  -- DESCRIPTION :
  --   Extends the coefficient vector of the natural number to store
  --   the carry-over.
 
  begin
    if carry /= 0
     then declare
            res : Natural_Number
                := new Natural_Number_Rep'(Extend(n.numb,carry));
          begin
            Clear(n);
            n := res;
          end;
    end if;
  end Extend;

  procedure Mul_Fact ( n : in out Natural_Number; f : in natural ) is

  -- DESCRIPTION :
  --   Multiplies the number n with 0 < f <= fact and stores the result in n.

    prod,carry : natural;

  begin
    if not Empty(n)
     then carry := 0;
          for i in n.numb'range loop
            prod := n.numb(i)*f + carry;
            n.numb(i) := prod mod the_base;
            carry := prod/the_base;
          end loop;
          Extend(n,carry);
    end if;
  end Mul_Fact;

  procedure Mul_Base ( n : in out Natural_Number; k : in natural ) is

  -- DESCRIPTION :
  --   Multiplies the number n with the_base**k.

  begin
    if k > 0 and not Empty(n)
     then declare
            npk : Natural_Number_Rep(n.size+k);
          begin
            npk.numb(0..k-1) := (0..k-1 => 0);
            npk.numb(k..npk.numb'last) := n.numb(n.numb'range);
            Clear(n);
            n := new Natural_Number_Rep'(npk);
          end;
    end if;
  end Mul_Base;

  procedure Div_Base ( n : in out Natural_Number; k : in natural ) is

  -- DESCRIPTION :
  --   Divides the number n by the_base**k.

  begin
    if k > 0 and not Empty(n)
     then if k > Size(n)
           then Clear(n);
           else declare
                  nmk : Natural_Number_Rep(n.size-k);
                begin
                  for i in nmk.numb'range loop
                    nmk.numb(i) := n.numb(k+i);
                  end loop;
                  Clear(n);
                  n := new Natural_Number_Rep'(nmk);
                end;
          end if;
    end if;
  end Div_Base;

-- ARITHMETIC OPERATIONS :

  function "+" ( n1 : Natural_Number; n2 : natural ) return Natural_Number is

    res : Natural_Number;

  begin
    if Empty(n1)
     then res := Create(n2);
     else declare
            cff : Array_of_Naturals(n1.numb'range);
            sum,carry : natural;
          begin
            carry := n2;
            for i in cff'range loop
              sum := n1.numb(i) + carry;
              cff(i) := sum mod the_base;
              carry := sum/the_base;
            end loop;
            res := new Natural_Number_Rep'(Extend(cff,carry));
          end;
    end if;
    return res;
  end "+";

  function "+" ( n1 : natural; n2 : Natural_Number ) return Natural_Number is
  begin
    return n2+n1;
  end "+";

  function "+" ( n1,n2 : Natural_Number ) return Natural_Number is

    res : Natural_Number;
    sum,carry : natural;

  begin
    if Empty(n1)
     then Copy(n2,res);
     elsif Empty(n2)
         then Copy(n1,res);
         else carry := 0;
              if Size(n1) > Size(n2)
               then declare
                      res_rep : Natural_Number_Rep(Size(n1)) := n1.all;
                    begin
                      for i in n2.numb'range loop
                        sum := res_rep.numb(i) + n2.numb(i) + carry;
                        res_rep.numb(i) := sum mod the_base;
                        carry := sum/the_base;
                      end loop;
                      if carry /= 0
                       then for i in n2.numb'last+1..n1.numb'last loop
                              sum := res_rep.numb(i) + carry;
                              res_rep.numb(i) := sum mod the_base;
                              carry := sum/the_base;
                              exit when (carry = 0);
                            end loop;
                      end if;
                      res := new Natural_Number_Rep'(res_rep);
                    end;
               else declare
                      res_rep : Natural_Number_Rep(Size(n2)) := n2.all;
                    begin
                      for i in n1.numb'range loop
                        sum := res_rep.numb(i) + n1.numb(i) + carry;
                        res_rep.numb(i) := sum mod the_base;
                        carry := sum/the_base;
                      end loop;
                      if carry /= 0
                       then for i in n1.numb'last+1..n2.numb'last loop
                              sum := res_rep.numb(i) + carry;
                              res_rep.numb(i) := sum mod the_base;
                              carry := sum/the_base;
                              exit when (carry = 0);
                            end loop;
                      end if;
                      res := new Natural_Number_Rep'(res_rep);
                    end;
             end if;
             Extend(res,carry);
    end if;
    return res;
  end "+";

  procedure Add ( n1 : in out Natural_Number; n2 : in natural ) is

    sum,carry : natural;

  begin
    if Empty(n1)
     then n1 := Create(n2);
     else carry := n2;
          for i in n1.numb'range loop
            sum := n1.numb(i) + carry;
            n1.numb(i) := sum mod the_base;
            carry := sum/the_base;
          end loop;
          Extend(n1,carry);
    end if;
  end Add;

  procedure Add ( n1 : in out Natural_Number; n2 : in Natural_Number ) is

    res : Natural_Number;
    sum,carry : natural;

  begin
    if Empty(n1)
     then Copy(n2,n1);
     elsif not Empty(n2)
         then carry := 0;
              if Size(n1) >= Size(n2)
               then for i in n2.numb'range loop
                      sum := n1.numb(i) + n2.numb(i) + carry;
                      n1.numb(i) := sum mod the_base;
                      carry := sum/the_base;
                    end loop;
                    if carry /= 0
                     then for i in n2.numb'last+1..n1.numb'last loop
                            sum := n1.numb(i) + carry;
                            n1.numb(i) := sum mod the_base;
                            carry := sum/the_base;
                            exit when (carry = 0);
                          end loop;
                          Extend(n1,carry);
                    end if;
               else declare
                      res_rep : Natural_Number_Rep(Size(n2)) := n2.all;
                    begin
                      for i in n1.numb'range loop
                        sum := res_rep.numb(i) + n1.numb(i) + carry;
                        res_rep.numb(i) := sum mod the_base;
                        carry := sum/the_base;
                      end loop;
                      if carry /= 0
                       then for i in n1.numb'last+1..n2.numb'last loop
                              sum := res_rep.numb(i) + carry;
                              res_rep.numb(i) := sum mod the_base;
                              carry := sum/the_base;
                              exit when (carry = 0);
                            end loop;
                      end if;
                      Clear(n1);
                      n1 := new Natural_Number_Rep'(res_rep);
                      if carry /= 0
                       then Extend(n1,carry);
                      end if;
                    end;
              end if;
    end if;
  end Add;

  function "-" ( n1 : Natural_Number; n2 : natural ) return Natural_Number is

    res : Natural_Number;
    diff,carry : integer;

  begin
    if not Empty(n1)
     then Copy(n1,res);
          declare
            n2cff : constant Array_of_Naturals := Create(n2);
            index : natural := n2cff'first;
          begin
            carry := n2cff(index);
            for i in n1.numb'range loop
              diff := n1.numb(i) - carry;
              if diff >= 0
               then res.numb(i) := diff mod the_base;
                    carry := diff/the_base;
               else diff := the_base + diff;
                    res.numb(i) := diff mod the_base;
                    carry := diff/the_base + (res.numb(i)+carry)/the_base;
              end if;
              if index < n2cff'last
               then index := index+1;
                    carry := carry + n2cff(index);
              end if;
              exit when (carry = 0);
            end loop;
          end;
    end if;
    return res;
  end "-";

  function "-" ( n1 : natural; n2 : Natural_Number ) return Natural_Number is
  
    tmp : Natural_Number := Create(n1);
    res : Natural_Number := tmp - n2;

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

  function "-" ( n1,n2 : Natural_Number ) return Natural_Number is

    res,diff : Natural_Number;

  begin
    Copy(n1,res);
    if not Empty(n2)
     then for i in reverse n2.numb'range loop
            Copy(res,diff);
            Div_Base(diff,i);
            if not Empty(diff)
             then Sub(diff,n2.numb(i));
                  for j in diff.numb'range loop
                    res.numb(i+j) := diff.numb(j);
                  end loop;
                  Clear(diff);
            end if;
          end loop;
    end if;
    return res;
  end "-";

  function "+" ( n : Natural_Number ) return Natural_Number is

    res : Natural_Number;

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

  function "-" ( n : Natural_Number ) return Natural_Number is

    res : Natural_Number;

  begin
    Copy(n,res);
    return res;
  end "-";

  procedure Min ( n : in out Natural_Number ) is
  begin
    null;
  end Min;

  procedure Sub ( n1 : in out Natural_Number; n2 : in natural ) is

    diff,carry : integer;

  begin
    if not Empty(n1) and n2 > 0
     then declare
            n2cff : constant Array_of_Naturals := Create(n2);
            index : natural := n2cff'first;
          begin
            carry := n2cff(index);
            for i in n1.numb'range loop
              diff := n1.numb(i) - carry;
              if diff >= 0
               then n1.numb(i) := diff mod the_base;
                    carry := diff/the_base;
               else diff := the_base + diff;
                    n1.numb(i) := diff mod the_base;
                    carry := diff/the_base + (n1.numb(i)+carry)/the_base;
              end if;
              if index < n2cff'last
               then index := index+1;
                    carry := carry + n2cff(index);
              end if;
              exit when (carry = 0);
            end loop;
          end;
    end if;
  end Sub;

  procedure Sub ( n1 : in out Natural_Number; n2 : in Natural_Number ) is

    res,diff : Natural_Number;

  begin
    if not Empty(n1) and not Empty(n2)
     then for i in reverse n2.numb'range loop
            Copy(n1,diff);
            Div_Base(diff,i);
            if not Empty(diff)
             then Sub(diff,n2.numb(i));
                  for j in diff.numb'range loop
                    n1.numb(i+j) := diff.numb(j);
                  end loop;
                  Clear(diff);
            end if;
          end loop;
    end if;
  end Sub;

  function "*" ( n1 : Natural_Number; n2 : natural ) return Natural_Number is

    res,prod,nres : Natural_Number;
    remainder,modfact : natural;
    cnt : natural := 0;

  begin
    if n2 /= 0
     then remainder := n2;
          while remainder /= 0 loop
            modfact := remainder mod fact;
            if modfact /= 0
             then Copy(n1,prod);
                  Mul_Fact(prod,modfact);
                  for i in 1..cnt loop
                    Mul_Fact(prod,fact);
                  end loop;
                  Add(res,prod);
                  Clear(prod);
            end if;
            cnt := cnt+1;
            remainder := remainder/fact;
          end loop;
    end if;
    return res;
  end "*";

  function "*" ( n1 : natural; n2 : Natural_Number ) return Natural_Number is
  begin
    return n2*n1;
  end "*";

  function "*" ( n1,n2 : Natural_Number ) return Natural_Number is

    res,prod : Natural_Number;

  begin
    if (Empty(n1) or else Empty(n2))
     then return res;
     else if Size(n1) >= Size(n2)
           then for i in n2.numb'range loop
                  prod := n1*n2.numb(i);
                  Mul_Base(prod,i);
                  Add(res,prod);
                  Clear(prod);
                end loop;
           else for i in n1.numb'range loop
                  prod := n2*n1.numb(i);
                  Mul_Base(prod,i);
                  Add(res,prod);
                  Clear(prod);
                end loop;
          end if;
    end if;
    return res;
  end "*";

  procedure Mul ( n1 : in out Natural_Number; n2 : in natural ) is

    res,prod : Natural_Number;
    remainder,modfact : natural;
    cnt : natural := 0;

  begin
    if n2 = 0
     then Clear(n1);
     else remainder := n2;
          while remainder /= 0 loop
            modfact := remainder mod fact;
            if modfact /= 0
             then Copy(n1,prod);
                  Mul_Fact(prod,modfact);
                  for i in 1..cnt loop
                    Mul_Fact(prod,fact);
                  end loop;
                  Add(res,prod);
                  Clear(prod);
            end if;
            cnt := cnt+1;
            remainder := remainder/fact;
          end loop;
          Clear(n1);
          n1 := res; 
    end if;
  end Mul;

  procedure Mul ( n1 : in out Natural_Number; n2 : in Natural_Number ) is
 
    res,prod : Natural_Number;

  begin
    if not Empty(n1)
     then if Empty(n2)
           then Clear(n1);
           else if Size(n1) >= Size(n2)
                 then for i in n2.numb'range loop
                        prod := n1*n2.numb(i);
                        Mul_Base(prod,i);
                        Add(res,prod);
                        Clear(prod);
                      end loop;
                 else for i in n1.numb'range loop
                        prod := n2*n1.numb(i);
                        Mul_Base(prod,i);
                        Add(res,prod);
                        Clear(prod);
                      end loop;
                end if;
                Clear(n1);
                n1 := res;
          end if;
    end if;
  end Mul;

  function "**" ( n1 : Natural_Number; n2 : natural ) return Natural_Number is

    res : Natural_Number;

  begin
    if n2 = 0
     then res := Create(1);
     else if not Empty(n1)
           then Copy(n1,res);
                for i in 1..n2-1 loop
                  Mul(res,n1);
                end loop;
          end if;
    end if;
    return res;
  end "**";

  function "**" ( n1 : natural; n2 : Natural_Number ) return Natural_Number is

    res,cnt : Natural_Number;

  begin
    if Equal(n2,0)
     then res := Create(1);
     else res := Create(n1);
          if n1 /= 0
           then cnt := Create(1);
                while not Equal(cnt,n2) loop
                  Mul(res,n1);
                  Add(cnt,1);
                end loop;
                Clear(cnt);
          end if;
    end if;
    return res;
  end "**";

  function "**" ( n1,n2 : Natural_Number ) return Natural_Number is

    res,cnt : Natural_Number;

  begin
    if Equal(n2,0)
     then res := Create(1);
     else if not Empty(n1)
           then Copy(n1,res);
                cnt := Create(1);
                while not Equal(cnt,n2) loop
                  Mul(res,n1);
                  Add(cnt,1);
                end loop;
                Clear(cnt);
          end if;
    end if;
    return res;
  end "**";

  procedure Small_Div ( n1 : in Natural_Number; n2 : in natural;
                        q : out Natural_Number; r : out natural ) is

  -- DESCRIPTION :
  --   n1 = n2*q + r, only applies when n2 <= fact.

   -- res : Natural_Number;
    carry : natural := 0;

  begin
    if n2 /= 0
     then if not Empty(n1)
           then -- res := new Natural_Number_Rep(Size(n1));
                q := new Natural_Number_Rep(Size(n1));
                for i in reverse 1..Size(n1) loop
                 -- res.numb(i) := (n1.numb(i)+carry)/n2;
                  q.numb(i) := (n1.numb(i)+carry)/n2;
                  carry := (n1.numb(i)+carry) mod n2;
                  carry := carry*the_base;
                end loop;
               -- res.numb(0) := (n1.numb(0)+carry)/n2;
                q.numb(0) := (n1.numb(0)+carry)/n2;
                carry := (n1.numb(0)+carry) mod n2;
          end if;
     else raise NUMERIC_ERROR;
    end if;
   -- q := res;
    r := carry;
  end Small_Div;

  procedure Small_Div ( n1 : in out Natural_Number; n2 : in natural;
                        r : out natural ) is

    q : Natural_Number;

  begin
    Small_Div(n1,n2,q,r);
    Copy(q,n1); Clear(q);
  end Small_Div;

  procedure Big_Div ( n1 : in Natural_Number; n2 : in natural;
                      q : out Natural_Number; r : out natural ) is

  -- DESCRIPTION :
  --   This division has to be applied when n2 > fact.

    wrkn1,shiftn1,remn1,acc,qres,qsub : Natural_Number;
    cntshf,cntsub,rres : natural;

  begin
    if n2 /= 0
     then
       if n1 < n2
        then 
          q := qres;
          r := Create(n1);
        else
          Copy(n1,shiftn1);                  -- we have n1 >= n2
          Small_Div(shiftn1,fact,wrkn1,rres); 
          cntshf := 0;
          while not (wrkn1 < n2) loop        -- invariant n2 <= shiftn1
            acc := Create(rres);
            for i in 1..cntshf loop
              Mul_Fact(acc,fact);
            end loop;
            Add(remn1,acc); Clear(acc);
            Copy(wrkn1,shiftn1);
            cntshf := cntshf + 1;
            Small_Div(wrkn1,fact,rres);      -- at end of loop
          end loop;                          --   shiftn1/fact < n2 <= shiftn1
          Clear(wrkn1);
          cntsub := 0;
          while not (shiftn1 < n2) loop      -- subtract n2 from shiftn1
            cntsub := cntsub + 1;
            Sub(shiftn1,n2);
          end loop;
          qsub := Create(cntsub);            -- intermediate quotient
          for i in 1..cntshf loop
            Mul_Fact(qsub,fact);
            Mul_Fact(shiftn1,fact);
          end loop;
          Add(shiftn1,remn1);
          Clear(remn1);
          Div(shiftn1,n2,qres,r);            -- recursive call
          Add(qsub,qres);                    -- collect results
          q := qsub;
          Clear(qres); Clear(shiftn1);       -- cleaning up
       end if;
     else
       raise NUMERIC_ERROR;
    end if;
  end Big_Div;

  function "/" ( n1 : Natural_Number; n2 : natural ) return Natural_Number is

    res : Natural_Number;
    r : natural;

  begin
    if not Empty(n1)
     then if n2 <= fact
           then Small_Div(n1,n2,res,r);
           else Big_Div(n1,n2,res,r);
          end if;
    end if;
    return res;
  end "/";

  function "/" ( n1 : natural; n2 : Natural_Number ) return natural is

    res : natural;

  begin
    if n1 < n2
     then res := 0;
     elsif not Empty(n2)
         then res := Create(n2);
              res := n1/res;
         else raise NUMERIC_ERROR;
    end if;
    return res;
  end "/";

  function "/" ( n1,n2 : Natural_Number ) return Natural_Number is

    res,rrr : Natural_Number;

  begin 
    Div(n1,n2,res,rrr);
    Clear(rrr);
    return res;
  end "/";

  function Rmd ( n1 : Natural_Number; n2 : natural ) return natural is

    res : natural;
    q : Natural_Number;

  begin
    if n2 <= fact
     then Small_Div(n1,n2,q,res);
     else Big_Div(n1,n2,q,res);
    end if;
    Clear(q);
    return res;
  end Rmd;

  function Rmd ( n1 : natural; n2 : Natural_Number ) return natural is

    res : natural;

  begin
    if n1 < n2
     then res := n1;
     elsif not Empty(n2)
         then res := Create(n2);
              res := n1 mod res;
         else raise NUMERIC_ERROR;
    end if;
    return res;
  end Rmd;

  function Rmd ( n1,n2 : Natural_Number ) return Natural_Number is

    res,q : Natural_Number;

  begin
    Div(n1,n2,q,res);
    Clear(q);
    return res;
  end Rmd;

  procedure Div ( n1 : in out Natural_Number; n2 : in natural ) is

    r : natural;

  begin
    Div(n1,n2,r);
  end Div;

  procedure Div ( n1 : in out Natural_Number; n2 : in Natural_Number ) is

    r : Natural_Number;

  begin
    Div(n1,n2,r);
    Clear(r);
  end Div;

  procedure Div ( n1 : in Natural_Number; n2 : in natural;
                  q : out Natural_Number; r : out natural ) is
  begin
    if n2 <= fact
     then Small_Div(n1,n2,q,r);
     else Big_Div(n1,n2,q,r);
    end if;
  end Div;

  procedure Div ( n1 : in out Natural_Number; n2 : in natural;
                  r : out natural ) is

    q : Natural_Number;

  begin
    Div(n1,n2,q,r);
    Copy(q,n1); Clear(q);
  end Div;

  procedure Div ( n1,n2 : in Natural_Number;
                  q : out Natural_Number; r : out Natural_Number ) is

    wrkn1,shiftn1,remn1,acc,qres,qsub : Natural_Number;
    cntshf,cntsub,rres : natural;

  begin
    if (Empty(n2) or else Equal(n2,0))
     then 
       raise NUMERIC_ERROR;
     else
       if n1 < n2
        then 
          q := qres;
          Copy(n1,r);
        else
          Copy(n1,shiftn1);                  -- we have n1 >= n2
          Small_Div(shiftn1,fact,wrkn1,rres);
          cntshf := 0;
          while not (wrkn1 < n2) loop        -- invariant n2 <= shiftn1
            acc := Create(rres);
            for i in 1..cntshf loop
              Mul_Fact(acc,fact);
            end loop;
            Add(remn1,acc); Clear(acc);
            Copy(wrkn1,shiftn1);
            cntshf := cntshf + 1;
            Small_Div(wrkn1,fact,rres);      -- at end of loop
          end loop;                          --   shiftn1/fact < n2 <= shiftn1
          Clear(wrkn1);
          cntsub := 0;
          while not (shiftn1 < n2) loop      -- subtract n2 from shiftn1
            cntsub := cntsub + 1;
            Sub(shiftn1,n2);
          end loop;
          qsub := Create(cntsub);            -- intermediate quotient
          for i in 1..cntshf loop
            Mul_Fact(qsub,fact);
            Mul_Fact(shiftn1,fact);
          end loop;
          Add(shiftn1,remn1);
          Clear(remn1);
          Div(shiftn1,n2,qres,r);            -- recursive call
          Add(qsub,qres);                    -- collect results
          q := qsub;
          Clear(qres); Clear(shiftn1);       -- cleaning up
       end if;
    end if;
  end Div;

  procedure Div ( n1 : in out Natural_Number; n2 : in Natural_Number;
                  r : out Natural_Number ) is

    q : Natural_Number;

  begin
    Div(n1,n2,q,r);
    Copy(q,n1); Clear(q);
  end Div;

-- DESTRUCTOR :

  procedure Clear ( n : in out Natural_Number ) is

    procedure free is
      new unchecked_deallocation(Natural_Number_Rep,Natural_Number);

  begin
    if not Empty(n)
     then free(n);
    end if;
  end Clear;

end Multprec_Natural_Numbers;