Removed MIT licensed Mathpaqs packages in favour of bindings to GMP

13 files changed, 382 insertions, 2513 deletions

diff --git a/cc0.txt b/license.txt
index 0e259d4..0e259d4 100644
--- a/cc0.txt
+++ b/license.txt
diff --git a/mit.txt b/mit.txt
deleted file mode 100644
index 76dce87..0000000
--- a/mit.txt
+++ /dev/null
@@ -1,19 +0,0 @@
-Copyright (c) <year> <copyright holders>
-
-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.
\ No newline at end of file
diff --git a/notes.txt b/notes.txt
index 19d02a7..5db08bc 100644
--- a/notes.txt
+++ b/notes.txt
@@ -10,8 +10,6 @@ Preference deduplication isn't worth it unless there exists some linear dedupe a
 future direction
 ----------------
 
-Replace Multi_Precision_Integers with own code, with consistent style and licensing
-
 Util to list paper ids that fit specific criteria to doublecheck potential errors.
 
 Add proper tiebreaker handling.
diff --git a/readme.txt b/readme.txt
index 601ddc1..74ddcbc 100644
--- a/readme.txt
+++ b/readme.txt
@@ -9,6 +9,7 @@ To compile this program, the following dependencies are needed:
 
     gnat (of course)
     gnu make
+    gmp
 
 Note however that make isn't strictly necessary if you have a look in the
 makefile to see the compilation commands required.
@@ -42,12 +43,3 @@ elect candidates. On the other hand, the AEC program is also a lot more
 verbose on the distribution of preferences, and doesn't do bulk exclusions.
 
 
-
-Support for bignums (the Multi_Precision_Integers package) was obtained from
-
-    https://sourceforge.net/projects/mathpaqs/
-
-Those source files are licensed under the MIT license. All other source is licensed
-under CC0. See cc0.txt and mit.txt for license details.
-
-
diff --git a/src/bundles.adb b/src/bundles.adb
index af11f37..a6e1472 100644
--- a/src/bundles.adb
+++ b/src/bundles.adb
@@ -56,7 +56,8 @@ package body Bundles is
            (This : in Bundle)
         return Natural is
     begin
-        return Rationals.Floor (Count_Papers (This) * This.Worth);
+        --  The number of votes under all normal circumstances should never exceed MAXINT
+        return Natural (Rationals.Floor (Count_Papers (This) * This.Worth));
     end Count_Votes;
 
 
@@ -78,7 +79,8 @@ package body Bundles is
             Papers :    out Natural) is
     begin
         Papers := Count_Papers (This);
-        Votes := Rationals.Floor (Papers * This.Worth);
+        --  The number of votes under all normal circumstances should never exceed MAXINT
+        Votes := Natural (Rationals.Floor (Papers * This.Worth));
     end Count_Both;
 
 
diff --git a/src/multi_precision_integers-check.adb b/src/multi_precision_integers-check.adb
deleted file mode 100644
index 237be3d..0000000
--- a/src/multi_precision_integers-check.adb
+++ /dev/null
@@ -1,170 +0,0 @@
-

-

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

---

---  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 Ada.Text_IO; use Ada.Text_IO;

--- !! will disappear in favour of Ada.Exceptions

-

-with Multi_precision_integers.IO;

-

-package body Multi_precision_integers.Check is

-

-  package IOi renames Multi_precision_integers.IO;

-

-  -- Don't be afraid by the bug reporting tools,

-  -- they are unlikely to pop out of a programme ;-)

-

-  Bug_file: Ada.Text_IO.File_type;

-  Bug_file_name: constant String:= "mupreint.bug";

-

-  -- The flaw detection (DEBUG mode) could suffer from

-  -- a chicken-and-egg problem: e.g. "*" works, but the "/"

-  -- to verify it doesn't !

-

-  Multiply_flawed_div1_has_rest       : exception;

-  Multiply_flawed_div1_wrong_quotient : exception;

-  Multiply_flawed_div2_has_rest       : exception;

-  Multiply_flawed_div2_wrong_quotient : exception;

-

-  Div_Rem_flawed : exception;

-

-  procedure Open_Bug_Report is

-  begin

-    Create( Bug_file, out_file, Bug_file_name );

-    Put_Line( Bug_file, "Bug in Multi_precision_integers");

-  end Open_Bug_Report;

-

-  procedure Close_Bug_Report is

-  begin

-    Close( Bug_file );

-    -- These console messages can provoke a silent quit in GUI apps,

-    -- so we put them after closing the report.

-    Put_Line( "Bug in Multi_precision_integers !");

-    Put_Line( "For details, read file: " & Bug_file_name);

-  end Close_Bug_Report;

-

-  procedure Test( m: multi_int; test_last: Boolean:= True ) is

-    last_nz: index_int:= 0;

-    Negative_block,

-    Last_index_has_zero,

-    Field_last_outside_range, Field_last_is_negative: exception;

-  begin

-    if m.zero then return; end if; -- 0, nothing to test

-    if m.last_used > m.n then raise Field_last_outside_range; end if;

-    if m.last_used <   0 then raise Field_last_is_negative; end if;

-    for i in 0 .. m.last_used loop

-      if m.blk(i) < 0 then

-        raise Negative_block;

-      end if;

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

-        last_nz:= i;

-      end if;

-    end loop;

-    if test_last and then 0 < last_nz and then last_nz < m.last_used then

-      raise Last_index_has_zero;

-    end if;

-  end Test;

-

-  procedure Check_Multiplication(i1,i2,i3: in multi_int) is

-    jeu: constant:= 5; -- 0 suffit

-    q1: Multi_int( i2.last_used + jeu );

-    r1: Multi_int( i1.last_used + i2.last_used + jeu );

-    q2: Multi_int( jeu );

-    r2: Multi_int( i2.last_used + jeu );

-

-    procedure Bug_Report is

-    begin

-      Open_Bug_Report;

-      Put_Line( Bug_file, "Multiply_and_verify");

-      Put( Bug_file, "i1 ="); IOi.Put_in_blocks(Bug_file, i1); New_Line(Bug_file);

