Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Product/m_homogeneous_bezout_numbers.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Integer_Vectors; use Standard_Integer_Vectors;
2: with Standard_Integer_Matrices; use Standard_Integer_Matrices;
3: with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
4: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
5: with Sets_of_Unknowns; use Sets_of_Unknowns;
6: with Degrees_in_Sets_of_Unknowns; use Degrees_in_Sets_of_Unknowns;
7:
8: package body m_Homogeneous_Bezout_Numbers is
9:
10: -- UTILITIES :
11:
12: function Create ( p : Poly_Sys ) return Set is
13:
14: -- DESCRIPTION :
15: -- Returns the set of the unknowns of the polynomial system p.
16:
17: s : Set := Create(p'length);
18:
19: begin
20: for i in p'range loop
21: Add(s,i);
22: end loop;
23: return s;
24: end Create;
25:
26: function Cardinalities ( z : Partition ) return Vector is
27:
28: -- DESCRIPTION :
29: -- Returns a vector which contains the cardinality of each set.
30:
31: res : Vector(z'range);
32:
33: begin
34: for i in z'range loop
35: res(i) := Extent_Of(z(i));
36: end loop;
37: return res;
38: end Cardinalities;
39:
40: -- TARGET ROUTINES :
41:
42: function Total_Degree ( p : Poly_Sys ) return natural is
43:
44: d : natural := 1;
45:
46: begin
47: for i in p'range loop
48: d := d*Degree(p(i));
49: end loop;
50: return d;
51: end Total_Degree;
52:
53: function Bezout_Number ( p : Poly_Sys; z : Partition ) return natural is
54:
55: k : constant vector := Cardinalities(z);
56: d : constant matrix := Degree_Table(p,z);
57:
58: begin
59: return Per(d,k); -- returns the permanent of the degree table
60: end Bezout_Number;
61:
62: function Bezout_Number ( p : Poly_Sys; z : Partition; max : natural )
63: return natural is
64:
65: -- DESCRIPTION :
66: -- Stops when the Bezout number becomes bigger than max.
67:
68: k : constant vector := Cardinalities(z);
69: d : constant matrix := Degree_Table(p,z);
70:
71: begin
72: return Per(d,k,max); -- returns the permanent of the degree table
73: end Bezout_Number;
74:
75: function Bezout_Number ( p : Poly_Sys ) return natural is
76:
77: s : Set := Create(p);
78: res : natural := Total_Degree(p);
79:
80: procedure Evaluate ( z : in Partition; cont : out boolean ) is
81: b : constant natural := Bezout_Number(p,z,res);
82: begin
83: if b < res
84: then res := b;
85: end if;
86: cont := true;
87: end Evaluate;
88: procedure Evaluate_Partitions is new Generate_Partitions(Evaluate);
89:
90: begin
91: Evaluate_Partitions(s);
92: Clear(s);
93: return res;
94: end Bezout_Number;
95:
96: function Bezout_Number ( max : natural; p : Poly_Sys ) return natural is
97:
98: s : Set := Create(p);
99: res : natural := Total_Degree(p);
100: cnt : natural := 0;
101:
102: procedure Evaluate ( z : in Partition; cont : out boolean ) is
103: b : constant natural := Bezout_Number(p,z,res);
104: begin
105: if b < res
106: then res := b;
107: end if;
108: cnt := cnt + 1;
109: if cnt < max
110: then cont := true;
111: else cont := false;
112: end if;
113: end Evaluate;
114: procedure Evaluate_Partitions is new Generate_Partitions(Evaluate);
115:
116: begin
117: Evaluate_Partitions(s);
118: Clear(s);
119: return res;
120: end Bezout_Number;
121:
122: function Bezout_Number ( p : Poly_Sys; min : natural ) return natural is
123:
124: s : Set := Create(p);
125: res : natural := Total_Degree(p);
126:
127: procedure Evaluate ( z : in Partition; cont : out boolean ) is
128: b : constant natural := Bezout_Number(p,z,res);
129: begin
130: if b < res
131: then res := b;
132: end if;
133: if res > min
134: then cont := true;
135: else cont := false;
136: end if;
137: end Evaluate;
138: procedure Evaluate_Partitions is new Generate_Partitions(Evaluate);
139:
140: begin
141: Evaluate_Partitions(s);
142: Clear(s);
143: return res;
144: end Bezout_Number;
145:
146: function Bezout_Number ( max : natural; p : Poly_Sys; min : natural )
147: return natural is
148:
149: s : Set := Create(p);
150: res : natural := Total_Degree(p);
151: cnt : natural := 0;
152:
153: procedure Evaluate ( z : in Partition; cont : out boolean ) is
154: b : constant natural := Bezout_Number(p,z,res);
155: begin
156: if b < res
157: then res := b;
158: end if;
159: cnt := cnt + 1;
160: if cnt < max and then res > min
161: then cont := true;
162: else cont := false;
163: end if;
164: end Evaluate;
165: procedure Evaluate_Partitions is new Generate_Partitions(Evaluate);
166:
167: begin
168: Evaluate_Partitions(s);
169: Clear(s);
170: return res;
171: end Bezout_Number;
172:
173: procedure Bezout_Number
174: ( p : in Poly_Sys; b,m : out natural; z : in out Partition ) is
175:
176: s : Set := Create(p);
177: tdg : constant natural := Total_Degree(p);
178: res : natural := tdg;
179:
180: procedure Evaluate ( nz : in Partition; cont : out boolean ) is
181: nb : constant natural := Bezout_Number(p,nz,res);
182: begin
183: if nb < res
184: then res := nb;
185: m := nz'length; Clear(z);
186: z(nz'range) := Create(nz);
187: end if;
188: cont := true;
189: end Evaluate;
190: procedure Evaluate_Partitions is new Generate_Partitions(Evaluate);
191:
192: begin
193: Evaluate_Partitions(s);
194: if res = tdg
195: then m := 1; z(1) := s;
196: else Clear(s);
197: end if;
198: b := res;
199: end Bezout_Number;
200:
201: procedure Bezout_Number
202: ( max : in natural; p : in Poly_Sys; b,m : out natural;
203: z : in out Partition ) is
204:
205: s : Set := Create(p);
206: tdg : constant natural := Total_Degree(p);
207: res : natural := tdg;
208: cnt : natural := 0;
209:
210: procedure Evaluate ( nz : in Partition; cont : out boolean ) is
211: nb : constant natural := Bezout_Number(p,nz,res);
212: begin
213: if nb < res
214: then res := nb;
215: m := nz'length; Clear(z);
216: z(nz'range) := Create(nz);
217: end if;
218: cnt := cnt + 1;
219: if cnt < max
220: then cont := true;
221: else cont := false;
222: end if;
223: end Evaluate;
224: procedure Evaluate_Partitions is new Generate_Partitions(Evaluate);
225:
226: begin
227: Evaluate_Partitions(s);
228: if res = tdg
229: then m := 1; z(1) := s;
230: else Clear(s);
231: end if;
232: b := res;
233: end Bezout_Number;
234:
235: procedure Bezout_Number
236: ( p : in Poly_Sys; min : in natural; b,m : out natural;
237: z : in out Partition ) is
238:
239: s : Set := Create(p);
240: tdg : constant natural := Total_Degree(p);
241: res : natural := tdg;
242:
243: procedure Evaluate ( nz : in Partition; cont : out boolean ) is
244: nb : constant natural := Bezout_Number(p,nz,res);
245: begin
246: if nb < res
247: then res := nb;
248: m := nz'length; Clear(z);
249: z(nz'range) := Create(nz);
250: end if;
251: if res > min
252: then cont := true;
253: else cont := false;
254: end if;
255: end Evaluate;
256: procedure Evaluate_Partitions is new Generate_Partitions(Evaluate);
257:
258: begin
259: Evaluate_Partitions(s);
260: if res = tdg
261: then m := 1; z(1) := s;
262: else Clear(s);
263: end if;
264: b := res;
265: end Bezout_Number;
266:
267: procedure Bezout_Number
268: ( max : in natural; p : in Poly_Sys; min : in natural;
269: b,m : out natural; z : in out Partition ) is
270:
271: s : Set := Create(p);
272: tdg : constant natural := Total_Degree(p);
273: res : natural := tdg;
274: cnt : natural := 0;
275:
276: procedure Evaluate ( nz : in Partition; cont : out boolean ) is
277: nb : constant natural := Bezout_Number(p,nz,res);
278: begin
279: if nb < res
280: then res := nb;
281: m := nz'length; Clear(z);
282: z(nz'range) := Create(nz);
283: end if;
284: cnt := cnt + 1;
285: if cnt < max and then res > min
286: then cont := true;
287: else cont := false;
288: end if;
289: end Evaluate;
290: procedure Evaluate_Partitions is new Generate_Partitions(Evaluate);
291:
292: begin
293: Evaluate_Partitions(s);
294: if res = tdg
295: then m := 1; z(1) := s;
296: else Clear(s);
297: end if;
298: b := res;
299: end Bezout_Number;
300:
301: function Evaluate ( z : partition; m : natural; p : Poly_Sys )
302: return natural is
303:
304: n : natural := p'length;
305: d : constant matrix := Degree_Table(p,z);
306:
307: function Bezout_number
308: ( n,m : natural; z: partition; d : matrix ) return natural is
309:
310: -- DESCRIPTION : the Bezout number is computed
311:
312: type boolean_array is array ( positive range <> ) of boolean;
313: Ii,Iacc : boolean_array(1..n) := (1..n => false);
314: b,b_mult : natural;
315:
316: procedure column ( j,start,number : in natural;
317: Ii,Iacc : in out boolean_array );
318:
319: -- DESCRIPTION : the computation of a term coming from a column
320: -- of the degree table;
321:
322: procedure row ( j : in natural; Ii,Iacc : in out boolean_array );
323:
324: -- DESCRIPTION : the computation of a row in the degree table
325:
326: procedure column ( j,start,number : in natural;
327: Ii,Iacc : in out boolean_array ) is
328: begin
329: if number > (n - start + 1)
330: then return;
331: elsif number = 0
332: then row(j,Ii,Iacc);
333: else for i in start..n loop
334: if not Ii(i)
335: then Ii(i) := true; Iacc(i) := true;
336: column(j,i+1,number-1,Ii,Iacc);
337: Ii(i) := false; Iacc(i) := false;
338: end if;
339: end loop;
340: end if;
341: end column;
342:
343: procedure row ( j : in natural; Ii,Iacc : in out boolean_array ) is
344: temp : natural := 1;
345: Iacc1 : boolean_array(1..n) := (1..n => false);
346: begin
347: for k in 1..n loop
348: if Iacc(k)
349: then temp := temp * d(k,j);
350: end if;
351: end loop;
352: if (j /= m) and (temp /= 0)
353: then b_mult := b_mult * temp;
354: column(j+1,1,Extent_Of(z(j+1)),Ii,Iacc1);
355: b_mult := b_mult / temp;
356: elsif j = m
357: then temp := temp * b_mult;
358: b := b + temp;
359: end if;
360: end row;
361:
362: begin
363: b := 0; b_mult := 1;
364: column(1,1,Extent_Of(z(1)),Ii,Iacc);
365: return b;
366: end Bezout_number;
367:
368: begin
369: return Bezout_number(n,m,z,d);
370: end Evaluate;
371:
372: procedure PB ( p : in Poly_Sys; b,m : out natural; z : in out Partition ) is
373:
374: n : natural := p'length;
375: wb,b_min : natural;
376: wz,z_min : partition(1..n);
377: wm,m_min : natural := 0;
378:
379: procedure pcopy ( p1,p2 : in out partition; p1n,p2n : in out natural ) is
380: -- DESCRIPTION : the partition p1 is copied to p2
381: begin
382: Clear(p2); p2 := Create(p1);
383: p2n := p1n;
384: end pcopy;
385:
386: begin
387: b_min := Total_Degree(p);
388: for i in 1..n loop
389: for k in 1..wm loop
390: Add(wz(k),i);
391: wb := Evaluate(wz,wm,p);
392: if (k = 1) or else (wb < b_min)
393: then pcopy(wz,z_min,wm,m_min);
394: b_min := wb;
395: end if;
396: Remove(wz(k),i);
397: end loop;
398: wm := wm + 1;
399: wz(wm) := Create(n);
400: Add(wz(wm),i);
401: wb := Evaluate(wz,wm,p);
402: if wb < b_min
403: then pcopy(wz,z_min,wm,m_min);
404: b_min := wb;
405: end if;
406: pcopy(z_min,wz,m_min,wm);
407: end loop;
408: b := b_min;
409: m := m_min;
410: pcopy(z_min,z,m_min,wm);
411: Clear(wz);
412: end PB;
413:
414: end m_Homogeneous_Bezout_Numbers;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to m_homogeneous_bezout_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 / Root_Counts / Product