summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/bundles.adb6
-rw-r--r--src/multi_precision_integers-check.adb170
-rw-r--r--src/multi_precision_integers-check.ads41
-rw-r--r--src/multi_precision_integers-io.adb216
-rw-r--r--src/multi_precision_integers-io.ads94
-rw-r--r--src/multi_precision_integers.adb1530
-rw-r--r--src/multi_precision_integers.ads266
-rw-r--r--src/rationals.adb482
-rw-r--r--src/rationals.ads59
9 files changed, 381 insertions, 2483 deletions
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;