Removed MIT licensed Mathpaqs packages in favour of bindings to GMP

1 files changed, 0 insertions, 1530 deletions

diff --git a/src/multi_precision_integers.adb b/src/multi_precision_integers.adb
deleted file mode 100644
index 70a6d28..0000000
--- a/src/multi_precision_integers.adb
+++ /dev/null
@@ -1,1530 +0,0 @@
------------------------------------------------------------------------------

---  File: mupreint.adb; see specification (mupreint.ads)

------------------------------------------------------------------------------

---    Aug-2007: - No more generics (Long_Block_type,

---                Block_type,... always the largest possible idea: J.C.)

---              - Fixed Basic(...) (based on J.C.'s remarks)

---    Nov-2006: - Multiply_internal with/without copy of result (automatic

---                detection of when it is needed)

---              - Explicit Multiply_internal for Multi_int * Basic_int

---              - Multiply(multi,basic,multi) available as procedure

---              - useless zeroing of quotient removed

---              - useless zeroing of blocks removed for indices

---                above last possible used in *

--- 24-Feb-2002: Div_Rem(u, v, v, r) also possible

--- 23-Feb-2002: DEBUG: +: multiplications are verified by dividing the result

---                     +: divisions are verified by comparing i2*q+r and i1

--- 15-Feb-2002: "zero" and 1st index in Divide_absolute_normalized

---                 bugs fixed by Duncan Sands (D.S.)

-

--- To-do/bug symbol: !!

-

-

-

-

-------------------------------------------------------------------------------

---

---  Copyright (c) 2007 G. de Montmollin, Univ. Neuchatel

---

---  Permission is hereby granted, free of charge, to any person obtaining a copy

---  of this software and associated documentation files (the "Software"), to deal

---  in the Software without restriction, including without limitation the rights

---  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell

---  copies of the Software, and to permit persons to whom the Software is

---  furnished to do so, subject to the following conditions:

---

---  The above copyright notice and this permission notice shall be included in all

---  copies or substantial portions of the Software.

---

---  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR

---  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

---  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE

---  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER

---  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,

---  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE

---  SOFTWARE.

---

-------------------------------------------------------------------------------

-

-

-

-

-with Multi_precision_integers.Check;

--- with Ada.Text_IO;

-with Ada.Unchecked_Deallocation;

-

-package body Multi_precision_integers is

-

-  function Shift_Left

-    (Value  : Block_type;

-     Amount : Natural) return Block_type;

-

-  function Shift_Right

-    (Value  : Block_type;

-     Amount : Natural) return Block_type;

-

-  function Shift_Left

-    (Value  : Long_Block_type;

-     Amount : Natural) return Long_Block_type;

-

-  function Shift_Right

-    (Value  : Long_Block_type;

-     Amount : Natural) return Long_Block_type;

-

-  pragma Import (Intrinsic, Shift_Left);

-  pragma Import (Intrinsic, Shift_Right);

-

-  package Check_internal renames Multi_precision_integers.Check;

-

-  -- Internal_error: exception;

-  -- Not_done: exception;

-

-  type compar is (smaller, equal, greater);

-

-  function Min (a,b: Index_int) return Index_int is

-  begin if a < b then return a; else return b; end if; end Min;

-

-  function Max (a,b: Index_int) return Index_int is

-  begin if a > b then return a; else return b; end if; end Max;

-

-  procedure Reduce_last_nonzero ( n: in out Multi_int ) is

-    old_last : constant Index_int:= n.last_used;

-  begin

-    if Debug then Check_internal.Test(n, test_last=> False); end if;

-

-    if n.zero then -- We avoid de-zeroing accidentally

-      return;      -- and returning a false non-zero with rubbish :-)

-    end if;

-

-    n.zero := True;

-    for i in 0 .. old_last loop -- after old_last it *is* rubbish anyway.

-      if n.blk(i) /= 0 then

-        n.zero := False;

-        n.last_used := i;

-      end if;

-    end loop;

-  end Reduce_last_nonzero;

-

-  function Compare_absolute (i1, i2: Multi_int) return compar is

-    l1, l2 : Index_int;

-  begin

-    -- On ne compare que ABS(i1) et ABS(i2)

-    l1:= i1.last_used;

-    l2:= i2.last_used;

-    if l1 > l2 then         -- i1 a plus de blocs non nuls

-      return greater;

-    elsif l1 < l2 then      -- i1 a moins de blocs non nuls

-      return smaller;

-    else                       -- i1 et i2 ont le meme nb de blocs

-      for i in reverse 0 .. l1 loop -- on parcourt du + signifiant au -

-        if    i1.blk(i) > i2.blk(i) then -- <<chiffre>> de i1 plus grand

-          return greater;

-        elsif i1.blk(i) < i2.blk(i) then -- <<chiffre>> de i1 plus petit

-          return smaller;

-        end if;

-        -- M\^emes chiffres -> au suivant!

-      end loop;

-      -- Bon, les 2 nombres sont egaux!

-      return equal;

-    end if;

-  end Compare_absolute;

-

-  ----- Informations, conversions

-

-  function Multi(small: Basic_int) return Multi_int is

-    abss: constant Long_Block_type:= Long_Block_type(abs small);

-    reste: Long_Block_type;

-    negs: constant Boolean:= small < 0;

-    Conversion_overflow : exception;

-

-  begin

-

-    if abss <= Long_Block_type(maxblock) then

-      return Multi_int'

-             ( n=>         0,          -- 1 bloc suffit

-               blk=>      (0=> Block_type(abss)),  -- le bloc contient le nombre

-               neg=>       negs,

-               zero=>      small = 0,

-               last_used=> 0

-             );

-    else

-      reste:= Shift_Right(abss, Block_type_bits);

-      if reste <= Long_Block_type(maxblock) then

-        return ( n=>         1,            -- il faut 2 blocs

-                 blk=>      (0=> Block_type(abss and maxblock), -- bloc 0

-                             1=> Block_type(reste)),  -- bloc 1

-                 neg=>       negs,

-                 zero=>      False,

-                 last_used=> 1

-               );

-      else

-        if Shift_Right(reste, Block_type_bits) > Long_Block_type(maxblock) then

-           raise Conversion_overflow;

-        end if;

-

-        return ( n=>     2,  -- il faut 3 blocs (e.g. 31 bits 15+15+1)

-                 blk=>  (0=> Block_type(abss and maxblock),   -- bloc 0

-                         1=> Block_type(reste and maxblock),  -- bloc 1

-                         2=> Block_type(Shift_Right(reste, Block_type_bits))  -- bloc 2

-                         ),

-                 neg=>   negs,

-                 zero=>  False,

-                 last_used=> 2

-               );

-      end if;

-    end if;

-  end Multi;

-

-  zero: constant Multi_int:= Multi(0);

-  one : constant Multi_int:= Multi(1);

-

-  Blocks_Per_Basic : constant Positive :=

-     (Basic_int'Size + Block_type'Size - 1) / Block_type'Size;

-

-  -- Convert Multi_int to Basic_int (when possible, else: Cannot_fit raised)

-  -- 2007:

