Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Matrices/greatest_common_divisors.adb, Revision 1.1.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>