-      Put( Bug_file, "i2 ="); IOi.Put_in_blocks(Bug_file, i2); New_Line(Bug_file);

-      Put( Bug_file, "i3 ="); IOi.Put_in_blocks(Bug_file, i3); New_Line(Bug_file);

-    end Bug_Report;

-

-  begin

-    Test(i1);

-    Test(i2);

-

-    if not (i1.zero or i2.zero) then

-      -- Now we divide i3 by i1, q1 should be = i2

-      Div_Rem_internal_both_export(i3,i1, q1,r1);

-      if not r1.zero then

-        Bug_Report;

-        Close_Bug_Report;

-        raise Multiply_flawed_div1_has_rest;

-      end if;

-      if not Equal( q1, i2 ) then

-        Bug_Report;

-        Put( Bug_file, "q1 ="); IOi.Put_in_blocks(Bug_file, q1); New_Line(Bug_file);

-        Close_Bug_Report;

-        raise Multiply_flawed_div1_wrong_quotient;

-      end if;

-      -- Now we divide q1 by i2, should be = 1

-      Div_Rem_internal_both_export(q1,i2, q2,r2);

-      if not r2.zero then

-        Bug_Report;

-        Close_Bug_Report;

-        raise Multiply_flawed_div2_has_rest;

-      end if;

-      if not Equal( q2, Multi(1) ) then

-        Bug_Report;

-        Put( Bug_file, "q2 ="); IOi.Put_in_blocks(Bug_file, q1); New_Line(Bug_file);

-        Close_Bug_Report;

-        raise Multiply_flawed_div2_wrong_quotient;

-      end if;

-    end if;

-

-  end Check_Multiplication;

-

-  procedure Check_Div_Rem(i1,i2,q,r: in multi_int) is

-

-    procedure Bug_Report is

-    begin

-      Open_Bug_Report;

-      Put_Line( Bug_file, "Div_Rem_and_verify");

-      Put( Bug_file, "i1 ="); IOi.Put_in_blocks(Bug_file, i1); New_Line(Bug_file);

-      Put( Bug_file, "i2 ="); IOi.Put_in_blocks(Bug_file, i2); New_Line(Bug_file);

-      Put( Bug_file, "q  ="); IOi.Put_in_blocks(Bug_file, q); New_Line(Bug_file);

-      Put( Bug_file, "r  ="); IOi.Put_in_blocks(Bug_file, r); New_Line(Bug_file);

-    end Bug_Report;

-

-  begin

-    Test(i1);

-    Test(i2);

-

-    if not Equal( i1, i2*q + r ) then

-      Bug_Report;

-      Close_Bug_Report;

-      raise Div_Rem_flawed;

-    end if;

-

-    Test(q);

-    Test(r);

-  end Check_Div_Rem;

-

-end Multi_precision_integers.Check;

diff --git a/src/multi_precision_integers-check.ads b/src/multi_precision_integers-check.ads
deleted file mode 100644
index 950a21f..0000000
--- a/src/multi_precision_integers-check.ads
+++ /dev/null
@@ -1,41 +0,0 @@
-

-

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

---

---  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.

---

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

-

-

-

-

-package Multi_precision_integers.Check is

-

-  -- check integrity

-  procedure Test (m: Multi_int; test_last: Boolean:= True );

-

-  -- i3 must be = i1 * i2

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

-

-  -- i1 must be = i2 * q + r

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

-

-end Multi_precision_integers.Check;

diff --git a/src/multi_precision_integers-io.adb b/src/multi_precision_integers-io.adb
deleted file mode 100644
index 5cdb1c4..0000000
--- a/src/multi_precision_integers-io.adb
+++ /dev/null
@@ -1,216 +0,0 @@
------------------------------------------------------------------------------

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

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

-

-

-

-

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

---

---  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.

---

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

-

-

-

-

-package body Multi_precision_integers.IO is

-

-  package IIO is new Integer_IO( index_int );

-

