![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to degree_structure.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Product|
Return to degree_structure.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Product |
1.1 maekawa 1: with unchecked_deallocation;
2: with text_io; use text_io;
3: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
4: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
5: with Standard_Complex_Numbers_io; use Standard_Complex_Numbers_io;
6: with Standard_Random_Numbers; use Standard_Random_Numbers;
7: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
8: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
9: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
10: with Degrees_in_Sets_of_Unknowns; use Degrees_in_Sets_of_Unknowns;
11: with Generate_Unions;
12:
13: package body Degree_Structure is
14:
15: -- DECLARATIONS :
16:
17: type zd(m : natural) is record
18: z : Partition(1..m);
19: d : Standard_Natural_Vectors.Vector(1..m);
20: end record;
21: type Link_To_zd is access zd;
22: procedure free is new unchecked_deallocation(zd,Link_To_zd);
23:
24: type dgst is array(positive range <>) of Link_To_zd;
25: type Link_To_dgst is access dgst;
26: procedure free is new unchecked_deallocation(dgst,Link_To_dgst);
27:
28: -- INTERNAL DATA :
29:
30: ds : Link_To_dgst;
31:
32: -- CREATORS :
33:
34: procedure Find_Partition ( p : in Poly;
35: z : in out Partition; m : in out natural;
36: dg : in out Standard_Natural_Vectors.Vector ) is
37:
38: -- DESCRIPTION :
39: -- This routine finds an good partition for the polynomial p
40:
41: -- ON ENTRY :
42: -- p a polynomial.
43:
44: -- ON RETURN :
45: -- z a partition of the set of unknowns of p;
46: -- m the number of sets in the partition z;
47: -- dg the degrees of the polynomial p in the sets of z,
48: -- dg(i) = Degree(p,z(i)).
49:
50: n : natural := Number_Of_Unknowns(p);
51: di : integer;
52: added : boolean;
53:
54: begin
55: for i in 1..n loop
56: di := Degree(p,i);
57: if di > 0
58: then added := false;
59: for j in 1..m loop
60: if di = dg(j)
61: then Add(z(j),i);
62: if Degree(p,z(j)) = dg(j)
63: then added := true;
64: else Remove(z(j),i);
65: end if;
66: end if;
67: exit when added;
68: end loop;
69: if not added
70: then m := m + 1;
71: Add(z(m),i);
72: dg(m) := di;
73: end if;
74: end if;
75: end loop;
76: end Find_Partition;
77:
78: procedure Create ( p : in Poly_Sys ) is
79:
80: n : natural := p'length;
81: z : Partition(1..n);
82: d : Standard_Natural_Vectors.Vector(1..n) := (1..n => 0);
83: m : natural;
84:
85: begin
86: if ds /= null
87: then Clear;
88: end if;
89: ds := new dgst(1..n);
90: for i in p'range loop
91: m := 0;
92: Create(z,n);
93: Find_Partition(p(i),z,m,d);
94: ds(i) := new zd(m);
95: for j in 1..m loop
96: ds(i).z(j) := Create(z(j));
97: ds(i).d(j) := d(j);
98: end loop;
99: Clear(z);
100: end loop;
101: end Create;
102:
103: procedure Put ( p : in Poly_Sys;
104: i,m : in natural; z : in Partition ) is
105:
106: n : natural := p'length;
107:
108: begin
109: if ds = null
110: then ds := new dgst(1..n);
111: end if;
112: ds(i) := new zd(m);
113: for j in 1..m loop
114: ds(i).z(j) := Create(z(j));
115: ds(i).d(j) := Degree(p(i),z(j));
116: end loop;
117: end Put;
118:
119: -- SELECTORS :
120:
121: function Empty return boolean is
122: begin
123: return (ds = null);
124: end Empty;
125:
126: function Get ( i : natural ) return natural is
127: begin
128: return ds(i).m;
129: end Get;
130:
131: procedure Get ( i : in natural; z : in out Partition;
132: d : out Standard_Natural_Vectors.Vector ) is
133: begin
134: for j in 1..ds(i).m loop
135: z(j) := Create(ds(i).z(j));
136: d(j) := ds(i).d(j);
137: end loop;
138: end Get;
139:
140: -- COMPUTING THE GENERALIZED BEZOUT NUMBER :
141:
142: function Matrix_Criterion ( z : Partition ) return boolean is
143:
144: -- DESCRIPTION :
145: -- This is the Matrix criterion for testing if
146: -- a product is admissible or not.
147:
148: n : natural := z'last - z'first +1;
149: mat : Matrix(1..n,1..n);
150: ipvt : Standard_Natural_Vectors.Vector(1..n);
151: eps : constant double_float := 10.0**(-10);
152: r : Complex_Number;
153: rcond : double_float;
154:
155: begin
156: for i in 1..n loop
157: r := Create(double_float(i+1));
158: for j in 1..n loop
159: if Is_In(z(i),j)
160: then mat(i,j) := r;
161: r := r*r;
162: else mat(i,j) := Create(0.0);
163: end if;
164: end loop;
165: end loop;
166: lufco(mat,n,ipvt,rcond);
167: return (abs(rcond) > eps);
168: exception
169: when others => return false;
170: end Matrix_Criterion;
171:
172: function Admissible ( z : Partition; n : natural ) return boolean is
173:
174: temp : Partition(1..n);
175: admis : boolean := true;
176:
177: begin
178: temp(1) := Create(z(1));
179: for i in 2..(n-1) loop
180: temp(i) := Create(z(i));
181: admis := Admissible(temp,i,z(i+1));
182: exit when not admis;
183: end loop;
184: Clear(temp);
185: return admis;
186: end Admissible;
187:
188: function Admissible ( z : Partition; n : natural; s : Set )
189: return boolean is
190: begin
191: for k in 1..n loop
192: if not Admissible(z,k,n,s)
193: then return false;
194: end if;
195: end loop;
196: return true;
197: end Admissible;
198:
199: function Admissible ( z : Partition; k,n : natural; s : Set )
200: return boolean is
201:
202: type arr is array (integer range <>) of boolean;
203: admis : boolean := true;
204:
205: procedure check (a : in arr; continue : out boolean) is
206:
207: u : Set := Create(s);
208:
209: begin
210: for i in a'range loop
211: if a(i)
212: then Union(u,z(i));
213: end if;
214: end loop;
215: admis := ( Extent_Of(u) >= k+1 );
216: continue := admis;
217: Clear(u);
218: end check;
219:
220: procedure gen is new Generate_Unions(arr,check);
221:
222: begin
223: gen(k,1,n);
224: return admis;
225: end Admissible;
226:
227: procedure Compute ( i,n,sum : in natural; res : in out natural;
228: z : in out Partition ) is
229: begin
230: if i > n
231: then res := res + sum;
232: else -- Pick out a set and check if it is allowed :
233: for j in 1..ds(i).m loop
234: if ds(i).d(j) /= 0 and then Admissible(z,i-1,ds(i).z(j))
235: then z(i) := Create(ds(i).z(j));
236: Compute(i+1,n,sum*ds(i).d(j),res,z);
237: Clear(z(i));
238: end if;
239: end loop;
240: end if;
241: end Compute;
242:
243: function Generalized_Bezout_Number return natural is
244:
245: res : natural := 0;
246: n : natural := ds'length;
247: z : Partition(1..n);
248:
249: begin
250: Compute(1,n,1,res,z);
251: return res;
252: end Generalized_Bezout_Number;
253:
254: function Generalized_Bezout_Number ( p : in Poly_Sys ) return natural is
255: begin
256: Create(p);
257: return Generalized_Bezout_Number;
258: end Generalized_Bezout_Number;
259:
260: -- DESTRUCTOR :
261:
262: procedure Clear is
263: begin
264: if ds /= null
265: then for i in ds'range loop
266: free(ds(i));
267: end loop;
268: free(ds);
269: end if;
270: end Clear;
271:
272: end Degree_Structure;