From 45c752c00a8befa4041b4dcd8243b0563779c573 Mon Sep 17 00:00:00 2001 From: Jed Barber Date: Sun, 2 Jul 2017 23:39:38 +1000 Subject: Removed MIT licensed Mathpaqs packages in favour of bindings to GMP --- cc0.txt | 121 --- license.txt | 121 +++ mit.txt | 19 - notes.txt | 2 - readme.txt | 10 +- src/bundles.adb | 6 +- src/multi_precision_integers-check.adb | 170 ---- src/multi_precision_integers-check.ads | 41 - src/multi_precision_integers-io.adb | 216 ----- src/multi_precision_integers-io.ads | 94 -- src/multi_precision_integers.adb | 1530 -------------------------------- src/multi_precision_integers.ads | 266 ------ src/rationals.adb | 482 +++++++--- src/rationals.ads | 59 +- 14 files changed, 503 insertions(+), 2634 deletions(-) delete mode 100644 cc0.txt create mode 100644 license.txt delete mode 100644 mit.txt delete mode 100644 src/multi_precision_integers-check.adb delete mode 100644 src/multi_precision_integers-check.ads delete mode 100644 src/multi_precision_integers-io.adb delete mode 100644 src/multi_precision_integers-io.ads delete mode 100644 src/multi_precision_integers.adb delete mode 100644 src/multi_precision_integers.ads diff --git a/cc0.txt b/cc0.txt deleted file mode 100644 index 0e259d4..0000000 --- a/cc0.txt +++ /dev/null @@ -1,121 +0,0 @@ -Creative Commons Legal Code - -CC0 1.0 Universal - - CREATIVE COMMONS CORPORATION IS NOT A LAW FIRM AND DOES NOT PROVIDE - LEGAL SERVICES. DISTRIBUTION OF THIS DOCUMENT DOES NOT CREATE AN - ATTORNEY-CLIENT RELATIONSHIP. CREATIVE COMMONS PROVIDES THIS - INFORMATION ON AN "AS-IS" BASIS. CREATIVE COMMONS MAKES NO WARRANTIES - REGARDING THE USE OF THIS DOCUMENT OR THE INFORMATION OR WORKS - PROVIDED HEREUNDER, AND DISCLAIMS LIABILITY FOR DAMAGES RESULTING FROM - THE USE OF THIS DOCUMENT OR THE INFORMATION OR WORKS PROVIDED - HEREUNDER. - -Statement of Purpose - -The laws of most jurisdictions throughout the world automatically confer -exclusive Copyright and Related Rights (defined below) upon the creator -and subsequent owner(s) (each and all, an "owner") of an original work of -authorship and/or a database (each, a "Work"). - -Certain owners wish to permanently relinquish those rights to a Work for -the purpose of contributing to a commons of creative, cultural and -scientific works ("Commons") that the public can reliably and without fear -of later claims of infringement build upon, modify, incorporate in other -works, reuse and redistribute as freely as possible in any form whatsoever -and for any purposes, including without limitation commercial purposes. -These owners may contribute to the Commons to promote the ideal of a free -culture and the further production of creative, cultural and scientific -works, or to gain reputation or greater distribution for their Work in -part through the use and efforts of others. - -For these and/or other purposes and motivations, and without any -expectation of additional consideration or compensation, the person -associating CC0 with a Work (the "Affirmer"), to the extent that he or she -is an owner of Copyright and Related Rights in the Work, voluntarily -elects to apply CC0 to the Work and publicly distribute the Work under its -terms, with knowledge of his or her Copyright and Related Rights in the -Work and the meaning and intended legal effect of CC0 on those rights. - -1. Copyright and Related Rights. A Work made available under CC0 may be -protected by copyright and related or neighboring rights ("Copyright and -Related Rights"). Copyright and Related Rights include, but are not -limited to, the following: - - i. the right to reproduce, adapt, distribute, perform, display, - communicate, and translate a Work; - ii. moral rights retained by the original author(s) and/or performer(s); -iii. publicity and privacy rights pertaining to a person's image or - likeness depicted in a Work; - iv. rights protecting against unfair competition in regards to a Work, - subject to the limitations in paragraph 4(a), below; - v. rights protecting the extraction, dissemination, use and reuse of data - in a Work; - vi. database rights (such as those arising under Directive 96/9/EC of the - European Parliament and of the Council of 11 March 1996 on the legal - protection of databases, and under any national implementation - thereof, including any amended or successor version of such - directive); and -vii. other similar, equivalent or corresponding rights throughout the - world based on applicable law or treaty, and any national - implementations thereof. - -2. Waiver. To the greatest extent permitted by, but not in contravention -of, applicable law, Affirmer hereby overtly, fully, permanently, -irrevocably and unconditionally waives, abandons, and surrenders all of -Affirmer's Copyright and Related Rights and associated claims and causes -of action, whether now known or unknown (including existing as well as -future claims and causes of action), in the Work (i) in all territories -worldwide, (ii) for the maximum duration provided by applicable law or -treaty (including future time extensions), (iii) in any current or future -medium and for any number of copies, and (iv) for any purpose whatsoever, -including without limitation commercial, advertising or promotional -purposes (the "Waiver"). Affirmer makes the Waiver for the benefit of each -member of the public at large and to the detriment of Affirmer's heirs and -successors, fully intending that such Waiver shall not be subject to -revocation, rescission, cancellation, termination, or any other legal or -equitable action to disrupt the quiet enjoyment of the Work by the public -as contemplated by Affirmer's express Statement of Purpose. - -3. Public License Fallback. Should any part of the Waiver for any reason -be judged legally invalid or ineffective under applicable law, then the -Waiver shall be preserved to the maximum extent permitted taking into -account Affirmer's express Statement of Purpose. In addition, to the -extent the Waiver is so judged Affirmer hereby grants to each affected -person a royalty-free, non transferable, non sublicensable, non exclusive, -irrevocable and unconditional license to exercise Affirmer's Copyright and -Related Rights in the Work (i) in all territories worldwide, (ii) for the -maximum duration provided by applicable law or treaty (including future -time extensions), (iii) in any current or future medium and for any number -of copies, and (iv) for any purpose whatsoever, including without -limitation commercial, advertising or promotional purposes (the -"License"). The License shall be deemed effective as of the date CC0 was -applied by Affirmer to the Work. Should any part of the License for any -reason be judged legally invalid or ineffective under applicable law, such -partial invalidity or ineffectiveness shall not invalidate the remainder -of the License, and in such case Affirmer hereby affirms that he or she -will not (i) exercise any of his or her remaining Copyright and Related -Rights in the Work or (ii) assert any associated claims and causes of -action with respect to the Work, in either case contrary to Affirmer's -express Statement of Purpose. - -4. Limitations and Disclaimers. - - a. No trademark or patent rights held by Affirmer are waived, abandoned, - surrendered, licensed or otherwise affected by this document. - b. Affirmer offers the Work as-is and makes no representations or - warranties of any kind concerning the Work, express, implied, - statutory or otherwise, including without limitation warranties of - title, merchantability, fitness for a particular purpose, non - infringement, or the absence of latent or other defects, accuracy, or - the present or absence of errors, whether or not discoverable, all to - the greatest extent permissible under applicable law. - c. Affirmer disclaims responsibility for clearing rights of other persons - that may apply to the Work or any use thereof, including without - limitation any person's Copyright and Related Rights in the Work. - Further, Affirmer disclaims responsibility for obtaining any necessary - consents, permissions or other rights required for any use of the - Work. - d. Affirmer understands and acknowledges that Creative Commons is not a - party to this document and has no duty or obligation with respect to - this CC0 or use of the Work. diff --git a/license.txt b/license.txt new file mode 100644 index 0000000..0e259d4 --- /dev/null +++ b/license.txt @@ -0,0 +1,121 @@ +Creative Commons Legal Code + +CC0 1.0 Universal + + CREATIVE COMMONS CORPORATION IS NOT A LAW FIRM AND DOES NOT PROVIDE + LEGAL SERVICES. DISTRIBUTION OF THIS DOCUMENT DOES NOT CREATE AN + ATTORNEY-CLIENT RELATIONSHIP. CREATIVE COMMONS PROVIDES THIS + INFORMATION ON AN "AS-IS" BASIS. CREATIVE COMMONS MAKES NO WARRANTIES + REGARDING THE USE OF THIS DOCUMENT OR THE INFORMATION OR WORKS + PROVIDED HEREUNDER, AND DISCLAIMS LIABILITY FOR DAMAGES RESULTING FROM + THE USE OF THIS DOCUMENT OR THE INFORMATION OR WORKS PROVIDED + HEREUNDER. + +Statement of Purpose + +The laws of most jurisdictions throughout the world automatically confer +exclusive Copyright and Related Rights (defined below) upon the creator +and subsequent owner(s) (each and all, an "owner") of an original work of +authorship and/or a database (each, a "Work"). + +Certain owners wish to permanently relinquish those rights to a Work for +the purpose of contributing to a commons of creative, cultural and +scientific works ("Commons") that the public can reliably and without fear +of later claims of infringement build upon, modify, incorporate in other +works, reuse and redistribute as freely as possible in any form whatsoever +and for any purposes, including without limitation commercial purposes. +These owners may contribute to the Commons to promote the ideal of a free +culture and the further production of creative, cultural and scientific +works, or to gain reputation or greater distribution for their Work in +part through the use and efforts of others. + +For these and/or other purposes and motivations, and without any +expectation of additional consideration or compensation, the person +associating CC0 with a Work (the "Affirmer"), to the extent that he or she +is an owner of Copyright and Related Rights in the Work, voluntarily +elects to apply CC0 to the Work and publicly distribute the Work under its +terms, with knowledge of his or her Copyright and Related Rights in the +Work and the meaning and intended legal effect of CC0 on those rights. + +1. Copyright and Related Rights. A Work made available under CC0 may be +protected by copyright and related or neighboring rights ("Copyright and +Related Rights"). Copyright and Related Rights include, but are not +limited to, the following: + + i. the right to reproduce, adapt, distribute, perform, display, + communicate, and translate a Work; + ii. moral rights retained by the original author(s) and/or performer(s); +iii. publicity and privacy rights pertaining to a person's image or + likeness depicted in a Work; + iv. rights protecting against unfair competition in regards to a Work, + subject to the limitations in paragraph 4(a), below; + v. rights protecting the extraction, dissemination, use and reuse of data + in a Work; + vi. database rights (such as those arising under Directive 96/9/EC of the + European Parliament and of the Council of 11 March 1996 on the legal + protection of databases, and under any national implementation + thereof, including any amended or successor version of such + directive); and +vii. other similar, equivalent or corresponding rights throughout the + world based on applicable law or treaty, and any national + implementations thereof. + +2. Waiver. To the greatest extent permitted by, but not in contravention +of, applicable law, Affirmer hereby overtly, fully, permanently, +irrevocably and unconditionally waives, abandons, and surrenders all of +Affirmer's Copyright and Related Rights and associated claims and causes +of action, whether now known or unknown (including existing as well as +future claims and causes of action), in the Work (i) in all territories +worldwide, (ii) for the maximum duration provided by applicable law or +treaty (including future time extensions), (iii) in any current or future +medium and for any number of copies, and (iv) for any purpose whatsoever, +including without limitation commercial, advertising or promotional +purposes (the "Waiver"). Affirmer makes the Waiver for the benefit of each +member of the public at large and to the detriment of Affirmer's heirs and +successors, fully intending that such Waiver shall not be subject to +revocation, rescission, cancellation, termination, or any other legal or +equitable action to disrupt the quiet enjoyment of the Work by the public +as contemplated by Affirmer's express Statement of Purpose. + +3. Public License Fallback. Should any part of the Waiver for any reason +be judged legally invalid or ineffective under applicable law, then the +Waiver shall be preserved to the maximum extent permitted taking into +account Affirmer's express Statement of Purpose. In addition, to the +extent the Waiver is so judged Affirmer hereby grants to each affected +person a royalty-free, non transferable, non sublicensable, non exclusive, +irrevocable and unconditional license to exercise Affirmer's Copyright and +Related Rights in the Work (i) in all territories worldwide, (ii) for the +maximum duration provided by applicable law or treaty (including future +time extensions), (iii) in any current or future medium and for any number +of copies, and (iv) for any purpose whatsoever, including without +limitation commercial, advertising or promotional purposes (the +"License"). The License shall be deemed effective as of the date CC0 was +applied by Affirmer to the Work. Should any part of the License for any +reason be judged legally invalid or ineffective under applicable law, such +partial invalidity or ineffectiveness shall not invalidate the remainder +of the License, and in such case Affirmer hereby affirms that he or she +will not (i) exercise any of his or her remaining Copyright and Related +Rights in the Work or (ii) assert any associated claims and causes of +action with respect to the Work, in either case contrary to Affirmer's +express Statement of Purpose. + +4. Limitations and Disclaimers. + + a. No trademark or patent rights held by Affirmer are waived, abandoned, + surrendered, licensed or otherwise affected by this document. + b. Affirmer offers the Work as-is and makes no representations or + warranties of any kind concerning the Work, express, implied, + statutory or otherwise, including without limitation warranties of + title, merchantability, fitness for a particular purpose, non + infringement, or the absence of latent or other defects, accuracy, or + the present or absence of errors, whether or not discoverable, all to + the greatest extent permissible under applicable law. + c. Affirmer disclaims responsibility for clearing rights of other persons + that may apply to the Work or any use thereof, including without + limitation any person's Copyright and Related Rights in the Work. + Further, Affirmer disclaims responsibility for obtaining any necessary + consents, permissions or other rights required for any use of the + Work. + d. Affirmer understands and acknowledges that Creative Commons is not a + party to this document and has no duty or obligation with respect to + this CC0 or use of the Work. 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) - -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 -- <> de i1 plus grand - return greater; - elsif i1.blk(i) < i2.blk(i) then -- <> 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 <> - 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 <>: 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; -- cgit