File: [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices / greatest_common_divisors.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:23 2000 UTC (23 years, 8 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.
|
package body Greatest_Common_Divisors is
use Euclidean_Domain;
function pos_gcd ( a,b : number ) return number is
res,r : number;
begin
if Equal(b,zero)
then Copy(a,res);
else r := Rmd(a,b);
if Equal(r,zero)
then Copy(b,res);
else res := pos_gcd(b,r);
end if;
Clear(r);
end if;
return res;
end pos_gcd;
function gcd ( a,b : number ) return number is
res,mina,minb : number;
begin
if a < zero
then if b < zero
then mina := -a; minb := -b;
res := pos_gcd(mina,minb);
Clear(mina); Clear(minb);
else mina := -a;
res := pos_gcd(mina,b);
Clear(mina);
end if;
elsif b < zero
then minb := -b;
res := pos_gcd(a,minb);
Clear(minb);
else res := pos_gcd(a,b);
end if;
return res;
end gcd;
function lcm ( a,b : number ) return number is
res : number;
g : number := gcd(a,b);
begin
if Equal(g,zero)
then res := zero;
else res := a/g;
Mul(res,b);
end if;
Clear(g);
return res;
end lcm;
procedure gcd ( a,b : in number; k,l,d : out number ) is
mina,minb,kk,ll : number;
procedure pos_gcd ( a,b : in number; k,l,d : out number ) is
r,aa,bb,k1,l1,k0,l0,h,tmp : number;
-- nk1,nl1 : number;
-- REQUIRED :
-- a >= b and a > 0 and b > 0
begin
r := Rmd(a,b);
if Equal(r,zero)
then k := zero; l := one; Copy(b,d);
else k0 := zero; l0 := one;
Copy(r,bb);
Copy(b,aa);
k1 := one; l1 := - a / b;
loop
r := Rmd(aa,bb);
exit when Equal(r,zero);
h := aa / bb;
-- nk1 := k0 - h*k1; Clear(k0); k0 := k1; k1 := nk1;
-- nl1 := l0 - h*l1; Clear(l0); l0 := l1; l1 := nl1;
tmp := k1; k1 := k0 - h*k1; k0 := tmp;
tmp := l1; l1 := l0 - h*l1; l0 := tmp;
Clear(h);
Copy(bb,aa);
Copy(r,bb);
Clear(r);
end loop;
k := k1; l := l1; d := bb;
end if;
end pos_gcd;
begin
if Equal(a,zero)
then if b < zero
then d := -b; k := zero; l := -one;
else Copy(b,d); k := zero; l := one;
end if;
return;
elsif Equal(b,zero)
then if a < zero
then d := -a; k := -one; l := zero;
else Copy(a,d); k := one; l := zero;
end if;
return;
end if;
if a < b
then if b < zero
then minb := -b; mina := -a;
pos_gcd(minb,mina,ll,kk,d);
Clear(minb); Clear(mina);
Min(kk); k := kk;
Min(ll); l := ll;
elsif a < zero
then mina := -a;
pos_gcd(b,mina,l,kk,d);
Clear(mina);
Min(kk); k := kk;
else pos_gcd(b,a,l,k,d);
end if;
else -- a >= b
if a < zero
then mina := -a; minb := -b;
pos_gcd(mina,minb,kk,ll,d);
Clear(mina); Clear(minb);
Min(kk); k := kk;
Min(ll); l := ll;
elsif b < zero
then minb := -b;
pos_gcd(a,minb,k,ll,d);
Clear(minb);
Min(ll); l := ll;
else pos_gcd(a,b,k,l,d);
end if;
end if;
end gcd;
end Greatest_Common_Divisors;