From 2b8b55de4a18757e8d6769e458c84f7c1df1e261 Mon Sep 17 00:00:00 2001 From: Jed Barber Date: Mon, 13 Feb 2017 18:27:13 +1100 Subject: Swapped out crypto package for something smaller, revised other code and readme/notes slightly --- src/multi_precision_integers.adb | 1500 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 1500 insertions(+) create mode 100644 src/multi_precision_integers.adb (limited to 'src/multi_precision_integers.adb') diff --git a/src/multi_precision_integers.adb b/src/multi_precision_integers.adb new file mode 100644 index 0000000..4d8d876 --- /dev/null +++ b/src/multi_precision_integers.adb @@ -0,0 +1,1500 @@ +----------------------------------------------------------------------------- +-- 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: !! + +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; -- cgit