File: [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Numbers / multprec_natural_numbers.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:26 2000 UTC (23 years, 7 months ago) by maekawa
Branch: PHC, MAIN
CVS Tags: v2, maekawa-ipv6, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, HEAD
Changes since 1.1: +0 -0
lines
Import the second public release of PHCpack.
OKed by Jan Verschelde.
|
with unchecked_deallocation;
package body Multprec_Natural_Numbers is
-- NOTES ON THE CHOICE OF REPRESENTATION AND IMPLEMENTATION :
-- 1) The base is decimal, for convenience of input/output.
-- Everything should remain correct if another base is chosen,
-- except of course the function "Decimal_Places".
-- The base is such that factor times the base is still a standard
-- natural number, so that standard arithmetic can be used.
-- 2) A natural number is implemented by means of a coefficient array.
-- The advantage over a linked-list representation is that the
-- coefficients can be accessed faster. The disadvantage of some
-- wasted storage is countered by the aim of implementing higher
-- precision floating-point numbers, where the mantissa is fixed.
-- 3) The usage of a pointer to access the coefficient vector permits
-- to have uniformity in input/output, since on input it is not
-- known in advance how long a natural number may be.
-- 4) A natural number n for which Empty(n) holds is considered as zero.
fact : constant natural := 10;
expo : constant natural := 8;
the_base : constant natural := fact**expo;
-- DATATRUCTURES :
type Natural_Number_Rep ( size : natural ) is record
numb : Array_of_Naturals(0..size);
end record;
-- AUXILIARY :
function Size ( n : natural ) return natural is
-- DESCRIPTION :
-- Determines the size of the natural number to represent n.
acc : natural := the_base;
begin
for i in 0..n loop
if n < acc
then return i;
else acc := acc*the_base;
end if;
end loop;
return n;
end Size;
-- CREATORS :
function Create ( n : natural ) return Array_of_Naturals is
res : Array_of_Naturals(0..Size(n));
remainder : natural := n;
begin
for i in res'range loop
res(i) := remainder mod the_base;
remainder := remainder/the_base;
end loop;
return res;
end Create;
function Create ( n : natural ) return Natural_Number is
res : Natural_Number := Create(Create(n));
begin
return res;
end Create;
function Create ( n : Array_of_Naturals ) return Natural_Number is
res : Natural_Number;
res_rep : Natural_Number_Rep(n'last);
begin
res_rep.numb := n;
res := new Natural_Number_Rep'(res_rep);
return res;
end Create;
function Create ( n : Natural_Number ) return natural is
res : natural;
begin
if Empty(n)
then res := 0;
else res := Coefficient(n,Size(n));
for i in reverse 0..Size(n)-1 loop
res := res*the_base + Coefficient(n,i);
end loop;
end if;
return res;
end Create;
-- SELECTORS :
function Base return natural is
begin
return the_base;
end Base;
function Exponent return natural is
begin
return expo;
end Exponent;
function Empty ( n : Natural_Number ) return boolean is
begin
return ( n = null );
end Empty;
function Size ( n : Natural_Number ) return natural is
begin
if n = null
then return 0;
else return n.size;
end if;
end Size;
function Decimal_Places ( n : natural ) return natural is
acc : natural := 1;
begin
for i in 0..n loop
if n < acc
then return i;
else acc := acc*10;
end if;
end loop;
return n;
end Decimal_Places;
function Decimal_Places ( n : Natural_Number ) return natural is
res : natural := 0;
begin
if not Empty(n)
then for i in reverse 0..Size(n) loop
if n.numb(i) /= 0
then res := Decimal_Places(n.numb(i)) + i*expo;
exit;
end if;
end loop;
end if;
return res;
end Decimal_Places;
function Coefficient ( n : Natural_Number; i : natural ) return natural is
begin
if (n = null or else i > n.size)
then return 0;
else return n.numb(i);
end if;
end Coefficient;
function Coefficients ( n : Natural_Number ) return Array_of_Naturals is
na : Array_of_Naturals(0..0) := (0..0 => 0);
begin
if Empty(n)
then return na;
else return n.numb;
end if;
end Coefficients;
-- COMPARISON AND COPYING :
-- Note that these implementations are all independent of the
-- actual representation of Natural_Number as they use the selectors.
function Equal ( n1,n2 : Array_of_Naturals ) return boolean is
begin
if n1'first /= n2'first
then return false;
elsif n1'last /= n2'last
then return false;
else for i in n1'range loop
if n1(i) /= n2(i)
then return false;
end if;
end loop;
return true;
end if;
end Equal;
function Equal ( n1 : Natural_Number; n2 : natural ) return boolean is
begin
if Empty(n1)
then if n2 = 0
then return true;
else return false;
end if;
else declare
n2cff : constant Array_of_Naturals := Create(n2);
begin
if n2cff'last > Size(n1)
then return false;
else for i in n2cff'range loop
if Coefficient(n1,i) /= n2cff(i)
then return false;
end if;
end loop;
for i in n2cff'last+1..Size(n1) loop
if Coefficient(n1,i) /= 0
then return false;
end if;
end loop;
end if;
end;
return true;
end if;
end Equal;
function Equal ( n1,n2 : Natural_Number ) return boolean is
min_size : natural;
begin
if Empty(n1)
then return Equal(n2,0);
else if Empty(n2)
then return Equal(n1,0);
else if Size(n1) < Size(n2)
then for i in Size(n1)+1..Size(n2) loop
if Coefficient(n2,i) /= 0
then return false;
end if;
end loop;
min_size := Size(n1);
elsif Size(n1) > Size(n2)
then for i in Size(n2)+1..Size(n1) loop
if Coefficient(n1,i) /= 0
then return false;
end if;
end loop;
min_size := Size(n2);
else min_size := Size(n1);
end if;
return Equal(Coefficients(n1)(0..min_size),
Coefficients(n2)(0..min_size));
end if;
end if;
end Equal;
function "<" ( n1 : Natural_Number; n2 : natural ) return boolean is
begin
if Empty(n1)
then if n2 > 0
then return true;
else return false;
end if;
else declare
n2cff : constant Array_of_Naturals := Create(n2);
begin
if n2cff'last > Size(n1)
then return true;
else if n2cff'last < Size(n1)
then for i in reverse n2cff'last+1..Size(n1) loop
if Coefficient(n1,i) /= 0
then return false;
end if;
end loop;
end if;
for i in reverse n2cff'range loop
if Coefficient(n1,i) >= n2cff(i)
then return false;
end if;
end loop;
end if;
end;
return true;
end if;
end "<";
function "<" ( n1 : natural; n2 : Natural_Number ) return boolean is
begin
if Empty(n2)
then return false;
else declare
n1cff : constant Array_of_Naturals := Create(n1);
begin
if n1cff'last > Size(n2)
then return false;
else if n1cff'last < Size(n2)
then for i in reverse n1cff'last+1..Size(n2) loop
if Coefficient(n2,i) /= 0
then return true;
end if;
end loop;
end if;
for i in reverse n1cff'range loop
if n1cff(i) >= Coefficient(n2,i)
then return false;
end if;
end loop;
end if;
end;
return true;
end if;
end "<";
function "<" ( n1,n2 : Natural_Number ) return boolean is
min_size : natural;
begin
if Empty(n1)
then if Empty(n2)
then return false;
else return true;
end if;
else if Empty(n2)
then return false;
else if Size(n1) < Size(n2)
then for i in Size(n1)+1..Size(n2) loop
if Coefficient(n2,i) /= 0
then return true;
end if;
end loop;
min_size := Size(n1);
elsif Size(n1) > Size(n2)
then for i in Size(n2)+1..Size(n1) loop
if Coefficient(n1,i) /= 0
then return false;
end if;
end loop;
min_size := Size(n2);
else min_size := Size(n1);
end if;
for i in reverse 0..min_size loop
if Coefficient(n1,i) < Coefficient(n2,i)
then return true;
elsif Coefficient(n1,i) > Coefficient(n2,i)
then return false;
end if;
end loop;
return false; -- n1 = n2
end if;
end if;
end "<";
function ">" ( n1 : Natural_Number; n2 : natural ) return boolean is
begin
if Empty(n1)
then return false;
else declare
n2cff : constant Array_of_Naturals := Create(n2);
begin
if n2cff'last > Size(n1)
then return false;
else if n2cff'last < Size(n1)
then for i in n2cff'last+1..Size(n1) loop
if Coefficient(n1,i) /= 0
then return true;
end if;
end loop;
end if;
for i in reverse n2cff'range loop
if Coefficient(n1,i) <= n2cff(i)
then return false;
end if;
end loop;
end if;
end;
return true;
end if;
end ">";
function ">" ( n1 : natural; n2 : Natural_Number ) return boolean is
begin
if Empty(n2)
then if n1 > 0
then return true;
else return false;
end if;
else declare
n1cff : constant Array_of_Naturals := Create(n1);
begin
if n1cff'last > Size(n2)
then return true;
else if n1cff'last < Size(n2)
then for i in n1cff'last+1..Size(n2) loop
if Coefficient(n2,i) /= 0
then return false;
end if;
end loop;
end if;
for i in reverse n1cff'range loop
if n1cff(i) <= Coefficient(n2,i)
then return false;
end if;
end loop;
end if;
end;
return true;
end if;
end ">";
function ">" ( n1,n2 : Natural_Number ) return boolean is
min_size : natural;
begin
if Empty(n2)
then if Empty(n1)
then return false;
else return true;
end if;
else if Empty(n1)
then return false;
else if Size(n1) < Size(n2)
then for i in Size(n1)+1..Size(n2) loop
if Coefficient(n2,i) /= 0
then return false;
end if;
end loop;
min_size := Size(n1);
elsif Size(n1) > Size(n2)
then for i in Size(n2)+1..Size(n1) loop
if Coefficient(n1,i) /= 0
then return true;
end if;
end loop;
min_size := Size(n2);
else min_size := Size(n1);
end if;
for i in reverse 0..min_size loop
if Coefficient(n1,i) > Coefficient(n2,i)
then return true;
elsif Coefficient(n1,i) < Coefficient(n2,i)
then return false;
end if;
end loop;
return false; -- n1 = n2
end if;
end if;
end ">";
procedure Copy ( n1 : in natural; n2 : in out Natural_Number ) is
begin
Clear(n2);
n2 := Create(n1);
end Copy;
procedure Copy ( n1 : in Natural_Number; n2 : in out Natural_Number ) is
begin
Clear(n2);
if not Empty(n1)
then declare
n2rep : Natural_Number_Rep(Size(n1));
begin
for i in n2rep.numb'range loop
n2rep.numb(i) := n1.numb(i);
end loop;
n2 := new Natural_Number_Rep'(n2rep);
end;
end if;
end Copy;
-- AUXILIARIES FOR ARITHMETIC OPERATIONS :
function Extend ( n : Array_of_Naturals; carry : natural )
return Natural_Number_Rep is
-- DESCRIPTION :
-- Extends the array of naturals as much as need to store the carry
-- obtained from adding a number to a natural number.
begin
if carry = 0
then declare
nrep : Natural_Number_Rep(n'last);
begin
nrep.numb := n;
return nrep;
end;
else declare
np1 : Array_of_Naturals(n'first..n'last+1);
carry1 : natural;
begin
np1(n'range) := n;
np1(np1'last) := carry mod the_base;
carry1 := carry/the_base;
return Extend(np1,carry1);
end;
end if;
end Extend;
procedure Extend ( n : in out Natural_Number; carry : in natural ) is
-- DESCRIPTION :
-- Extends the coefficient vector of the natural number to store
-- the carry-over.
begin
if carry /= 0
then declare
res : Natural_Number
:= new Natural_Number_Rep'(Extend(n.numb,carry));
begin
Clear(n);
n := res;
end;
end if;
end Extend;
procedure Mul_Fact ( n : in out Natural_Number; f : in natural ) is
-- DESCRIPTION :
-- Multiplies the number n with 0 < f <= fact and stores the result in n.
prod,carry : natural;
begin
if not Empty(n)
then carry := 0;
for i in n.numb'range loop
prod := n.numb(i)*f + carry;
n.numb(i) := prod mod the_base;
carry := prod/the_base;
end loop;
Extend(n,carry);
end if;
end Mul_Fact;
procedure Mul_Base ( n : in out Natural_Number; k : in natural ) is
-- DESCRIPTION :
-- Multiplies the number n with the_base**k.
begin
if k > 0 and not Empty(n)
then declare
npk : Natural_Number_Rep(n.size+k);
begin
npk.numb(0..k-1) := (0..k-1 => 0);
npk.numb(k..npk.numb'last) := n.numb(n.numb'range);
Clear(n);
n := new Natural_Number_Rep'(npk);
end;
end if;
end Mul_Base;
procedure Div_Base ( n : in out Natural_Number; k : in natural ) is
-- DESCRIPTION :
-- Divides the number n by the_base**k.
begin
if k > 0 and not Empty(n)
then if k > Size(n)
then Clear(n);
else declare
nmk : Natural_Number_Rep(n.size-k);
begin
for i in nmk.numb'range loop
nmk.numb(i) := n.numb(k+i);
end loop;
Clear(n);
n := new Natural_Number_Rep'(nmk);
end;
end if;
end if;
end Div_Base;
-- ARITHMETIC OPERATIONS :
function "+" ( n1 : Natural_Number; n2 : natural ) return Natural_Number is
res : Natural_Number;
begin
if Empty(n1)
then res := Create(n2);
else declare
cff : Array_of_Naturals(n1.numb'range);
sum,carry : natural;
begin
carry := n2;
for i in cff'range loop
sum := n1.numb(i) + carry;
cff(i) := sum mod the_base;
carry := sum/the_base;
end loop;
res := new Natural_Number_Rep'(Extend(cff,carry));
end;
end if;
return res;
end "+";
function "+" ( n1 : natural; n2 : Natural_Number ) return Natural_Number is
begin
return n2+n1;
end "+";
function "+" ( n1,n2 : Natural_Number ) return Natural_Number is
res : Natural_Number;
sum,carry : natural;
begin
if Empty(n1)
then Copy(n2,res);
elsif Empty(n2)
then Copy(n1,res);
else carry := 0;
if Size(n1) > Size(n2)
then declare
res_rep : Natural_Number_Rep(Size(n1)) := n1.all;
begin
for i in n2.numb'range loop
sum := res_rep.numb(i) + n2.numb(i) + carry;
res_rep.numb(i) := sum mod the_base;
carry := sum/the_base;
end loop;
if carry /= 0
then for i in n2.numb'last+1..n1.numb'last loop
sum := res_rep.numb(i) + carry;
res_rep.numb(i) := sum mod the_base;
carry := sum/the_base;
exit when (carry = 0);
end loop;
end if;
res := new Natural_Number_Rep'(res_rep);
end;
else declare
res_rep : Natural_Number_Rep(Size(n2)) := n2.all;
begin
for i in n1.numb'range loop
sum := res_rep.numb(i) + n1.numb(i) + carry;
res_rep.numb(i) := sum mod the_base;
carry := sum/the_base;
end loop;
if carry /= 0
then for i in n1.numb'last+1..n2.numb'last loop
sum := res_rep.numb(i) + carry;
res_rep.numb(i) := sum mod the_base;
carry := sum/the_base;
exit when (carry = 0);
end loop;
end if;
res := new Natural_Number_Rep'(res_rep);
end;
end if;
Extend(res,carry);
end if;
return res;
end "+";
procedure Add ( n1 : in out Natural_Number; n2 : in natural ) is
sum,carry : natural;
begin
if Empty(n1)
then n1 := Create(n2);
else carry := n2;
for i in n1.numb'range loop
sum := n1.numb(i) + carry;
n1.numb(i) := sum mod the_base;
carry := sum/the_base;
end loop;
Extend(n1,carry);
end if;
end Add;
procedure Add ( n1 : in out Natural_Number; n2 : in Natural_Number ) is
res : Natural_Number;
sum,carry : natural;
begin
if Empty(n1)
then Copy(n2,n1);
elsif not Empty(n2)
then carry := 0;
if Size(n1) >= Size(n2)
then for i in n2.numb'range loop
sum := n1.numb(i) + n2.numb(i) + carry;
n1.numb(i) := sum mod the_base;
carry := sum/the_base;
end loop;
if carry /= 0
then for i in n2.numb'last+1..n1.numb'last loop
sum := n1.numb(i) + carry;
n1.numb(i) := sum mod the_base;
carry := sum/the_base;
exit when (carry = 0);
end loop;
Extend(n1,carry);
end if;
else declare
res_rep : Natural_Number_Rep(Size(n2)) := n2.all;
begin
for i in n1.numb'range loop
sum := res_rep.numb(i) + n1.numb(i) + carry;
res_rep.numb(i) := sum mod the_base;
carry := sum/the_base;
end loop;
if carry /= 0
then for i in n1.numb'last+1..n2.numb'last loop
sum := res_rep.numb(i) + carry;
res_rep.numb(i) := sum mod the_base;
carry := sum/the_base;
exit when (carry = 0);
end loop;
end if;
Clear(n1);
n1 := new Natural_Number_Rep'(res_rep);
if carry /= 0
then Extend(n1,carry);
end if;
end;
end if;
end if;
end Add;
function "-" ( n1 : Natural_Number; n2 : natural ) return Natural_Number is
res : Natural_Number;
diff,carry : integer;
begin
if not Empty(n1)
then Copy(n1,res);
declare
n2cff : constant Array_of_Naturals := Create(n2);
index : natural := n2cff'first;
begin
carry := n2cff(index);
for i in n1.numb'range loop
diff := n1.numb(i) - carry;
if diff >= 0
then res.numb(i) := diff mod the_base;
carry := diff/the_base;
else diff := the_base + diff;
res.numb(i) := diff mod the_base;
carry := diff/the_base + (res.numb(i)+carry)/the_base;
end if;
if index < n2cff'last
then index := index+1;
carry := carry + n2cff(index);
end if;
exit when (carry = 0);
end loop;
end;
end if;
return res;
end "-";
function "-" ( n1 : natural; n2 : Natural_Number ) return Natural_Number is
tmp : Natural_Number := Create(n1);
res : Natural_Number := tmp - n2;
begin
Clear(tmp);
return res;
end "-";
function "-" ( n1,n2 : Natural_Number ) return Natural_Number is
res,diff : Natural_Number;
begin
Copy(n1,res);
if not Empty(n2)
then for i in reverse n2.numb'range loop
Copy(res,diff);
Div_Base(diff,i);
if not Empty(diff)
then Sub(diff,n2.numb(i));
for j in diff.numb'range loop
res.numb(i+j) := diff.numb(j);
end loop;
Clear(diff);
end if;
end loop;
end if;
return res;
end "-";
function "+" ( n : Natural_Number ) return Natural_Number is
res : Natural_Number;
begin
Copy(n,res);
return res;
end "+";
function "-" ( n : Natural_Number ) return Natural_Number is
res : Natural_Number;
begin
Copy(n,res);
return res;
end "-";
procedure Min ( n : in out Natural_Number ) is
begin
null;
end Min;
procedure Sub ( n1 : in out Natural_Number; n2 : in natural ) is
diff,carry : integer;
begin
if not Empty(n1) and n2 > 0
then declare
n2cff : constant Array_of_Naturals := Create(n2);
index : natural := n2cff'first;
begin
carry := n2cff(index);
for i in n1.numb'range loop
diff := n1.numb(i) - carry;
if diff >= 0
then n1.numb(i) := diff mod the_base;
carry := diff/the_base;
else diff := the_base + diff;
n1.numb(i) := diff mod the_base;
carry := diff/the_base + (n1.numb(i)+carry)/the_base;
end if;
if index < n2cff'last
then index := index+1;
carry := carry + n2cff(index);
end if;
exit when (carry = 0);
end loop;
end;
end if;
end Sub;
procedure Sub ( n1 : in out Natural_Number; n2 : in Natural_Number ) is
res,diff : Natural_Number;
begin
if not Empty(n1) and not Empty(n2)
then for i in reverse n2.numb'range loop
Copy(n1,diff);
Div_Base(diff,i);
if not Empty(diff)
then Sub(diff,n2.numb(i));
for j in diff.numb'range loop
n1.numb(i+j) := diff.numb(j);
end loop;
Clear(diff);
end if;
end loop;
end if;
end Sub;
function "*" ( n1 : Natural_Number; n2 : natural ) return Natural_Number is
res,prod,nres : Natural_Number;
remainder,modfact : natural;
cnt : natural := 0;
begin
if n2 /= 0
then remainder := n2;
while remainder /= 0 loop
modfact := remainder mod fact;
if modfact /= 0
then Copy(n1,prod);
Mul_Fact(prod,modfact);
for i in 1..cnt loop
Mul_Fact(prod,fact);
end loop;
Add(res,prod);
Clear(prod);
end if;
cnt := cnt+1;
remainder := remainder/fact;
end loop;
end if;
return res;
end "*";
function "*" ( n1 : natural; n2 : Natural_Number ) return Natural_Number is
begin
return n2*n1;
end "*";
function "*" ( n1,n2 : Natural_Number ) return Natural_Number is
res,prod : Natural_Number;
begin
if (Empty(n1) or else Empty(n2))
then return res;
else if Size(n1) >= Size(n2)
then for i in n2.numb'range loop
prod := n1*n2.numb(i);
Mul_Base(prod,i);
Add(res,prod);
Clear(prod);
end loop;
else for i in n1.numb'range loop
prod := n2*n1.numb(i);
Mul_Base(prod,i);
Add(res,prod);
Clear(prod);
end loop;
end if;
end if;
return res;
end "*";
procedure Mul ( n1 : in out Natural_Number; n2 : in natural ) is
res,prod : Natural_Number;
remainder,modfact : natural;
cnt : natural := 0;
begin
if n2 = 0
then Clear(n1);
else remainder := n2;
while remainder /= 0 loop
modfact := remainder mod fact;
if modfact /= 0
then Copy(n1,prod);
Mul_Fact(prod,modfact);
for i in 1..cnt loop
Mul_Fact(prod,fact);
end loop;
Add(res,prod);
Clear(prod);
end if;
cnt := cnt+1;
remainder := remainder/fact;
end loop;
Clear(n1);
n1 := res;
end if;
end Mul;
procedure Mul ( n1 : in out Natural_Number; n2 : in Natural_Number ) is
res,prod : Natural_Number;
begin
if not Empty(n1)
then if Empty(n2)
then Clear(n1);
else if Size(n1) >= Size(n2)
then for i in n2.numb'range loop
prod := n1*n2.numb(i);
Mul_Base(prod,i);
Add(res,prod);
Clear(prod);
end loop;
else for i in n1.numb'range loop
prod := n2*n1.numb(i);
Mul_Base(prod,i);
Add(res,prod);
Clear(prod);
end loop;
end if;
Clear(n1);
n1 := res;
end if;
end if;
end Mul;
function "**" ( n1 : Natural_Number; n2 : natural ) return Natural_Number is
res : Natural_Number;
begin
if n2 = 0
then res := Create(1);
else if not Empty(n1)
then Copy(n1,res);
for i in 1..n2-1 loop
Mul(res,n1);
end loop;
end if;
end if;
return res;
end "**";
function "**" ( n1 : natural; n2 : Natural_Number ) return Natural_Number is
res,cnt : Natural_Number;
begin
if Equal(n2,0)
then res := Create(1);
else res := Create(n1);
if n1 /= 0
then cnt := Create(1);
while not Equal(cnt,n2) loop
Mul(res,n1);
Add(cnt,1);
end loop;
Clear(cnt);
end if;
end if;
return res;
end "**";
function "**" ( n1,n2 : Natural_Number ) return Natural_Number is
res,cnt : Natural_Number;
begin
if Equal(n2,0)
then res := Create(1);
else if not Empty(n1)
then Copy(n1,res);
cnt := Create(1);
while not Equal(cnt,n2) loop
Mul(res,n1);
Add(cnt,1);
end loop;
Clear(cnt);
end if;
end if;
return res;
end "**";
procedure Small_Div ( n1 : in Natural_Number; n2 : in natural;
q : out Natural_Number; r : out natural ) is
-- DESCRIPTION :
-- n1 = n2*q + r, only applies when n2 <= fact.
-- res : Natural_Number;
carry : natural := 0;
begin
if n2 /= 0
then if not Empty(n1)
then -- res := new Natural_Number_Rep(Size(n1));
q := new Natural_Number_Rep(Size(n1));
for i in reverse 1..Size(n1) loop
-- res.numb(i) := (n1.numb(i)+carry)/n2;
q.numb(i) := (n1.numb(i)+carry)/n2;
carry := (n1.numb(i)+carry) mod n2;
carry := carry*the_base;
end loop;
-- res.numb(0) := (n1.numb(0)+carry)/n2;
q.numb(0) := (n1.numb(0)+carry)/n2;
carry := (n1.numb(0)+carry) mod n2;
end if;
else raise NUMERIC_ERROR;
end if;
-- q := res;
r := carry;
end Small_Div;
procedure Small_Div ( n1 : in out Natural_Number; n2 : in natural;
r : out natural ) is
q : Natural_Number;
begin
Small_Div(n1,n2,q,r);
Copy(q,n1); Clear(q);
end Small_Div;
procedure Big_Div ( n1 : in Natural_Number; n2 : in natural;
q : out Natural_Number; r : out natural ) is
-- DESCRIPTION :
-- This division has to be applied when n2 > fact.
wrkn1,shiftn1,remn1,acc,qres,qsub : Natural_Number;
cntshf,cntsub,rres : natural;
begin
if n2 /= 0
then
if n1 < n2
then
q := qres;
r := Create(n1);
else
Copy(n1,shiftn1); -- we have n1 >= n2
Small_Div(shiftn1,fact,wrkn1,rres);
cntshf := 0;
while not (wrkn1 < n2) loop -- invariant n2 <= shiftn1
acc := Create(rres);
for i in 1..cntshf loop
Mul_Fact(acc,fact);
end loop;
Add(remn1,acc); Clear(acc);
Copy(wrkn1,shiftn1);
cntshf := cntshf + 1;
Small_Div(wrkn1,fact,rres); -- at end of loop
end loop; -- shiftn1/fact < n2 <= shiftn1
Clear(wrkn1);
cntsub := 0;
while not (shiftn1 < n2) loop -- subtract n2 from shiftn1
cntsub := cntsub + 1;
Sub(shiftn1,n2);
end loop;
qsub := Create(cntsub); -- intermediate quotient
for i in 1..cntshf loop
Mul_Fact(qsub,fact);
Mul_Fact(shiftn1,fact);
end loop;
Add(shiftn1,remn1);
Clear(remn1);
Div(shiftn1,n2,qres,r); -- recursive call
Add(qsub,qres); -- collect results
q := qsub;
Clear(qres); Clear(shiftn1); -- cleaning up
end if;
else
raise NUMERIC_ERROR;
end if;
end Big_Div;
function "/" ( n1 : Natural_Number; n2 : natural ) return Natural_Number is
res : Natural_Number;
r : natural;
begin
if not Empty(n1)
then if n2 <= fact
then Small_Div(n1,n2,res,r);
else Big_Div(n1,n2,res,r);
end if;
end if;
return res;
end "/";
function "/" ( n1 : natural; n2 : Natural_Number ) return natural is
res : natural;
begin
if n1 < n2
then res := 0;
elsif not Empty(n2)
then res := Create(n2);
res := n1/res;
else raise NUMERIC_ERROR;
end if;
return res;
end "/";
function "/" ( n1,n2 : Natural_Number ) return Natural_Number is
res,rrr : Natural_Number;
begin
Div(n1,n2,res,rrr);
Clear(rrr);
return res;
end "/";
function Rmd ( n1 : Natural_Number; n2 : natural ) return natural is
res : natural;
q : Natural_Number;
begin
if n2 <= fact
then Small_Div(n1,n2,q,res);
else Big_Div(n1,n2,q,res);
end if;
Clear(q);
return res;
end Rmd;
function Rmd ( n1 : natural; n2 : Natural_Number ) return natural is
res : natural;
begin
if n1 < n2
then res := n1;
elsif not Empty(n2)
then res := Create(n2);
res := n1 mod res;
else raise NUMERIC_ERROR;
end if;
return res;
end Rmd;
function Rmd ( n1,n2 : Natural_Number ) return Natural_Number is
res,q : Natural_Number;
begin
Div(n1,n2,q,res);
Clear(q);
return res;
end Rmd;
procedure Div ( n1 : in out Natural_Number; n2 : in natural ) is
r : natural;
begin
Div(n1,n2,r);
end Div;
procedure Div ( n1 : in out Natural_Number; n2 : in Natural_Number ) is
r : Natural_Number;
begin
Div(n1,n2,r);
Clear(r);
end Div;
procedure Div ( n1 : in Natural_Number; n2 : in natural;
q : out Natural_Number; r : out natural ) is
begin
if n2 <= fact
then Small_Div(n1,n2,q,r);
else Big_Div(n1,n2,q,r);
end if;
end Div;
procedure Div ( n1 : in out Natural_Number; n2 : in natural;
r : out natural ) is
q : Natural_Number;
begin
Div(n1,n2,q,r);
Copy(q,n1); Clear(q);
end Div;
procedure Div ( n1,n2 : in Natural_Number;
q : out Natural_Number; r : out Natural_Number ) is
wrkn1,shiftn1,remn1,acc,qres,qsub : Natural_Number;
cntshf,cntsub,rres : natural;
begin
if (Empty(n2) or else Equal(n2,0))
then
raise NUMERIC_ERROR;
else
if n1 < n2
then
q := qres;
Copy(n1,r);
else
Copy(n1,shiftn1); -- we have n1 >= n2
Small_Div(shiftn1,fact,wrkn1,rres);
cntshf := 0;
while not (wrkn1 < n2) loop -- invariant n2 <= shiftn1
acc := Create(rres);
for i in 1..cntshf loop
Mul_Fact(acc,fact);
end loop;
Add(remn1,acc); Clear(acc);
Copy(wrkn1,shiftn1);
cntshf := cntshf + 1;
Small_Div(wrkn1,fact,rres); -- at end of loop
end loop; -- shiftn1/fact < n2 <= shiftn1
Clear(wrkn1);
cntsub := 0;
while not (shiftn1 < n2) loop -- subtract n2 from shiftn1
cntsub := cntsub + 1;
Sub(shiftn1,n2);
end loop;
qsub := Create(cntsub); -- intermediate quotient
for i in 1..cntshf loop
Mul_Fact(qsub,fact);
Mul_Fact(shiftn1,fact);
end loop;
Add(shiftn1,remn1);
Clear(remn1);
Div(shiftn1,n2,qres,r); -- recursive call
Add(qsub,qres); -- collect results
q := qsub;
Clear(qres); Clear(shiftn1); -- cleaning up
end if;
end if;
end Div;
procedure Div ( n1 : in out Natural_Number; n2 : in Natural_Number;
r : out Natural_Number ) is
q : Natural_Number;
begin
Div(n1,n2,q,r);
Copy(q,n1); Clear(q);
end Div;
-- DESTRUCTOR :
procedure Clear ( n : in out Natural_Number ) is
procedure free is
new unchecked_deallocation(Natural_Number_Rep,Natural_Number);
begin
if not Empty(n)
then free(n);
end if;
end Clear;
end Multprec_Natural_Numbers;
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to multprec_natural_numbers.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Numbers