-  table: constant array(basic_int'(0)..15) of Character:=

-         ('0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F');

-

-   -- 15-Feb-2002: Bugfix case i=0. Spotted by Duncan Sands

-

-  function Chiffres_i_non_nul(i: multi_int; base: number_base:= 10) return Natural is

-    nombre: multi_int(i.last_used);

-    la_base    : constant basic_int :=  basic_int(base);

-    nchiffres: Natural:= 1;

-

-    procedure Comptage_rapide( C: Positive ) is

-      test  : multi_int(i.n);

-      base_puiss_C: constant multi_int:= Multi( basic_int(base) ) ** C;

-    begin

-      loop

-        Fill(test, nombre / base_puiss_C );

-        exit when test.zero;

-        -- quotient non nul, donc on a au moins C chiffres

-        Fill(nombre, test);

-        nchiffres:= nchiffres + C;

-      end loop;

-    end Comptage_rapide;

-

-  begin

-    Fill(nombre, i);

-    Comptage_rapide( 400 );

-    Comptage_rapide( 20 );

-    loop

-      Fill(nombre, nombre / la_base);

-      exit when nombre.zero;

-      nchiffres:= nchiffres + 1;

-    end loop;

-    return nchiffres;

-  end Chiffres_i_non_nul;

-

-  function Number_of_digits(i: multi_int; base: number_base:= 10) return Natural is

-  begin

-    if i.zero then

-      return 1;

-    else

-      return Chiffres_i_non_nul(i,base);

-    end if;

-  end Number_of_digits;

-

-  function Str(i: multi_int; base: number_base:= 10) return String is

-    res: String(1..1 + Number_of_digits(i,base)):= (others=> 'x');

-    nombre : multi_int(i.n):= i;

-    chiffre: basic_int;

-    la_base: constant basic_int :=  basic_int(base);

-

-  begin

-    if nombre.zero or else not nombre.neg then

-      res(1):= ' ';

-    else

-      res(1):= '-';

-    end if;

-    nombre.neg:= False;

-

-    -- maintenant nombre et base sont >=0, MOD=REM

-    for k in reverse 2 .. res'Last loop

-      Div_Rem( nombre, la_base, nombre, chiffre );

-      res(k):= table( chiffre );

-      exit when nombre.zero;

-    end loop;

-    return res;

-

-  end Str;

-

-

--- !!! recursion !!!

-

-  function Val(s: String) return multi_int is

-    formatting_error: exception;

-  begin

-    if s="" then

-      return Multi(0);

-    elsif s(s'First)='-' then

-      return -Val(s(s'First+1..s'Last));

-    elsif s(s'First)='+' then

-      return  Val(s(s'First+1..s'Last));

-    elsif s(s'Last) in '0'..'9' then

-      return basic_int'Value(s(s'Last..s'Last)) + 10 *

-             Val(s(s'First..s'Last-1));

-    else

-      raise formatting_error;

-    end if;

-  end Val;

-

-  procedure Put_in_blocks(File  : in File_Type;

-                          Item  : in multi_int) is

-  begin

-    if Item.neg then put(File,'-'); else put(File,'+'); end if;

-    Put(File, " [ ");

-    IIO.Put(File, 1+Item.n , 3);

-    Put(File, " blocks ]: ");

-    Put(File,'{');

-    if Item.n > Item.last_used then

-      IIO.Put(File, Item.n - Item.last_used, 3);

-      Put(File, " unused |");

-    end if;

-    for k in reverse 0 .. Item.last_used loop

-      Put(File, Block_type'Image(Item.blk(k)));

-      if k>0 then Put(File,'|'); end if;

-    end loop;

-    Put(File,'}');

-  end Put_in_blocks;

-

-  procedure Put_in_blocks(Item  : in multi_int) is

-  begin

-    Put_in_blocks( Standard_Output, Item );

-  end Put_in_blocks;

-

-  procedure Get(File  : in  File_Type;

-                Item  : out multi_int;

-                Width : in Field := 0) is

-  begin

-    null; -- !!!

-  end Get;

-

-  procedure Get(Item  : out multi_int;

-                Width : in  Field := 0) is

-

-  begin

-    Get(Standard_Input, Item, Width);

-  end Get;

-

-

-  procedure Put(File  : in File_Type;

-                Item  : in multi_int;

-                Width : in Field := 0;

-                Base  : in Number_Base := Default_Base) is

-

-  begin

-    if Width = 0 then       -- No padding required (default)

-      Put(File, Str(Item, Base));

-    else                    -- Left padding required -> slow

-      declare

-        la_chaine: String(1..Width);

-      begin

-        Put(la_chaine, Item, Base);

-        Put(File, la_chaine);

-      end;

-    end if;

-  end Put;

-

-  procedure Put(Item  : in multi_int;

-                Width : in Field := 0;

-                Base  : in Number_Base := Default_Base) is

-

-  begin

-    Put(Standard_Output, Item, Width, Base);

-  end Put;

-

-  procedure Get(From : in  String;

-                Item : out multi_int;

-                Last : out Positive) is

-  begin

-    Last:= 1;

-    null; -- !!!

-  end Get;

-

-

-  procedure Put(To   : out String;

-                Item : in multi_int;

-                Base : in Number_Base := Default_Base) is

-

-    nchiffres: constant Natural:= Number_of_digits(Item, Base);

-    blancs: constant String(To'Range):= (others=> ' ');

-  begin

-    if nchiffres > To'Length then

-      raise Layout_Error;

-    else

-      To:= blancs;

-      To( To'Last - nchiffres .. To'Last ):= Str(Item, Base);

-    end if;

-  end Put;

-

-end Multi_precision_integers.IO;

diff --git a/src/multi_precision_integers-io.ads b/src/multi_precision_integers-io.ads
deleted file mode 100644
index 558a62e..0000000
--- a/src/multi_precision_integers-io.ads
+++ /dev/null
@@ -1,94 +0,0 @@
-------------------------------------------------------------------------------

---  File:            Multi_precision_integers-IO.ads

---  Description:     Child of package 'Multi_precision_integers: I/O

---  Date/version:    2006 ; 15-Feb-2002 / 22.XI.1999 / 22.12.1996

---  Author:          Gautier de Montmollin

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

-

-

-

-

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

---

---  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 Text_IO;     use Text_IO;

-

-package Multi_precision_integers.IO is

-

-  Default_Base: Number_Base := 10;

-

-  -- Returns the number of digits in the specified base:

-  function Number_of_digits(i: Multi_int; base: Number_Base:= 10) return Natural;

-

-  -- Returns the image of i in the specified base:

-  function Str(i: Multi_int; base: Number_Base:= 10) return String;

-

-  -- Returns the value of number in string:

-  function Val(s: String) return Multi_int;

-

-  -- Output to file, in block format:

-  procedure Put_in_blocks(File  : in File_Type;

-                          Item  : in Multi_int);

-

-  -- Output to standard input, in block format:

-  procedure Put_in_blocks(Item  : in Multi_int);

-

-  ---- The following mimic the Text_IO.Integer_IO

-

-  -- Get from file:

-  procedure Get(File  : in  File_Type;

-                Item  : out Multi_int;

-                Width : in Field := 0);

-

-  -- Get from standard input:

-  procedure Get(Item  : out Multi_int;

-                Width : in  Field := 0);

-

-  -- Put to file:

-  procedure Put(File  : in File_Type;

-                Item  : in Multi_int;

-                Width : in Field := 0;

-                Base  : in Number_Base := Default_Base);

-  -- Width=0 : no default formatting, since inpredicatble length

-

-  -- Put to standard output:

-  procedure Put(Item  : in Multi_int;

-                Width : in Field := 0;

-                Base  : in Number_Base := Default_Base);

-  -- Width=0 : no default formatting, since inpredicatble length

-

-  -- Get from string:

-  procedure Get(From : in  String;

-                Item : out Multi_int;

-                Last : out Positive);

-

-  -- Put to string:

-  procedure Put(To   : out String;

-                Item : in Multi_int;

-                Base : in Number_Base := Default_Base);

-

-end Multi_precision_integers.IO;

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;

diff --git a/src/multi_precision_integers.ads b/src/multi_precision_integers.ads
deleted file mode 100644
index 7035924..0000000
--- a/src/multi_precision_integers.ads
+++ /dev/null
@@ -1,266 +0,0 @@
-------------------------------------------------------------------------------

---  File:            multi_precision_integers.ads

---

---  Description:     Multiple precision integers package

---

---  Date/version:       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: - unsigned types for blocks

---                                - a block uses 2 bits more

---                                - Ada95+ only

---                      Mar-2002: Bugs fixed related to zero field, division

---                      Nov-1999: - procedures (no stack, less copies !)

---                                - new data structure

---                      Dec-1996: First version (operators only)

---

---  Author:          G. de Montmollin, Univ. Neuchatel

---

---  Thanks to:       Duncan Sands, Univ. Paris-Sud, CNRS

---                   Jeffrey R. Carter

---

---  Tested on:       Intel 586 (32 bit) - Windows 98, NT4+ - GNAT 3.13p+

---                   Intel 586 (32 bit) - Windows 98, NT4+ - ObjectAda 7.2.1+

---                   Alpha-AXP (64 bit) - OpenVMS 7.1      - Compaq Ada

---

---  Division algorithm adaptated from BigInt 1.0 library,

---  by Stephen Adams, that refers to

---  D. E. Knuth, the Art of computer programming

---  volume 2, "Seminumerical Algorithms"

---  section 4.3.1, "Multiple-Precision Arithmetic"

---

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

-

-

-

-

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

---

---  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 System;

-

-package Multi_precision_integers is

-

-  -- Integers for array indexing --

-

-  subtype Index_int is Integer;

-

-  -- THE multi-precision integer type --

-

-  type Multi_int (n: Index_int) is private;

-

-  -- Integer type for small values --

-

-  subtype Basic_int is Integer; -- the "normal" signed integer

-

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

-  -- Debug mode: checks results of arithmetic operations and    --

-  -- Multi_int variables' integrity.                            --

-  -- CAUTION: Debug = True reduces monstruously the performance --

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

-

-  Debug : constant Boolean:= False;

-

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

-  ----- Informations, conversions, filling ----

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

-

-  -- Convert Basic_int to Multi_int

-  function Multi (small: Basic_int) return Multi_int;

-

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

-  function Basic (large: Multi_int) return Basic_int;

-

-  -- Fill an Multi_int of greater array dimension with a smaller one

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

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

-

-  -- Comparisons

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

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

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

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

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

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

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

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

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

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

-

-  -- Other informations

-  function Bits_per_block return Positive;

-

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

-  -------- Arithmetic operators.                                   ----------

-  -------- For speed, the "procedure" variants should be preferred ----------

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

-

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

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

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

-

-  procedure Opp (i: in out Multi_int);

-  function "+" (i : Multi_int) return Multi_int;

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

-

-  procedure Abso (i: in out Multi_int);

-  function "ABS" (i : Multi_int) return Multi_int;

-

-  function Sign (i: Multi_int) return Basic_int;

-  function Even (i: Multi_int) return Boolean;

-  function Odd (i : Multi_int) return Boolean;

-

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

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

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

-

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

-  -- Addition, subtraction --

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

-

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

-

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

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

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

-

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

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

-    renames Sub;

-

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

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

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

-

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

-  -- Multiplication --

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

-

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

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

-    renames Multiply;

-

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

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

-    renames Multiply;

-

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

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

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

-

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

-  -- Division, Remainder --

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

-

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

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

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

-

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

-

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

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

-  function "Rem" (i1, i2: Multi_int) return Multi_int;

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

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

-  function "Mod" (i1, i2: Multi_int) return Multi_int;

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

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

-

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

-  -- Power --

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

-

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

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

-  -- + 26-Mar-2002 :

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

-                   modulo : Multi_int);

-

-  Cannot_fit, Empty_multi_int : exception;

-

-  Array_too_small : exception;

-

-  Result_undersized,

-  Quotient_undersized,

-  Remainder_undersized : exception;

-

-  Division_by_zero : exception;

-

-  Zero_power_zero, Power_negative : exception;

-  Power_modulo_non_positive : exception;

-

-private

-

-  --> Long_Block_type is used for + and * of blocks.

-  -- It is by design

-  --   a/ the largest possible modular integer, and

-  --   b/ twice the size of Block_type defined below.

-  -- With the double of bits (2n) one can store m**2 + 2m without overflow

-  --   where m = 2**n - 1 is the largest value possible on n bits.

-  type Long_Block_type is mod System.Max_Binary_Modulus;

-

-  --> Same size as Long_Block_type, but signed:

-  type Long_Block_type_signed is

-    range -2**(Long_Block_type'Size-1)..2**(Long_Block_type'Size-1)-1;

-

-  --> Block_type: unsigned integer used to store a chunk of a Multi_int.

-  type Block_type is mod 2 ** (Long_Block_type'Size / 2);

-

-  Block_type_bits: constant:= Block_type'Size;

-

-  cardblock: constant:= 2 ** Block_type_bits;

-  -- Number of possible values

-  maxblock: constant:= Block_type'Last;

-  -- NB: GNAT (2006) optimizes out correctly the Block_type(l and maxblock_long)

-  --     to Block_type(l), the latter form being caught when range checks are on.

-

-  type Block_array is array( Index_int range <> ) of Block_type;

-  -- 2006: was "of Basic_int" for obscure reasons...

-

-  type Multi_int(n: Index_int) is record

-    blk:       Block_array( 0..n ); -- the n blocks with ABSOLUTE value

-    neg:       Boolean;             -- negative flag

-    zero:      Boolean:=True;       -- zero flag (supercedes the other fields)

-    last_used: Index_int;           -- the others blocks are supposed 0

-  end record;

-

-  -- NB the `zero' field supercedes EVERY other information (last_used, neg)

-

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

-  --   Format of type Multi_int.blk: ( i_0, i_1, ..., i_k, *, ..., * )      --

-  --   i_0..i_k are >=0 ; others (*) are treated as 0                       --

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

-

-  -- Some internal procedures use by check:

-

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

-

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

-

-end Multi_precision_integers;

diff --git a/src/rationals.adb b/src/rationals.adb
index 3e3cd88..80e070a 100644
--- a/src/rationals.adb
+++ b/src/rationals.adb
@@ -1,10 +1,5 @@
 
 
-with Ada.Strings.Fixed;
-with Multi_Precision_Integers.IO;
-use Multi_Precision_Integers.IO;
-
-
 --  This source is licensed under Creative Commons CC0 v1.0.
 --
 --  To read the full text, see license.txt in the main directory of this repository
@@ -13,31 +8,173 @@ use Multi_Precision_Integers.IO;
 --  For a human readable summary, go to https://creativecommons.org/publicdomain/zero/1.0/
 
 
+with
+
+    Interfaces.C;
+
+use type
+
+    Interfaces.C.int;
+
+
 package body Rationals is
 
 
-    function Reduce
-           (Numerator, Denominator : in Multi_Int)
-        return Fraction
-    is
-        A, B, Temp : Multi_Int (M_Size);
-        Ret_N, Ret_D : Multi_Int (M_Size);
+    procedure mpq_canonicalize
+           (Op : in out mpq_t)
+        with Inline => True;
+    pragma Import (C, mpq_canonicalize, "__gmpq_canonicalize");
+
+    procedure mpq_init
+           (X : out mpq_t)
+        with Inline => True;
+    pragma Import (C, mpq_init, "__gmpq_init");
+
+    procedure mpq_clear
+           (X : in out mpq_t)
+        with Inline => True;
+    pragma Import (C, mpq_clear, "__gmpq_clear");
+
+    procedure mpq_set
+           (Rop : in out mpq_t;
+            Op  : in     mpq_t)
+        with Inline => True;
+    pragma Import (C, mpq_set, "__gmpq_set");
+
+    procedure mpq_set_si
+           (Rop : in out mpq_t;
+            Op1 : in     Interfaces.C.long;
+            Op2 : in     Interfaces.C.unsigned_long)
+        with Inline => True;
+    pragma Import (C, mpq_set_si, "__gmpq_set_si");
+
+
+
+
+    procedure mpq_add
+           (Sum              : in out mpq_t;
+            Addend1, Addend2 : in     mpq_t)
+        with Inline => True;
+    pragma Import (C, mpq_add, "__gmpq_add");
+
+    procedure mpq_sub
+           (Difference          : in out mpq_t;
+            Minuend, Subtrahend : in     mpq_t)
+        with Inline => True;
+    pragma Import (C, mpq_sub, "__gmpq_sub");
+
+    procedure mpq_mul
+           (Product                  : in out mpq_t;
+            Multiplier, Multiplicand : in     mpq_t)
+        with Inline => True;
+    pragma Import (C, mpq_mul, "__gmpq_mul");
+
+    procedure mpq_div
+           (Quotient          : in out mpq_t;
+            Dividend, Divisor : in     mpq_t)
+        with Inline => True;
+    pragma Import (C, mpq_div, "__gmpq_div");
+
+    procedure mpq_neg
+           (Negated : in out mpq_t;
+            Operand : in     mpq_t)
+        with Inline => True;
+    pragma Import (C, mpq_neg, "__gmpq_neg");
+
+
+
+
+    function mpq_cmp
+           (Op1, Op2 : in mpq_t)
+        return Interfaces.C.int
+        with Inline => True;
+    pragma Import (C, mpq_cmp, "__gmpq_cmp");
+
+    function mpq_equal
+           (Op1, Op2 : in mpq_t)
+        return Interfaces.C.int
+        with Inline => True;
+    pragma Import (C, mpq_equal, "__gmpq_equal");
+
+
+
+
+    procedure mpz_init
+           (X : out mpz_t)
+        with Inline => True;
+    pragma Import (C, mpz_init, "__gmpz_init");
+
+    procedure mpz_clear
+           (X : in out mpz_t)
+        with Inline => True;
+    pragma Import (C, mpz_clear, "__gmpz_clear");
+
+    function mpz_get_si
+           (Op : in mpz_t)
+        return Interfaces.C.long
+        with Inline => True;
+    pragma Import (C, mpz_get_si, "__gmpz_get_si");
+
+    procedure mpz_set_si
+           (Rop : in out mpz_t;
+            Op  : in     Interfaces.C.long)
+        with Inline => True;
+    pragma Import (C, mpz_set_si, "__gmpz_set_si");
+
+
+
+
+    procedure mpz_cdiv_q
+           (Q    : in out mpz_t;
+            N, R : in     mpz_t)
+        with Inline => True;
+    pragma Import (C, mpz_cdiv_q, "__gmpz_cdiv_q");
+
+    procedure mpz_fdiv_q
+           (Q    : in out mpz_t;
+            N, R : in     mpz_t)
+        with Inline => True;
+    pragma Import (C, mpz_fdiv_q, "__gmpz_fdiv_q");
+
+    procedure mpz_mod
+           (R    : in out mpz_t;
+            N, D : in     mpz_t)
+        with Inline => True;
+    pragma Import (C, mpz_mod, "__gmpz_mod");
+
+
+
+
+    function mpz_cmp
+           (Op1, Op2 : in mpz_t)
+        return Interfaces.C.int
+        with Inline => True;
+    pragma Import (C, mpz_cmp, "__gmpz_cmp");
+
+
+
+
+    procedure Initialize
+           (This : in out Fraction) is
     begin
-        Fill (A, Numerator);
-        Fill (B, Denominator);
+        mpq_init (This.Data);
+    end Initialize;
 
-        --  Euclid's algorithm
-        loop
-            Temp := A;
-            A := B;
-            Fill (B, Temp mod B);
-            exit when Equal (B, 0);
-        end loop;
+    procedure Adjust
+           (This : in out Fraction)
+    is
+        Temp : mpq_t;
+    begin
+        mpq_init (Temp);
+        mpq_set (Temp, This.Data);
+        This.Data := Temp;
+    end Adjust;
 
-        Fill (Ret_N, Numerator / A);
-        Fill (Ret_D, Denominator / A);
-        return (Num => Ret_N, Den => Ret_D);
-    end Reduce;
+    procedure Finalize
+           (This : in out Fraction) is
+    begin
+        mpq_clear (This.Data);
+    end Finalize;
 
 
 
@@ -46,29 +183,35 @@ package body Rationals is
            (Left, Right : in Fraction)
         return Fraction is
     begin
-        return Reduce
-               (Left.Num * Right.Den + Left.Den * Right.Num,
-                Left.Den * Right.Den);
+        return Result : Fraction do
+            mpq_add (Result.Data, Left.Data, Right.Data);
+        end return;
     end "+";
 
     function "+"
            (Left  : in Fraction;
             Right : in Integer)
-        return Fraction is
+        return Fraction
+    is
+        Temp : Fraction;
     begin
-        return Reduce
-               (Left.Num + Left.Den * Multi (Right),
-                Left.Den);
+        return Result : Fraction do
+            mpq_set_si (Temp.Data, Interfaces.C.long (Right), 1);
+            mpq_add (Result.Data, Left.Data, Temp.Data);
+        end return;
     end "+";
 
     function "+"
            (Left  : in Integer;
             Right : in Fraction)
-        return Fraction is
+        return Fraction
+    is
+        Temp : Fraction;
     begin
-        return Reduce
-               (Multi (Left) * Right.Den + Right.Num,
-                Right.Den);
+        return Result : Fraction do
+            mpq_set_si (Temp.Data, Interfaces.C.long (Left), 1);
+            mpq_add (Result.Data, Temp.Data, Right.Data);
+        end return;
     end "+";
 
 
@@ -78,39 +221,44 @@ package body Rationals is
            (Left, Right : in Fraction)
         return Fraction is
     begin
-        return Reduce
-               (Left.Num * Right.Den - Left.Den * Right.Num,
-                Left.Den * Right.Den);
+        return Result : Fraction do
+            mpq_sub (Result.Data, Left.Data, Right.Data);
+        end return;
     end "-";
 
     function "-"
            (Left  : in Fraction;
             Right : in Integer)
-        return Fraction is
+        return Fraction
+    is
+        Temp : Fraction;
     begin
-        return Reduce
-               (Left.Num - Left.Den * Multi (Right),
-                Left.Den);
+        return Result : Fraction do
+            mpq_set_si (Temp.Data, Interfaces.C.long (Right), 1);
+            mpq_sub (Result.Data, Left.Data, Temp.Data);
+        end return;
     end "-";
 
     function "-"
            (Left  : in Integer;
             Right : in Fraction)
-        return Fraction is
+        return Fraction
+    is
+        Temp : Fraction;
     begin
-        return Reduce
-               (Multi (Left) * Right.Den - Right.Num,
-                Right.Den);
+        return Result : Fraction do
+            mpq_set_si (Temp.Data, Interfaces.C.long (Left), 1);
+            mpq_sub (Result.Data, Temp.Data, Right.Data);
+        end return;
     end "-";
 
-
-
-
     function "-"
            (Right : in Fraction)
         return Fraction is
     begin
-        return (Num => - Right.Num, Den => Right.Den);
+        return Result : Fraction do
+            mpq_neg (Result.Data, Right.Data);
+        end return;
     end "-";
 
 
@@ -120,29 +268,35 @@ package body Rationals is
            (Left, Right : in Fraction)
         return Fraction is
     begin
-        return Reduce
-               (Left.Num * Right.Num,
-                Left.Den * Right.Den);
+        return Result : Fraction do
+            mpq_mul (Result.Data, Left.Data, Right.Data);
+        end return;
     end "*";
 
     function "*"
            (Left  : in Fraction;
             Right : in Integer)
-        return Fraction is
+        return Fraction
+    is
+        Temp : Fraction;
     begin
-        return Reduce
-               (Left.Num * Multi (Right),
-                Left.Den);
+        return Result : Fraction do
+            mpq_set_si (Temp.Data, Interfaces.C.long (Right), 1);
+            mpq_mul (Result.Data, Left.Data, Temp.Data);
+        end return;
     end "*";
 
     function "*"
            (Left  : in Integer;
             Right : in Fraction)
-        return Fraction is
+        return Fraction
+    is
+        Temp : Fraction;
     begin
-        return Reduce
-               (Multi (Left) * Right.Num,
-                Right.Den);
+        return Result : Fraction do
+            mpq_set_si (Temp.Data, Interfaces.C.long (Left), 1);
+            mpq_mul (Result.Data, Temp.Data, Right.Data);
+        end return;
     end "*";
 
 
@@ -152,36 +306,48 @@ package body Rationals is
            (Left, Right : in Fraction)
         return Fraction is
     begin
-        return Reduce
-               (Left.Num * Right.Den,
-                Left.Den * Right.Num);
+        return Result : Fraction do
+            mpq_div (Result.Data, Left.Data, Right.Data);
+        end return;
     end "/";
 
     function "/"
            (Left  : in Fraction;
             Right : in Integer)
-        return Fraction is
+        return Fraction
+    is
+        Temp : Fraction;
     begin
-        return Reduce
-               (Left.Num,
-                Left.Den * Multi (Right));
+        return Result : Fraction do
+            mpq_set_si (Temp.Data, Interfaces.C.long (Right), 1);
+            mpq_div (Result.Data, Left.Data, Temp.Data);
+        end return;
     end "/";
 
     function "/"
            (Left  : in Integer;
             Right : in Fraction)
-        return Fraction is
+        return Fraction
+    is
+        Temp : Fraction;
     begin
-        return Reduce
-               (Right.Num,
-                Multi (Left) * Right.Den);
+        return Result : Fraction do
+            mpq_set_si (Temp.Data, Interfaces.C.long (Left), 1);
+            mpq_div (Result.Data, Temp.Data, Right.Data);
+        end return;
     end "/";
 
     function "/"
            (Left, Right : in Integer)
         return Fraction is
     begin
-        return Reduce (Multi (Left), Multi (Right));
+        return Result : Fraction do
+            mpq_set_si
+                   (Result.Data,
+                    Interfaces.C.long (Left),
+                    Interfaces.C.unsigned_long (Right));
+            mpq_canonicalize (Result.Data);
+        end return;
     end "/";
 
 
@@ -191,24 +357,29 @@ package body Rationals is
            (Left, Right : in Fraction)
         return Boolean is
     begin
-        return  Equal (Left.Num, Right.Num) and
-                Equal (Left.Den, Right.Den);
+        return mpq_equal (Left.Data, Right.Data) /= 0;
     end "=";
 
     function "="
            (Left  : in Fraction;
             Right : in Integer)
-        return Boolean is
+        return Boolean
+    is
+        Temp : Fraction;
     begin
-        return Equal (Left.Num, Right) and Equal (Left.Den, 1);
+        mpq_set_si (Temp.Data, Interfaces.C.long (Right), 1);
+        return mpq_equal (Left.Data, Temp.Data) /= 0;
     end "=";
 
     function "="
            (Left  : in Integer;
             Right : in Fraction)
-        return Boolean is
+        return Boolean
+    is
+        Temp : Fraction;
     begin
-        return Equal (Right.Num, Left) and Equal (Right.Den, 1);
+        mpq_set_si (Temp.Data, Interfaces.C.long (Left), 1);
+        return mpq_equal (Temp.Data, Right.Data) /= 0;
     end "=";
 
 
@@ -218,23 +389,29 @@ package body Rationals is
            (Left, Right : in Fraction)
         return Boolean is
     begin
-        return Left.Num * Right.Den <= Left.Den * Right.Num;
+        return mpq_cmp (Left.Data, Right.Data) <= 0;
     end "<=";
 
     function "<="
            (Left  : in Fraction;
             Right : in Integer)
-        return Boolean is
+        return Boolean
+    is
+        Temp : Fraction;
     begin
-        return Left.Num <= Left.Den * Multi (Right);
+        mpq_set_si (Temp.Data, Interfaces.C.long (Right), 1);
+        return mpq_cmp (Left.Data, Temp.Data) <= 0;
     end "<=";
 
     function "<="
            (Left  : in Integer;
             Right : in Fraction)
-        return Boolean is
+        return Boolean
+    is
+        Temp : Fraction;
     begin
-        return Multi (Left) * Right.Den <= Right.Num;
+        mpq_set_si (Temp.Data, Interfaces.C.long (Left), 1);
+        return mpq_cmp (Temp.Data, Right.Data) <= 0;
     end "<=";
 
 
@@ -244,23 +421,29 @@ package body Rationals is
            (Left, Right : in Fraction)
         return Boolean is
     begin
-        return Left.Num * Right.Den < Left.Den * Right.Num;
+        return mpq_cmp (Left.Data, Right.Data) < 0;
     end "<";
 
     function "<"
            (Left  : in Fraction;
             Right : in Integer)
-        return Boolean is
+        return Boolean
+    is
+        Temp : Fraction;
     begin
-        return Left.Num < Left.Den * Multi (Right);
+        mpq_set_si (Temp.Data, Interfaces.C.long (Right), 1);
+        return mpq_cmp (Left.Data, Temp.Data) < 0;
     end "<";
 
     function "<"
            (Left  : in Integer;
             Right : in Fraction)
-        return Boolean is
+        return Boolean
+    is
+        Temp : Fraction;
     begin
-        return Multi (Left) * Right.Den < Right.Num;
+        mpq_set_si (Temp.Data, Interfaces.C.long (Left), 1);
+        return mpq_cmp (Temp.Data, Right.Data) < 0;
     end "<";
 
 
@@ -270,23 +453,29 @@ package body Rationals is
            (Left, Right : in Fraction)
         return Boolean is
     begin
-        return Left.Num * Right.Den >= Left.Den * Right.Num;
+        return mpq_cmp (Left.Data, Right.Data) >= 0;
     end ">=";
 
     function ">="
            (Left  : in Fraction;
             Right : in Integer)
-        return Boolean is
+        return Boolean
+    is
+        Temp : Fraction;
     begin
-        return Left.Num >= Left.Den * Multi (Right);
+        mpq_set_si (Temp.Data, Interfaces.C.long (Right), 1);
+        return mpq_cmp (Left.Data, Temp.Data) >= 0;
     end ">=";
 
     function ">="
            (Left  : in Integer;
             Right : in Fraction)
-        return Boolean is
+        return Boolean
+    is
+        Temp : Fraction;
     begin
-        return Multi (Left) * Right.Den >= Right.Num;
+        mpq_set_si (Temp.Data, Interfaces.C.long (Left), 1);
+        return mpq_cmp (Temp.Data, Right.Data) >= 0;
     end ">=";
 
 
@@ -296,23 +485,29 @@ package body Rationals is
            (Left, Right : in Fraction)
         return Boolean is
     begin
-        return Left.Num * Right.Den > Left.Den * Right.Num;
+        return mpq_cmp (Left.Data, Right.Data) > 0;
     end ">";
 
     function ">"
            (Left  : in Fraction;
             Right : in Integer)
-        return Boolean is
+        return Boolean
+    is
+        Temp : Fraction;
     begin
-        return Left.Num > Left.Den * Multi (Right);
+        mpq_set_si (Temp.Data, Interfaces.C.long (Right), 1);
+        return mpq_cmp (Left.Data, Temp.Data) > 0;
     end ">";
 
     function ">"
            (Left  : in Integer;
             Right : in Fraction)
-        return Boolean is
+        return Boolean
+    is
+        Temp : Fraction;
     begin
-        return Multi (Left) * Right.Den > Right.Num;
+        mpq_set_si (Temp.Data, Interfaces.C.long (Left), 1);
+        return mpq_cmp (Temp.Data, Right.Data) > 0;
     end ">";
 
 
@@ -320,72 +515,73 @@ package body Rationals is
 
     function Numerator
            (Item : in Fraction)
-        return Integer is
+        return Long_Integer is
     begin
-        return Basic (Item.Num);
+        return Long_Integer (mpz_get_si (Item.Data.mp_num));
     end Numerator;
 
+
     function Denominator
            (Item : in Fraction)
-        return Integer is
+        return Long_Integer is
     begin
-        return Basic (Item.Den);
+        return Long_Integer (mpz_get_si (Item.Data.mp_den));
     end Denominator;
 
+
     function Floor
            (Item : in Fraction)
-        return Integer is
-    begin
-        return Basic (Item.Num / Item.Den);
+        return Long_Integer
+    is
+        Temp   : mpz_t;
+        Result : Long_Integer;
+    begin
+        mpz_init (Temp);
+        mpz_fdiv_q (Temp, Item.Data.mp_num, Item.Data.mp_den);
+        Result := Long_Integer (mpz_get_si (Temp));
+        mpz_clear (Temp);
+        return Result;
     end Floor;
 
+
     function Ceiling
            (Item : in Fraction)
-        return Integer is
-    begin
-        if Equal (Item.Num mod Item.Den, 0) then
-            return Basic (Item.Num / Item.Den);
-        else
-            return 1 + Basic (Item.Num / Item.Den);
-        end if;
+        return Long_Integer
+    is
+        Temp   : mpz_t;
+        Result : Long_Integer;
+    begin
+        mpz_init (Temp);
+        mpz_cdiv_q (Temp, Item.Data.mp_num, Item.Data.mp_den);
+        Result := Long_Integer (mpz_get_si (Temp));
+        mpz_clear (Temp);
+        return Result;
     end Ceiling;
 
+
     function Round
            (Item : in Fraction)
-        return Integer is
-    begin
-        if Item.Num mod Item.Den >= Item.Den / 2 then
-            return 1 + Basic (Item.Num / Item.Den);
-        else
-            return Basic (Item.Num / Item.Den);
+        return Long_Integer
+    is
+        Temp1, Temp2, Temp3 : mpz_t;
+        Result : Long_Integer;
+    begin
+        mpz_init (Temp1); mpz_init (Temp2); mpz_init (Temp3);
+        --  This method of calculating whether Num mod Den >= Den / 2 is messy
+        mpz_mod (Temp1, Item.Data.mp_num, Item.Data.mp_den);
+        mpz_set_si (Temp3, 2);
+        mpz_fdiv_q (Temp2, Item.Data.mp_den, Temp3);
+        --  Here the actual Num / Den division takes place
+        mpz_fdiv_q (Temp3, Item.Data.mp_num, Item.Data.mp_den);
+        Result := Long_Integer (mpz_get_si (Temp3));
+        if mpz_cmp (Temp1, Temp2) >= 0 then
+            Result := Result + 1;
         end if;
+        mpz_clear (Temp1); mpz_clear (Temp2); mpz_clear (Temp3);
+        return Result;
     end Round;
 
 
-
-
-    function Image
-           (Item : in Fraction)
-        return String is
-    begin
-        return Str (Item.Num) & '/' & Str (Item.Den);
-    end Image;
-
-    function Value
-           (Item : in String)
-        return Fraction
-    is
-        use Ada.Strings;
-        use Ada.Strings.Fixed;
-        A, B, S : Integer;
-    begin
-        S := Index (Item, "/");
-        A := Integer'Value (Item (Item'First .. S - 1));
-        B := Integer'Value (Item (S + 1 .. Item'Last));
-        return Reduce (Multi (A), Multi (B));
-    end Value;
-
-
 end Rationals;
 
 
diff --git a/src/rationals.ads b/src/rationals.ads
index f12f161..f907633 100644
--- a/src/rationals.ads
+++ b/src/rationals.ads
@@ -1,8 +1,5 @@
 
 
-private with Multi_Precision_Integers;
-
-
 --  This source is licensed under Creative Commons CC0 v1.0.
 --
 --  To read the full text, see license.txt in the main directory of this repository
@@ -11,6 +8,12 @@ private with Multi_Precision_Integers;
 --  For a human readable summary, go to https://creativecommons.org/publicdomain/zero/1.0/
 
 
+private with
+
+    Ada.Finalization,
+    Interfaces.C;
+
+
 package Rationals is
 
 
@@ -185,49 +188,63 @@ package Rationals is
 
     function Numerator
            (Item : in Fraction)
-        return Integer;
+        return Long_Integer;
 
     function Denominator
            (Item : in Fraction)
-        return Integer;
+        return Long_Integer;
 
     function Floor
            (Item : in Fraction)
-        return Integer;
+        return Long_Integer;
 
     function Ceiling
            (Item : in Fraction)
-        return Integer;
+        return Long_Integer;
 
     function Round
            (Item : in Fraction)
-        return Integer;
+        return Long_Integer;
 
 
+private
 
 
-    function Image
-           (Item : in Fraction)
-        return String;
+    --  Types necessary to bind to GMP
+    type mp_limb_t is new Interfaces.C.unsigned;
+    type mp_ptr is access mp_limb_t;
 
-    function Value
-           (Item : in String)
-        return Fraction;
+    type mpz_t is record
+        mp_alloc, mp_size : Interfaces.C.int;
+        mp_d              : mp_ptr;
+    end record;
 
+    type mpq_t is record
+        mp_num, mp_den : mpz_t;
+    end record;
 
-private
 
 
-    use Multi_Precision_Integers;
 
+    --  The actual type exported by this package, which
+    --  has to be Controlled so it gets deallocated properly
+    type Fraction is new Ada.Finalization.Controlled with record
+        Data : mpq_t;
+    end record;
 
-    M_Size : constant Basic_Int := 4;
+    overriding procedure Initialize
+           (This : in out Fraction);
 
+    overriding procedure Adjust
+           (This : in out Fraction);
 
-    type Fraction is record
-        Num : Multi_Int (M_Size);
-        Den : Multi_Int (M_Size);
-    end record;
+    overriding procedure Finalize
+           (This : in out Fraction);
+
+
+
+
+    pragma Linker_Options ("-lgmp");
 
 
 end Rationals;