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 --- src/rationals.adb | 482 ++++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 339 insertions(+), 143 deletions(-) (limited to 'src/rationals.adb') 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; -- cgit