-  --  - correct code for block sizes smaller than Basic_int

-  --  - fixed usage of negative flag

-  function Basic(large: Multi_int) return Basic_int is

-    type Same_as_Basic_natural is mod 2 ** Basic_int'Size;

-    function Shift_Left

-      (Value  : Same_as_Basic_natural;

-       Amount : Natural) return Same_as_Basic_natural;

-    pragma Import (Intrinsic, Shift_Left);

-    result : Same_as_Basic_natural;

-    block_value: Block_type;

-    type Huge_int is mod System.Max_Binary_Modulus;

-    last_bit: Natural;

-  begin

-    if large.zero then -- <- 17-Feb-2002

-      return 0;

-    end if;

-    -- Case: too many blocks (whatever sizes)

-    if 1 + large.last_used > Blocks_Per_Basic then

-      raise Cannot_fit;

-    end if;

-    -- Case: block size and contents larger than basic

-    block_value:= large.blk(large.last_used);

-    if Huge_int(block_value) > Huge_int(Basic_int'Last) then

-      raise Cannot_fit;

-    end if;

-    declare

-      tmp: Block_type:= block_value;

-    begin

-      last_bit:= 0;

-      while tmp /= 0 loop

-        tmp:= tmp / 2;

-        last_bit:= last_bit + 1;

-      end loop;

-    end;

-    result:= Same_as_Basic_natural(block_value);

-    -- If the following loop was on all blocks,

-    -- the shift by Block_type_bits in the loop could do meaningless

-    -- things the case Basic_int'Size <= Block_Type'Size

-    for b in reverse 0 .. large.last_used-1 loop

-      result:= Shift_Left(result, Block_type_bits);

-      -- An overflow is not detected by shifting (it's the way we want it!)

-      -- so we need to detect the overall overflow by locating the

-      -- last bit.

-      last_bit:= last_bit + Block_type_bits;

-      if last_bit > Basic_int'Size - 1 then

-        -- ^ "- 1" because of sign bit in Basic_int

-        raise Cannot_fit;

-      end if;

-      result:= result + Same_as_Basic_natural(large.blk(b));

-    end loop;

-    if large.neg then

-      return -Basic_int(result);

-    else

-      return Basic_int(result);

-    end if;

-  end Basic;

-

-  -- 14-Feb-2002: "zero" bug fixed by Duncan Sands

-  procedure Fill(what: out Multi_int; with_smaller:Multi_int) is

-    l: constant Index_int:= with_smaller.last_used;

-  begin

-    if Debug then Check_internal.Test(with_smaller); end if;

-    what.zero:= with_smaller.zero;

-

-    if with_smaller.zero then

-      return;

-    end if;

-

-    if what.n < l then

-      raise Array_too_small;   -- contenant trop petit

-    end if;

-

-    what.blk(0..l):= with_smaller.blk(0..l); -- copy contents

-    what.neg:= with_smaller.neg;

-    what.last_used:= l;

-  end Fill;

-

-  procedure Fill(what:out Multi_int; with_basic: Basic_int) is

-  begin

-    Fill( what, Multi(with_basic) );

-  end Fill;

-

-  function Bits_per_block return Positive is

-  begin

-    return Block_type_bits;

-  end Bits_per_block;

-

-  ---------------------------

-  ----- Unary operators -----

-  ---------------------------

-

-  function "+" (i: Multi_int) return Multi_int is begin return i; end "+";

-

-  procedure Opp(i: in out Multi_int) is

-  begin

-    i.neg:= not i.neg; -- -0 possible, anyway i.zero = True in such a case

-  end Opp;

-

-  function "-" (i: Multi_int) return Multi_int is

-    res: Multi_int(i.n):= i; -- copy + stack :-(

-  begin

-    Opp(res);

-    return res;

-  end "-";

-

-  procedure Abso(i: in out Multi_int) is

-  begin

-    i.neg:= False;

-  end Abso;

-

-  function "Abs" (i: Multi_int) return Multi_int is

-    abs_i: Multi_int(i.n):= i; -- copy + stack :-(

-  begin

-    if Debug then Check_internal.Test(i); end if;

-    abs_i.neg:= False;

-    return abs_i;

-  end "Abs";

-

-  function Sign(i: Multi_int) return Basic_int is

-  begin

-    if    i.zero then return  0;

-    elsif i.neg  then return -1;

-    else              return +1;

-    end if;

-  end Sign;

-

-  function Even(i: Multi_int) return Boolean is

-  begin

-    return i.zero or i.blk(0) mod 2 = 0;

-  end Even;

-

-  function Odd (i: Multi_int) return Boolean is

-  begin

-    return (not i.zero) and i.blk(0) mod 2 = 1;

-  end Odd;

-

-  ----------------------------

-  ----- Binary operators -----

-  ----------------------------

-

-  -- Internal algorithm to add two numbers AS POSITIVE ( > 0 ) !

-

-  procedure Add_absolute(i1,i2: in Multi_int; i3: out Multi_int) is

-    l1: constant Index_int:= i1.last_used;

-    l2: constant Index_int:= i2.last_used;

-    min_ind: constant Index_int:= Min( l1, l2 );

-    max_ind: constant Index_int:= Max( l1, l2 );

-    s: Long_Block_type:= 0;

-    retenue_finale: Block_type;

-  begin

-    if Debug then Check_internal.Test(i1); Check_internal.Test(i2); end if;

-

-    if max_ind > i3.n then

-      raise Result_undersized;

-    end if; -- 17-Feb-2002

-

-    -- (1) On additionne sur le <<support commun>>

-    for ind in 0 .. min_ind loop

-      s:= Long_Block_type(i1.blk(ind)) + Long_Block_type(i2.blk(ind)) +

-              Shift_Right(s, Block_type_bits); --  (retenue)

-      i3.blk(ind):= Block_type(s and maxblock);

-      -- NB: dans un cas de Add(a,b,a) ou Add(a,b,b),

-      -- i1.blk(ind) ou i2.blk(ind) est modifie en meme temps!

-    end loop;

-

-    -- (2) On poursuit au besoin si i1 a plus de blocs...

-    if l1 > min_ind then

-      for ind in min_ind+1 .. max_ind loop

-        s:= Long_Block_type(i1.blk(ind)) +

-              Shift_Right(s, Block_type_bits); --  (retenue)

-        i3.blk(ind):= Block_type(s and maxblock);

-      end loop;

-    -- ... ou bien si i2 en a plus.

-    elsif l2 > min_ind then

-      for ind in min_ind+1 .. max_ind loop

-        s:= Long_Block_type(i2.blk(ind)) +

-              Shift_Right(s, Block_type_bits); --  (retenue)

-        i3.blk(ind):= Block_type(s and maxblock);

-      end loop;

-    end if;

-

-    -- (3) Il peut rester une retenue

-    retenue_finale:= Block_type(Shift_Right(s, Block_type_bits));

-    if retenue_finale /= 0 then

-      if max_ind+1 > i3.n then

-        raise Result_undersized;

-      end if; -- 17-Feb-2002

-      i3.blk(max_ind+1):= retenue_finale;

-      i3.last_used:= max_ind+1;

-    else

-      i3.last_used:= max_ind;

-    end if;

-

-    -- (4) i3 = i1+i2 > 0

-    i3.neg:= False;

-    i3.zero:= False;

-

-  end Add_absolute;

-

-  -- Internal algorithm to subtract two numbers AS POSITIVE ( > 0 ) !

-

-  procedure Sub_absolute(i1,i2: in Multi_int; i3: in out Multi_int;

-                         sgn: out Boolean) is

-    l1: constant Index_int:= i1.last_used;

-    l2: constant Index_int:= i2.last_used;

-    max_ind: constant Index_int:= Max( l1, l2 );

-    ai, bi: Long_Block_type;

-    s: Block_type;

-    retenue_finale: Long_Block_type;

-  begin

-    if Debug then Check_internal.Test(i1); Check_internal.Test(i2); end if;

-

-    if max_ind > i3.n then raise Result_undersized; end if; -- 17-Feb-2002

-

-    i3.last_used:= 0;

-    i3.zero:= True;

-    s:= 0;

-

-    -- (1) Soustraction avec retenue

-    for ind in 0 .. max_ind loop

-      if ind <= l1 then

-        ai:= Long_Block_type(i1.blk(ind));

-      else

-        ai:= 0;

-      end if;

-      if ind <= l2 then

-        bi:= Long_Block_type(i2.blk(ind)) + Long_Block_type(s);

-      else

-        bi:= Long_Block_type(s);

-      end if;

-

-      if ai < bi then

-        ai:= ai + cardblock;

-        s:= 1;

-      else

-        s:= 0;

-      end if;

-

-      i3.blk(ind):= Block_type(ai-bi);

-      -- NB: dans un cas de Sub(a,b,a) ou Sub(a,b,b),

-      -- i1.blk(ind) ou i2.blk(ind) est modifie en meme temps!

-

-      if i3.blk(ind) /= 0 then -- au passage, on corrige .last_used et .zero

-        i3.last_used:= ind;

-        i3.zero:= False;

-      end if;

-    end loop;

-

-    -- (2) Traitement de la derni\`ere retenue

-    if s = 0 then

-      i3.neg := False;

-      sgn    := False;

-    else

-      i3.neg := True;

-      sgn    := True;

-      i3.last_used:= 0;

-      retenue_finale:= 1; -- on fait "9-chaque chiffre" et on ajoute 1 au tout

-      for i in 0 .. max_ind loop

-        retenue_finale:=

-          Long_Block_type(maxblock) -

-          Long_Block_type(i3.blk(i)) + retenue_finale;

-        i3.blk(i):= Block_type(retenue_finale and maxblock);

-        if i3.blk(i) /= 0 then

-          i3.last_used:= i;

-        end if;

-        retenue_finale:= Shift_Right(retenue_finale, Block_type_bits);

-      end loop;

-    end if;

-

-  end Sub_absolute;

-

-  procedure Add(i1,i2: in Multi_int; i3: in out Multi_int) is

-    sgn: Boolean;

-  begin

-    -- (1) Les cas o\`u i1 ou i2 = 0

-    if i1.zero and i2.zero then

-      i3.zero:= True;

-    elsif i1.zero then

-      Fill( i3, i2 );

-    elsif i2.zero then

-      Fill( i3, i1 );

-    -- (2) Maintenant: i1 /= 0 et i2 /= 0; on regarde les signes

-    -- (2.1) Facile: i1 et i2 de m\^eme signe

-    elsif i1.neg = i2.neg then

-      Add_absolute( i1,i2, i3 ); -- On fait comme si i1>0 et i2>0

-      i3.neg:= i1.neg;           -- et on met le bon signe

-    -- (2.2) i1 < 0, i2 > 0, donc i3 = i2 - abs(i1)

-    elsif i1.neg and not i2.neg then

-      Sub_absolute( i2,i1, i3, sgn);

-    -- (2.3) i1 > 0, i2 < 0, donc i3 = i1 - abs(i2)

-    elsif i2.neg and not i1.neg then

-      Sub_absolute( i1,i2, i3, sgn );

-    end if;

-  end Add;

-

-  function "+" (i1,i2: Multi_int) return Multi_int is

-    somme: Multi_int( Max(i1.n, i2.n) + 1 );

-  begin

-    Add( i1,i2, somme );

-    return somme;

-  end "+";

-

-  procedure Sub(i1,i2: in Multi_int; i3: in out Multi_int) is

-    sgn: Boolean;

-  begin

-    -- (1) Les cas o\`u i1 ou i2 = 0

-    if    i1.zero and i2.zero then i3.zero:= True;

-    elsif i1.zero then Fill( i3, i2 ); i3.neg:= not i2.neg;

-    elsif i2.zero then Fill( i3, i1 );

-

-    -- (2) Maintenant: i1 /= 0 et i2 /= 0; on regarde les signes

-

-    -- (2.1) Facile: i1 et i2 de m\^eme signe

-    elsif i1.neg = i2.neg then

-      Sub_absolute( i1,i2, i3, sgn ); -- On fait comme si i1>0 et i2>0

-                                      -- et on met le bon signe

-    i3.neg:= i1.neg xor sgn;

-    -- 26-Mar-2002: equivalent a:

-    --      if i1.neg then

-    --        i3.neg:= NOT sgn;

-    --      else

-    --        i3.neg:= sgn;

-    --      end if;

-

-    -- (2.2) i1 < 0, i2 > 0, donc i3 = i1-i2 = - (abs(i1) + abs(i2))

-    elsif i1.neg and not i2.neg then

-      Add_absolute( i1,i2, i3 );

-      i3.neg:= True;

-

-    -- (2.3) i1 > 0, i2 < 0, donc i3 = i1-i2 = i1 + (-i2) = i1 + abs(i2)

-    elsif i2.neg and not i1.neg then

-      Add_absolute( i1,i2, i3 );

-    end if;

-

-  end Sub;

-

-  function "-" (i1,i2: Multi_int) return Multi_int is

-    diff: Multi_int( Max(i1.n, i2.n) + 1); -- +1: retenue possible (add_abs.)

-  begin

-    Sub( i1,i2, diff );

-    return diff;

-  end "-";

-

-  function "+" (i1: Multi_int; i2: Basic_int) return Multi_int is

-  begin return i1 + Multi(i2); end "+";

-

-  function "+" (i1: Basic_int; i2: Multi_int) return Multi_int is

-  begin return Multi(i1) + i2; end "+";

-

-  function "-" (i1: Multi_int; i2: Basic_int) return Multi_int is

-  begin return i1 - Multi(i2); end "-";

-

-  function "-" (i1: Basic_int; i2: Multi_int) return Multi_int is

-  begin return Multi(i1) - i2; end "-";

-

-  ----- Begin of MULTIPLICATION part -----

-

-  -- Added 2006: choice to copy result into i3 or write directly into i3

-  generic

-    copy: Boolean;

-  procedure Multiply_internal_m_m(i1,i2: in Multi_int; i3: in out Multi_int);

-

-  type p_Block_array is access Block_array;

-  procedure Dispose is new Ada.Unchecked_Deallocation(Block_array,p_Block_array);

-

-  -------------------

-  -- Multi * Multi --

-  -------------------

-

-  -- To do: implement a faster algorithm.

-  -- 1) Karatsuba's algorithm

-  -- Ada code for string arithm exists !!

-  -- http://www.csc.liv.ac.uk/~ped/teachadmin/algor/karatsuba.ada

-  -- 2) Better: Schönhage-Strassen algorithm (no Ada code)

-

-  procedure Multiply_internal_m_m(i1,i2: in Multi_int; i3: in out Multi_int) is

-    l1: constant Index_int:= i1.last_used;

-    l2: constant Index_int:= i2.last_used;

-    last_max: constant Index_int:= l1 + l2 + 2;

-    prod,sum_carry,rk,i1j : Long_Block_type;

-    i,k: Index_int;

-    res: p_Block_array;

-    -- res: buffer used in the "copy" variant to avoid

-    -- problems with Multiply(i,j,i) or Multiply(j,i,i)

-  begin

-    if i1.zero or i2.zero then

-      i3.zero:= True;

-      return;

-    end if;

-

-    if last_max > i3.n then

-      raise Result_undersized;

-    end if;

-

-    if copy then

-      res:= new Block_array( 0..last_max );

-      for k in res'Range loop res(k):= 0; end loop;

-      -- Seems slower :-( :  res:= new Block_array'( 0..last_max => 0);

-    else

-      for k in 0..last_max loop i3.blk(k):= 0; end loop;

-      -- Slower :-( :  i3.blk(0..last_max):= (others => 0);

-    end if;

-

-    i3.zero:= False;

-    i3.last_used:= last_max;

-    -- NB: va changer i1.last_used ou i2.last_used si

-    -- i1 ou i2 et i3 sont les memes

-

-    for j in 0..l1 loop

-      i1j:= Long_Block_type(i1.blk(j));

-      sum_carry:= 0;

-      i:= 0;

-      k:= j;

-      loop

-        if i <= l2 then

-          prod:= i1j * Long_Block_type(i2.blk(i));

-        else

-          exit when sum_carry = 0; -- nothing more to add

-          prod:= 0;

-        end if;

-        if copy then

-          rk:= Long_Block_type(res(k));

-        else

-          rk:= Long_Block_type(i3.blk(k));

-        end if;

-        sum_carry:= rk + prod + sum_carry;

-        if copy then

-          res(k):= Block_type(sum_carry and maxblock); -- somme

-        else

-          i3.blk(k):= Block_type(sum_carry and maxblock); -- somme

-        end if;

-        sum_carry:= Shift_Right(sum_carry, Block_type_bits); -- retenue

-        i:= i + 1;

-        k:= k + 1;

-      end loop;

-    end loop;

-

-    if copy then

-      i3.blk(res'Range):= res.all;

-      Dispose(res);

-    end if;

-

-    Reduce_last_nonzero( i3 );

-

-    i3.neg:= i1.neg /= i2.neg;

-

-  end Multiply_internal_m_m;

-

-  procedure Multiply_internal_copy is

-    new Multiply_internal_m_m( copy => True );

-  procedure Multiply_internal_copy_export(i1,i2: in Multi_int; i3: in out Multi_int)

-    renames Multiply_internal_copy;

-  -- ^ At least GNAT <= GPL 2006 requires the trick with renames...

-  --   ObjectAda 7.2.2 too -> there must be a good reason...

-

-  procedure Multiply_internal_no_copy is

-    new Multiply_internal_m_m( copy => False );

-

-  -------------------

-  -- Multi * Basic --

-  -- added 2006    --

-  -------------------

-

-  generic

-    copy: Boolean;

-  procedure Multiply_internal_m_b(i1: in Multi_int; i2: Basic_int; i3: in out Multi_int);

-

-  procedure Multiply_internal_m_b(i1: in Multi_int; i2: Basic_int; i3: in out Multi_int) is

-    l1: constant Index_int:= i1.last_used;

-    last_max: constant Index_int:= l1 + 2;

-    prod,sum_carry,rk,i2a : Long_Block_type;

-    k: Index_int;

-    res: p_Block_array;

-    -- res: buffer used in the "copy" variant to avoid

-    -- problems with Multiply(i,j,i) or Multiply(j,i,i)

-  begin

-    if i1.zero or i2=0 then

-      i3.zero:= True;

-      return;

-    end if;

-

-    if last_max > i3.n then

-      raise Result_undersized;

-    end if;

-

-    if copy then

-      res:= new Block_array( 0..last_max );

-      for k in res'Range loop res(k):= 0; end loop;

-      -- Seems slower :-( :  res:= new Block_array'( 0..last_max => 0);

-    else

-      for k in 0..last_max loop i3.blk(k):= 0; end loop;

-      -- Slower :-( :  i3.blk(0..last_max):= (others => 0);

-    end if;

-

-    i3.zero:= False;

-    i3.last_used:= last_max;

-    -- NB: va changer i1.last_used ou i2.last_used si i1 ou i2 et i3 sont les memes

-    i2a:= Long_Block_type(abs i2);

-

-    for j in 0..l1 loop

-      k:= j;

-      sum_carry:= 0;

-      prod:= Long_Block_type(i1.blk(j)) * i2a;

-      loop

-        if copy then

-          rk:= Long_Block_type(res(k));

-        else

-          rk:= Long_Block_type(i3.blk(k));

-        end if;

-        sum_carry:= rk + prod + sum_carry;

-        if copy then

-          res(k):= Block_type(sum_carry and maxblock); -- somme

-        else

-          i3.blk(k):= Block_type(sum_carry and maxblock); -- somme

-        end if;

-        sum_carry:= Shift_Right(sum_carry, Block_type_bits); -- retenue

-        exit when sum_carry = 0; -- nothing more to add

-        prod:= 0;

-        k:= k + 1;

-      end loop;

-    end loop;

-

-    if copy then

-      i3.blk(res'Range):= res.all;

-      Dispose(res);

-    end if;

-

-    Reduce_last_nonzero( i3 );

-

-    i3.neg:= i1.neg /= (i2 < 0);

-

-  end Multiply_internal_m_b;

-

-  procedure Multiply_internal_copy is

-    new Multiply_internal_m_b( copy => True );

-

-  procedure Multiply_internal_no_copy is

-    new Multiply_internal_m_b( copy => False );

-

-  procedure Multiply(i1,i2: in Multi_int; i3: in out Multi_int) is

-    use System;

-  begin

-    if Debug then

-      declare

-        m1: constant Multi_int:= i1;

-        m2: constant Multi_int:= i2;

-      begin

-        Multiply_internal_no_copy(m1,m2,i3);

-        Check_internal.Check_Multiplication(m1,m2,i3);

-      end;

-    else

-      if i1'Address = i3'Address or i2'Address = i3'Address then

-        -- Ada.Text_IO.Put_Line("* with copy");

-        Multiply_internal_copy(i1,i2,i3);

-      else

-        -- Ada.Text_IO.Put_Line("* without copy");

-        Multiply_internal_no_copy(i1,i2,i3);

-      end if;

-    end if;

-  end Multiply;

-

-  procedure Multiply(i1: in Multi_int; i2: Basic_int; i3: in out Multi_int) is

-    use System;

-  begin

-    if Debug then

-      declare

-        m1: constant Multi_int:= i1;

-        m2: constant Basic_int:= i2;

-      begin

-        Multiply_internal_no_copy(m1,m2,i3);

-        Check_internal.Check_Multiplication(m1,Multi(m2),i3);

-      end;

-    else

-      if i1'Address = i3'Address or i2'Address = i3'Address then

-        -- Ada.Text_IO.Put_Line("* with copy");

-        Multiply_internal_copy(i1,i2,i3);

-      else

-        -- Ada.Text_IO.Put_Line("* without copy");

-        Multiply_internal_no_copy(i1,i2,i3);

-      end if;

-    end if;

-  end Multiply;

-

-  function "*" (i1,i2: Multi_int) return Multi_int is

-  begin

-    if i1.zero or i2.zero then

-      return zero;

-    else

-      declare

-        prod: Multi_int( i1.last_used + i2.last_used + 2 );

-      begin

-        Multiply( i1,i2, prod );

-        return prod;

-      end;

-    end if;

-  end "*";

-

-  function "*" (i1: Multi_int; i2: Basic_int) return Multi_int is

-  begin

-    if i1.zero or i2=0 then

-      return zero;

-    else

-      declare

-        prod: Multi_int( i1.last_used + 4 );

-      begin

-        Multiply( i1,i2, prod );

-        return prod;

-      end;

-    end if;

-  end "*";

-

-  function "*" (i1: Basic_int; i2: Multi_int) return Multi_int is

-  begin

-    if i2.zero or i1=0 then

-      return zero;

-    else

-      declare

-        prod: Multi_int( i2.last_used + 4 );

-      begin

-        Multiply( i2,i1, prod );

-        return prod;

-      end;

-    end if;

-  end "*";

-

-  ----- Begin of DIVISION part -----

-

-  -- Interne: Division et reste en 1 coup

-

-  procedure Div_Rem(a,b: Long_Block_type; q,r: out Long_Block_type) is

-    Conflict_with_REM: exception;

-  begin

-    q:= a / b;

-    r:= a - b*q;

-    if Debug and then r /= (a rem b) then

-      raise Conflict_with_REM;

-    end if;

-  end Div_Rem;

-

-  procedure Divide_absolute_normalized ( u: in out Multi_int; -- output: u = r

-                                         v: in     Multi_int;

-                                         q: in out Multi_int  ) is

-    qi: Index_int:= u.last_used - v.last_used - 1; -- was: q.n; D.S. Feb-2002

-    v1: constant Long_Block_type:= Long_Block_type(v.blk(v.last_used  ));

-    v2: constant Long_Block_type:= Long_Block_type(v.blk(v.last_used-1));

-

-    vlast     : constant Index_int:= v.last_used;

-    v1L       : constant Long_Block_type := v1;

-    guess,

-    comparand : Long_Block_type ;

-

-    function Divide_subtract ( ustart: Index_int ) return Block_type is

-      ui    : Index_int;

-      carry : Long_Block_type;

-    begin

-      if guess = 0 then

-        return 0;

-      end if;

-      ui:= ustart;

-      carry:= 0;

-

-      -- On soustrait (le chiffre du quotient) * diviseur au dividende

-

-      for vi in 0 .. vlast loop

-        declare

-          prod: constant Long_Block_type   := Long_Block_type(v.blk(vi)) * guess + carry;

-          bpro: constant Block_type:= Block_type(prod and maxblock);

-          diff: constant Long_Block_type_signed   := Long_Block_type_signed(u.blk(ui)) - Long_Block_type_signed(bpro);

-        begin

-          if diff < 0 then

-            u.blk(ui) := Block_type(diff + cardblock);

-            carry := Shift_Right(prod, Block_type_bits) + 1;

-          else

-            u.blk(ui) := Block_type(diff);

-            carry := Shift_Right(prod, Block_type_bits);

-          end if;

-          ui:= ui + 1;

-        end;

-      end loop;

-

-      if carry = 0 then

-        return Block_type(guess and maxblock);

-      end if;

-

-      declare

-        diff: constant Long_Block_type_signed :=

-          Long_Block_type_signed(u.blk(ui)) - Long_Block_type_signed(carry and maxblock);

-      begin

-          if diff < 0 then

-            u.blk(ui) := Block_type(diff + cardblock); -- carry generated

-          else

-            u.blk(ui) := Block_type(diff);

-            return Block_type(guess and maxblock);

-          end if;

-      end;

-

-      -- Carry was generated

-      declare

-        icarry: Block_type := 0;

-      begin

-        ui := ustart;

-        for vi in 0 .. vlast loop

-          declare

-            sum: constant Long_Block_type :=

-              Long_Block_type(v.blk(vi)) +

-              Long_Block_type(u.blk(ui)) +

-              Long_Block_type(icarry);

-          begin

-            u.blk(ui) := Block_type(sum and maxblock);

-            ui:= ui + 1;

-            icarry := Block_type(Shift_Right(sum, Block_type_bits));

-          end;

-        end loop;

-

-        if icarry = 1 then

-          u.blk(ui) := Block_type((Long_Block_type(u.blk(ui))+1) and maxblock);

-        end if;

-      end;

-

-      return Block_type((guess-1) and maxblock);

-

-    end Divide_subtract;

-

-    is_q_zero: Boolean:= True;

-

-  begin -- Divide_absolute_normalized

-    -- for i in q.blk'Range loop q.blk(i):= 0; end loop;

-    --

-    -- ^ zeroing useless: q.last_used = u.last_used-v.last_used-1

-    --   and q.blk(0..q.last_used) is written below q.blk(qi) := ...

-    --   GM 4-nov-2006

-

-    q.last_used:= qi; -- was: q.n; D.S. Feb-2002

-

-    for j in reverse vlast+1 .. u.last_used loop

-      declare

-        uj : constant Long_Block_type := Long_Block_type(u.blk(j));

-        uj1: constant Long_Block_type := Long_Block_type(u.blk(j-1));

-        uj2: constant Long_Block_type := Long_Block_type(u.blk(j-2));

-        ujL: Long_Block_type;

-        rmL: Long_Block_type;

-      begin

-        ujL := Shift_Left(uj, Block_type_bits) + uj1;

-        Div_Rem( ujL, v1L, guess, rmL );

-        comparand := Shift_Left(rmL, Block_type_bits) + uj2;

-

-        while comparand < v2 * guess loop

-          guess:= guess - 1;

-          comparand:= comparand + Shift_Left(v1L, Block_type_bits);

-          exit when comparand > cardblock * cardblock;

-        end loop;

-

-        q.blk(qi) := Divide_subtract( j - vlast - 1 );

-

-        if q.blk(qi) /= 0 and then is_q_zero then -- n'arrive que 0 ou 1 fois

-          is_q_zero:= False;

-          q.last_used:= qi;

-        end if;

-

-        qi:= qi - 1;

-      end;

-

-    end loop; -- j

-

-    q.zero:= is_q_zero;

-

-    if Debug then Check_internal.Test(q); end if;

-

-  end Divide_absolute_normalized;

-

-  procedure Divide_absolute_big_small ( u:   in     Multi_int;

-                                        v:   in     Long_Block_type;

-                                        q:      out Multi_int;

-                                        r:      out Long_Block_type ) is

-    n: Long_Block_type;

-    Quotient_constraint_error: exception;

-    last_u_nz:  constant Index_int:= u.last_used;

-    u_zero: constant Boolean:= u.zero;

-    -- in case u and q are the same variables

-    is_q_zero: Boolean:= True;

-  begin

-    if q.n < last_u_nz then raise Quotient_constraint_error; end if;

-    q.last_used:= 0;

-    q.neg:= False;

-    r:= 0;

-    if not u_zero then

-      for i in reverse 0 .. last_u_nz loop

-        n:= Long_Block_type(u.blk(i)) + Shift_Left(r, Block_type_bits);

-        r:= n mod v;

-        q.blk(i):= Block_type(n  /  v);

-        if q.blk(i)/= 0 and then is_q_zero then

-          is_q_zero:= False;

-          q.last_used:= i;

-        end if;

-      end loop;

-      q.zero:= is_q_zero;

-    end if;

-  end Divide_absolute_big_small;

-

-  procedure Solve_signs_for_Div_Rem (i1n,i2n: in Boolean; qn,rn: out Boolean) is

-  begin

-    -- Invariant: i1= i2*q+r   on cherche (pos) = (pos)*(pos)+(pos)

-

-    if i1n and i2n then        -- i1<0;  i2<0  (-i1) = (-i2) *  q  + (-r)

-      qn:= False; -- Quotient > 0

-    --      rn:= True;  -- Reste    < 0

-    elsif i1n then             -- i1<0;  i2>0  (-i1) =   i2  *(-q) + (-r)

-      qn:= True;  -- Quotient < 0

-    --      rn:= True;  -- Reste    < 0

-    elsif i2n then             -- i1>0;  i2<0    i1  = (-i2) *(-q) +   r

-      qn:= True;  -- Quotient < 0

-    --      rn:= False; -- Reste    > 0

-    else                       -- i1>0;  i2>0    i1  =   i2  *  q  +   r

-      qn:= False; -- Quotient > 0

-    --      rn:= False; -- Reste    > 0

-    end if;

-    -- on observe que... "(A rem B) has the sign of A " ARM 4.5.5

-    -- en effet on peut mettre:

-    rn:= i1n;

-  end Solve_signs_for_Div_Rem;

-

-  procedure Div_Rem (i1: in     Multi_int; i2: in     Basic_int;

-                     q :    out Multi_int;  r:    out Basic_int) is

-    i1_neg: constant Boolean:= i1.neg;

-    -- in case i1 and q are the same variables

-    rneg: Boolean;

-    lai2, lr: Long_Block_type;

-  begin

-    if Debug then Check_internal.Test(i1); end if;

-    if i2=0 then raise Division_by_zero; end if;

-

-    if i1.zero then -- 15-Feb-2002: 0/i2

-      q.zero:= True;

-      r:= 0;

-      return;

-    end if;

-

-    lai2:= Long_Block_type(abs i2);

-    Divide_absolute_big_small( i1, lai2, q, lr );

-    r:= Basic_int(lr);

-

-    Solve_signs_for_Div_Rem( i1_neg,i2<0, q.neg, rneg );

-    if rneg then r:= -r; end if;

-

-  end Div_Rem;

-

-  type Div_Rem_mode is (div_only, both);

-

-  generic

-    div_rem_output: Div_Rem_mode;

-  procedure Div_Rem_internal (i1,i2: in Multi_int; q,r: in out Multi_int);

-

-  procedure Div_Rem_internal (i1,i2: in Multi_int; q,r: in out Multi_int) is

-

-    -- Calculate u/v

-

-    procedure Divide_absolute ( u,v: in     Multi_int;

-                                q,r: in out Multi_int  ) is

-      shift: Integer:= 0;

-      v1: Block_type:= v.blk(v.last_used);

-      v_zero, v1_zero: exception;

-      u_work: Multi_int(u.last_used+2);

-      use System;

-

-      procedure Normalization ( source: in     Multi_int;

-                                target: in out Multi_int ) is

-        carry: Block_type:= 0;

-        tl: constant Index_int:= target.last_used;

-        blk:  Block_type;

-      begin

-        for i in 0 .. source.last_used loop

-          blk:= source.blk(i);

-          target.blk(i) := Shift_Left(blk, shift) + carry;

-          carry         := Shift_Right(blk, Block_type_bits - shift);

-        end loop;

-        if source.last_used < tl then

-          target.blk(source.last_used+1):= carry;

-        end if;

-        for i in source.last_used+2 .. tl  loop

-          target.blk(i):= 0;

-        end loop;

-      end Normalization;

-

-      procedure Unnormalization ( m: in out Multi_int) is

-        carry: Block_type:= 0;

-        blk:  Block_type;

-      begin

-        for i in reverse 0 .. m.last_used loop

-          blk:= m.blk(i);

-          m.blk(i) := Shift_Right(blk, shift) + carry;

-          carry    := Shift_Left(blk, Block_type_bits - shift);

-        end loop;

-      end Unnormalization;

-

-    begin -- Divide_absolute (multi u / multi v)

-

-      if Debug then

-        if v.zero then raise v_zero; end if;

-        if v1=0 then raise v1_zero; end if;

-      end if;

-

-      -- Calculate shift needed to normalize

-      u_work.last_used:= u_work.n;

-      u_work.zero:= False;

-      while v1 < 2**(Block_type_bits-1) loop

-        shift:= shift + 1;

-        v1:= v1 * 2;

-      end loop;

-      if shift = 0 then                  -- no shift needed

-        u_work.blk( 0..u.last_used ):= u.blk( 0..u.last_used );

-        u_work.blk( u.last_used+1 .. u_work.last_used):= (0,0);

-        -- Now, u is copied, so a Div_Rem(u, v, u, r) won't crash

-

-        if v'Address = q'Address then

-          declare

-            v_work: Multi_int(v.last_used);

-          begin

-            -- 23-Feb-2002: also copy v, in case of a Div_Rem(u, v, v, r)

-            v_work.blk( 0..v.last_used ):= v.blk( 0..v.last_used );

-            v_work.neg      := v.neg;

-            v_work.zero     := v.zero;

-            v_work.last_used:= v.last_used;

-            -- Now, u is copied, so a Div_Rem(u, v, v, r) won't crash

-            -- Ada.Text_IO.Put_Line("* divisor with copy");

-            Divide_absolute_normalized( u_work,v_work, q );

-          end;

-        else

-          -- Ada.Text_IO.Put_Line("* divisor without copy");

-          Divide_absolute_normalized( u_work,v, q );

-        end if;

-

-      else  -- shift needed

-        declare

-          v_work: Multi_int(v.last_used);

-        begin

-          v_work.last_used:= v_work.n;

-          Normalization( u, u_work );

-          Normalization( v, v_work );

-          Reduce_last_nonzero( v_work );

-

-          Divide_absolute_normalized( u_work,v_work, q );

-        end;

-

-        if div_rem_output /= div_only then

-          Unnormalization( u_work );

-        end if;

-      end if;

-      q.neg:= False; -- check friendly

-      if div_rem_output /= div_only then

-        u_work.neg:= False; -- check friendly

-        Reduce_last_nonzero( u_work );

-        Fill( r, u_work );

-      end if;

-

-    end Divide_absolute;

-

-    l1: constant Index_int:= i1.last_used;

-    l2: constant Index_int:= i2.last_used;

-    rl: Long_Block_type;

-  begin -- Div_Rem_internal

-    if i2.zero then raise Division_by_zero; end if;

-

-    if i1.zero then -- 15-Feb-2002: 0/i2

-      q.zero:= True;

-      r.zero:= True;

-      return;

-    end if;

-

-    if q.n < l1 - l2 then

-      -- 17-Feb-2002

-      raise Quotient_undersized;

-    end if;

-

-    if div_rem_output /= div_only and then r.n < Max( l1, l2 ) then

-      -- 17-Feb-2002

-      raise Remainder_undersized;

-    end if;

-

-    if l2 = 0 then

-      if l1 = 0 then      -- On a affaire a une ridicule division d'entiers

-        q.blk(0):= i1.blk(0) / i2.blk(0);

-        if div_rem_output /= div_only then

-          r.blk(0):= Block_type(

-            abs(

-                 Long_Block_type_signed(i1.blk(0))

-               - Long_Block_type_signed(i2.blk(0))

-               * Long_Block_type_signed(q.blk(0))

-            )

-          );

-        end if;

-        q.zero:= q.blk(0) = 0;

-        q.last_used:= 0;

-      else                -- multi / entier

-        Divide_absolute_big_small( i1, Long_Block_type(i2.blk(0)), q, rl );

-        if div_rem_output /= div_only then

-          r.blk(0):= Block_type(rl);

-        end if;

-      end if;

-      if div_rem_output /= div_only then

-        r.zero:= r.blk(0) = 0;

-        r.last_used:= 0;

-      end if;

-

-    else  -- multi / multi

-

-      case Compare_absolute(i2 , i1) is

-

-        when greater =>

-          q.zero:= True;    -- q:=  0;

-          q.last_used:= 0;

-          q.neg:= False;

-

-          if div_rem_output /= div_only then

-            Fill( r, i1 );  -- r:= i1, q:=0 car i1 = 0 * i2 (>i1 en v.abs) + r

-          end if;

-          return;

-

-        when equal =>

-          Fill( q, one ); -- Fill( q, Multi(1) );

-          r.zero:= True;  -- Fill( r, Multi(0) );

-

-        when smaller => -- cas <<normal>>: diviseur < dividende

-

-          Divide_absolute( i1,i2, q,r );

-

-      end case;

-    end if;

-

-    Solve_signs_for_Div_Rem( i1.neg,i2.neg, q.neg,r.neg );

-  end Div_Rem_internal;

-

-  procedure Div_Rem_internal_div_only is

-    new Div_Rem_internal( div_rem_output => div_only );

-

-  procedure Div_Rem_internal_both is

-    new Div_Rem_internal( div_rem_output => both );

-

-  procedure Div_Rem_internal_both_export(i1,i2: in Multi_int; q,r: in out Multi_int)

-    renames Div_Rem_internal_both;

-

-  procedure Div_Rem (i1,i2: in Multi_int; q,r: out Multi_int) is

-  begin

-    if Debug then

-      declare

-        m1: constant Multi_int:= i1;

-        m2: constant Multi_int:= i2;

-      begin

-        Div_Rem_internal_both(m1,m2,q,r);

-        Check_internal.Check_Div_Rem(m1,m2,q,r);

-      end;

-    else

-      Div_Rem_internal_both(i1,i2,q,r);

-    end if;

-  end Div_Rem;

-

-  procedure Divide (i1,i2: in Multi_int; q: out Multi_int) is

-  begin

-    if Debug then

-      declare

-        m1: constant Multi_int:= i1;

-        m2: constant Multi_int:= i2;

-        r: Multi_int( Max( i1.last_used, i2.last_used) + 2 );

-      begin

-        Div_Rem_internal_both(m1,m2,q,r);

-        Check_internal.Check_Div_Rem(m1,m2,q,r);

-      end;

-    else

-      declare

-        r: Multi_int(0); -- Fake

-      begin

-        Div_Rem_internal_div_only(i1,i2,q,r);

-      end;

-    end if;

-  end Divide;

-

-  function "/" (i1,i2: Multi_int) return Multi_int is

-    q: Multi_int( Max( 0, i1.last_used - i2.last_used + 1) );

-    r: Multi_int( Max( i1.last_used, i2.last_used) + 2 );

-  begin

-    Div_Rem(i1,i2, q,r);

-    return q;

-  end "/";

-

-  function "/" (i1: Multi_int; i2: Basic_int) return Multi_int is

-    q: Multi_int(i1.last_used + 1);

-    r: Basic_int;

-  begin

-    Div_Rem(i1,i2, q,r);

-    return q;

-  end "/";

-

-  function "rem" (i1,i2: Multi_int) return Multi_int is

-    q: Multi_int(Max(0,i1.last_used - i2.last_used + 1));

-    r: Multi_int(Max(i1.last_used,i2.last_used) + 2);

-  begin

-    Div_Rem(i1,i2, q,r);

-    return r;

-  end "rem";

-

-  function "rem" (i1: Multi_int; i2: Basic_int) return Multi_int is

-  begin return i1 rem Multi(i2); end "rem";

-

-  function "rem" (i1: Multi_int; i2: Basic_int) return Basic_int is

-    q: Multi_int(i1.last_used + 1);

-    r: Basic_int;

-  begin

-    Div_Rem(i1,i2, q,r);

-    return r;

-  end "rem";

-

-  function "mod" (i1,i2: Multi_int) return Multi_int is

-    q: Multi_int(Max(0,i1.last_used - i2.last_used + 1));

-    r: Multi_int(Max(i1.last_used,i2.last_used) + 2);

-  begin

-    -- Ada RM, 4.5.5 Multiplying Operators

-    -- (8)

-    -- The signed integer modulus operator is defined such that

-    -- the result of A mod B has the sign of B and an absolute value

-    -- less than the absolute value of B; in addition, for some signed

-    -- integer value N, this result satisfies the relation:

-    -- (9) A = B*N + (A mod B)

-

-    Div_Rem(i1,i2, q,r);

-    if r.zero or else i2.neg = r.neg then  --  (A rem B) est nul ou

-      return r;     -- a le meme signe que B, donc (A mod B) = (A rem B)

-    else  -- signe opposes

-      return i2+r;  -- alors (B + (A rem B)) est le bon candidat

-    end if;

-  end "mod";

-

-  function "mod" (i1: Multi_int; i2: Basic_int) return Multi_int is

-  begin return i1 mod Multi(i2); end "mod";

-

-  function "mod" (i1: Multi_int; i2: Basic_int) return Basic_int is

-    r: constant Basic_int:= i1 rem i2;

-  begin

-    if r=0 or else (i2<0) = (r<0) then  --  (A rem B) est nul ou

-      return r;     -- a le meme signe que B, donc (A mod B) = (A rem B)

-    else  -- signe opposes

-      return i2+r;  -- alors (B + (A rem B)) est le bon candidat

-    end if;

-  end "mod";

-

------ End of DIVISION part ------

-

------ Begin of POWER part -------

-

-  procedure Power (i: Multi_int; n: Natural; ipn: out Multi_int) is

-    max_ipn_last: Index_int; -- 17-Feb-2002

-  begin

-    if i.zero then

-      if n=0 then

-        raise Zero_power_zero;

-      else

-        -- The 0**n = 0 case (17-Feb-2002).

-        ipn.zero:= True; -- 4-Nov-2006, was: Fill( ipn, Multi(0) );

-        return;

-      end if;

-    end if;

-

-    max_ipn_last:= ((1+i.last_used) * Index_int(n)-1)+2;

-    if ipn.n < max_ipn_last then

-      raise Result_undersized;

-    end if;

-

-    case n is

-      when 0 => Fill( ipn, one ); -- the i**0 = 1 case

-      when 1 => Fill( ipn, i);    -- the i**1 = i case

-      when others =>

-        declare

-          nn: Natural:= n-1;

-          i0, ii: Multi_int( max_ipn_last );

-        begin

-          Fill(i0, i);

-          Fill(ii, i0 );

-

-          while nn > 0 loop

-            if nn mod 2 = 0 then -- x^(2 c) = (x^2) ^c

-              Mult(i0,i0, i0);

-              nn:= nn / 2;

-            else

-              Mult(i0,ii, ii);

-              nn:= nn - 1;

-            end if;

-          end loop;

-          Fill( ipn, ii);

-        end;

-    end case;

-  end Power;

-

-  function "**" (i: Multi_int; n: Natural) return Multi_int is

-    ipn: Multi_int( (1+i.last_used) * Index_int(n)+2 );

-  begin

-    Power(i,n,ipn);

-    return ipn;

-  end "**";

-

-  procedure Power (i: Multi_int; n: Multi_int; ipn: out Multi_int;

-                   modulo: Multi_int) is

-    max_ipn_last: Index_int;

-  begin

-    if i.zero then

-      if n.zero then

-        raise Zero_power_zero;

-      else

-        -- The 0**n = 0 case (17-Feb-2002).

-        ipn.zero:= True; -- 4-Nov-2006, was: Fill( ipn, Multi(0) );

-        return;

-      end if;

-    end if;

-

-    if n.neg then

-     raise Power_negative;

-    end if;

-

-    if modulo.zero or else (i.neg or modulo.neg) then

-      raise Power_modulo_non_positive;

-    end if;

-

-    max_ipn_last:= 2*modulo.last_used+2;

-    if ipn.n < max_ipn_last then

-      raise Result_undersized;

-    end if;

-

-    if n.zero then

-      Fill( ipn, one ); -- the i**0 = 1 case

-    elsif Equal( n, one ) then

-      Fill( ipn, i);    -- the i**1 = i case

-    else

-      declare

-        nn: Multi_int(n.n):= n;

-        i0, ii, dummy: Multi_int( max_ipn_last );

-        dummy_b: Basic_int;

-      begin

-        Subtract( nn, one, nn ); -- nn:= nn - 1;

-        Fill(i0, i);

-        Fill(ii, i0 );

-

-        while nn > 0 loop

-          if Even(nn) then                -- x^(2 c) = (x^2) ^c

-            Mult(i0,i0, i0);

-            Div_Rem(nn, 2, nn, dummy_b);  -- nn:= nn/2

-            Div_Rem(i0,modulo,dummy,i0);  -- i0:= i0 mod modulo

-          else

-            Mult(i0,ii, ii);

-            Subtract( nn, one, nn );      -- nn:= nn - 1;

-            Div_Rem(ii,modulo,dummy,ii);  -- ii:= ii mod modulo

-          end if;

-        end loop;

-        Fill( ipn, ii);

-      end;

-    end if;

-  end Power;

-

------ End of POWER part ---------

-

------ Comparisons

-

-  function Equal (i1,i2: Multi_int) return Boolean is

-  begin

-    if i1.zero and then i2.zero then

-      return True;

-    end if;

-

-    if i1.zero = i2.zero and then

-       i1.neg  = i2.neg  and then

-       i1.last_used = i2.last_used then

-      return i1.blk(0..i1.last_used) = i2.blk(0..i2.last_used);

-    else

-      return False;

-    end if;

-  end Equal;

-

-  function Equal (i1: Multi_int; i2:Basic_int) return Boolean is

-  begin

-    return Equal( i1, Multi(i2) );

-  end Equal;

-

-  function ">" (i1,i2: Multi_int) return Boolean is

-  begin

-    -- (1) Cas \'evident o\`u:         i1 <= i2

-    if (i1.zero or i1.neg) and then             -- i1 <= 0 et

-       (i2.zero or not i2.neg) then             -- i2 >= 0

-        return False;

-    end if;

-

-    -- (2.1) Cas \'evident o\`u:       i1 > i2

-    if ((not i1.zero) and not i1.neg) and then  -- i1 > 0 et

-       (i2.zero or i2.neg) then                 -- i2 <= 0

-        return True;

-    end if;

-

-    -- (2.2) Cas \'evident o\`u:       i1 > i2

-    if (i1.zero or not i1.neg) and then         -- i1 >= 0 et

-       ((not i2.zero) and i2.neg) then          -- i2 < 0

-        return True;

-    end if;

-

-    -- Cas faciles resolus:

-    -- i1 > i2  -  0  +

-    -------------------

-    --  -       #  F  F

-    --  0       T  F  F

-    --  +       T  T  #

-

-    -- On a les cas avec "#", o\`u i1 et i2 ont le meme signe

-

-    if i1.neg then

-      return not (Compare_absolute (i1,i2) = greater);

-    else

-      return     (Compare_absolute (i1,i2) = greater);

-    end if;

-

-  end ">";

-

-  function ">" (i1: Multi_int; i2:Basic_int) return Boolean is

-  begin

-    return i1 > Multi(i2);

-  end ">";

-

-  function "<" (i1,i2: Multi_int) return Boolean is

-  begin return i2>i1; end "<";

-

-  function "<" (i1: Multi_int; i2:Basic_int) return Boolean is

-  begin

-    return i1 < Multi(i2);

-  end "<";

-

-  function ">=" (i1,i2: Multi_int) return Boolean is

-  begin return not (i2>i1); end ">=";

-

-  function ">=" (i1: Multi_int; i2:Basic_int) return Boolean is

-  begin

-    return i1 >= Multi(i2);

-  end ">=";

-

-  function "<=" (i1,i2: Multi_int) return Boolean is

-  begin return not (i1>i2); end "<=";

-

-  function "<=" (i1: Multi_int; i2:Basic_int) return Boolean is

-  begin

-    return i1 <= Multi(i2);

-  end "<=";

-

-end Multi_precision_integers;