Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Matrices/greatest_common_divisors.adb, Revision 1.1
1.1 ! maekawa 1: package body Greatest_Common_Divisors is
! 2:
! 3: use Euclidean_Domain;
! 4:
! 5: function pos_gcd ( a,b : number ) return number is
! 6:
! 7: res,r : number;
! 8:
! 9: begin
! 10: if Equal(b,zero)
! 11: then Copy(a,res);
! 12: else r := Rmd(a,b);
! 13: if Equal(r,zero)
! 14: then Copy(b,res);
! 15: else res := pos_gcd(b,r);
! 16: end if;
! 17: Clear(r);
! 18: end if;
! 19: return res;
! 20: end pos_gcd;
! 21:
! 22: function gcd ( a,b : number ) return number is
! 23:
! 24: res,mina,minb : number;
! 25:
! 26: begin
! 27: if a < zero
! 28: then if b < zero
! 29: then mina := -a; minb := -b;
! 30: res := pos_gcd(mina,minb);
! 31: Clear(mina); Clear(minb);
! 32: else mina := -a;
! 33: res := pos_gcd(mina,b);
! 34: Clear(mina);
! 35: end if;
! 36: elsif b < zero
! 37: then minb := -b;
! 38: res := pos_gcd(a,minb);
! 39: Clear(minb);
! 40: else res := pos_gcd(a,b);
! 41: end if;
! 42: return res;
! 43: end gcd;
! 44:
! 45: function lcm ( a,b : number ) return number is
! 46:
! 47: res : number;
! 48: g : number := gcd(a,b);
! 49:
! 50: begin
! 51: if Equal(g,zero)
! 52: then res := zero;
! 53: else res := a/g;
! 54: Mul(res,b);
! 55: end if;
! 56: Clear(g);
! 57: return res;
! 58: end lcm;
! 59:
! 60: procedure gcd ( a,b : in number; k,l,d : out number ) is
! 61:
! 62: mina,minb,kk,ll : number;
! 63:
! 64: procedure pos_gcd ( a,b : in number; k,l,d : out number ) is
! 65:
! 66: r,aa,bb,k1,l1,k0,l0,h,tmp : number;
! 67: -- nk1,nl1 : number;
! 68:
! 69: -- REQUIRED :
! 70: -- a >= b and a > 0 and b > 0
! 71:
! 72: begin
! 73: r := Rmd(a,b);
! 74: if Equal(r,zero)
! 75: then k := zero; l := one; Copy(b,d);
! 76: else k0 := zero; l0 := one;
! 77: Copy(r,bb);
! 78: Copy(b,aa);
! 79: k1 := one; l1 := - a / b;
! 80: loop
! 81: r := Rmd(aa,bb);
! 82: exit when Equal(r,zero);
! 83: h := aa / bb;
! 84: -- nk1 := k0 - h*k1; Clear(k0); k0 := k1; k1 := nk1;
! 85: -- nl1 := l0 - h*l1; Clear(l0); l0 := l1; l1 := nl1;
! 86: tmp := k1; k1 := k0 - h*k1; k0 := tmp;
! 87: tmp := l1; l1 := l0 - h*l1; l0 := tmp;
! 88: Clear(h);
! 89: Copy(bb,aa);
! 90: Copy(r,bb);
! 91: Clear(r);
! 92: end loop;
! 93: k := k1; l := l1; d := bb;
! 94: end if;
! 95: end pos_gcd;
! 96:
! 97: begin
! 98: if Equal(a,zero)
! 99: then if b < zero
! 100: then d := -b; k := zero; l := -one;
! 101: else Copy(b,d); k := zero; l := one;
! 102: end if;
! 103: return;
! 104: elsif Equal(b,zero)
! 105: then if a < zero
! 106: then d := -a; k := -one; l := zero;
! 107: else Copy(a,d); k := one; l := zero;
! 108: end if;
! 109: return;
! 110: end if;
! 111: if a < b
! 112: then if b < zero
! 113: then minb := -b; mina := -a;
! 114: pos_gcd(minb,mina,ll,kk,d);
! 115: Clear(minb); Clear(mina);
! 116: Min(kk); k := kk;
! 117: Min(ll); l := ll;
! 118: elsif a < zero
! 119: then mina := -a;
! 120: pos_gcd(b,mina,l,kk,d);
! 121: Clear(mina);
! 122: Min(kk); k := kk;
! 123: else pos_gcd(b,a,l,k,d);
! 124: end if;
! 125: else -- a >= b
! 126: if a < zero
! 127: then mina := -a; minb := -b;
! 128: pos_gcd(mina,minb,kk,ll,d);
! 129: Clear(mina); Clear(minb);
! 130: Min(kk); k := kk;
! 131: Min(ll); l := ll;
! 132: elsif b < zero
! 133: then minb := -b;
! 134: pos_gcd(a,minb,k,ll,d);
! 135: Clear(minb);
! 136: Min(ll); l := ll;
! 137: else pos_gcd(a,b,k,l,d);
! 138: end if;
! 139: end if;
! 140: end gcd;
! 141:
! 142: end Greatest_Common_Divisors;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>