![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to set_structures_and_volumes.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 / Implift|
Return to set_structures_and_volumes.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift |
1.1 maekawa 1: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
2: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
3: with Standard_Natural_Vectors;
4: with Standard_Complex_Vectors; use Standard_Complex_Vectors;
5: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
6: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
7: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
8: with Standard_Complex_Poly_Randomizers; use Standard_Complex_Poly_Randomizers;
9: with Standard_Complex_Substitutors; use Standard_Complex_Substitutors;
10: with Standard_Complex_Laur_Systems; use Standard_Complex_Laur_Systems;
11: with Standard_Poly_Laur_Convertors; use Standard_Poly_Laur_Convertors;
12: with Arrays_of_Integer_Vector_Lists; use Arrays_of_Integer_Vector_Lists;
13: with Set_Structure;
14: with Random_Product_System;
15: with Random_Product_Start_Systems;
16: with Power_Lists; use Power_Lists;
17: with Volumes; use Volumes;
18: with Mixed_Homotopy_Continuation; use Mixed_Homotopy_Continuation;
19:
20: package body Set_Structures_and_Volumes is
21:
22: -- AUXILIARIES :
23:
24: procedure Build_RPS ( k,n : in natural; p : in Poly_Sys ) is
25:
26: -- DESCRIPTION :
27: -- If the set structure is empty, then a set structure will
28: -- be constructed for the first k polynomials of p;
29: -- the data in Random_Product_System will be build.
30:
31: begin
32: if Set_Structure.Empty
33: then declare
34: tmp : Poly_Sys(p'range);
35: t : Term;
36: cnst : Poly;
37: begin
38: tmp(1..k) := p(1..k);
39: t.dg := new Standard_Natural_Vectors.Vector'(1..n => 0);
40: t.cf := Create(1.0);
41: cnst := Create(t);
42: tmp(k+1..n) := (k+1..n => cnst);
43: Random_Product_Start_Systems.Build_Set_Structure(tmp);
44: Standard_Natural_Vectors.Clear
45: (Standard_Natural_Vectors.Link_to_Vector(t.dg));
46: Clear(cnst);
47: end;
48: end if;
49: Random_Product_System.Init(n);
50: Random_Product_Start_Systems.Build_Random_Product_System(n);
51: end Build_RPS;
52:
53: procedure Build_Elimination_Matrix
54: ( k,n : in natural; ind : in Standard_Natural_Vectors.Vector;
55: m : in out Matrix;
56: pivots : in out Standard_Natural_Vectors.Vector;
57: degenerate : out boolean ) is
58:
59: -- DESCRIPTION :
60: -- builds the elimination matrix m.
61:
62: -- ON ENTRY :
63: -- k the number of unknowns to be eliminated;
64: -- n the number of polynomials and unknowns;
65: -- ind entries in Random_Product_System;
66: -- ind(l) indicates lth hyperplane;
67:
68: -- ON RETURN :
69: -- m contains all k hyperplanes to use for elimination;
70: -- pivots if not degenerate, then m(i,pivots(i)) /= 0;
71: -- degenerate is true if the first k hyperplanes were inconsistent.
72:
73: degen : boolean;
74:
75: begin
76: for i in ind'range loop -- build the matrix
77: declare
78: h : Vector(0..n) := Random_Product_System.Get_Hyperplane(i,ind(i));
79: begin
80: for j in 1..n loop
81: m(i,j) := h(j);
82: end loop;
83: m(i,n+1) := h(0);
84: end;
85: end loop;
86: diagonalize(m,k,n+1);
87: degen := false; -- check degeneracy
88: declare
89: eps : constant double_float := 10.0**(-10);
90: begin
91: for i in pivots'range loop
92: for j in 1..n loop
93: if AbsVal(m(i,j)) > eps
94: then pivots(i) := j;
95: end if;
96: exit when (pivots(i) /= 0);
97: end loop;
98: degen := (pivots(i) = 0);
99: exit when degen;
100: end loop;
101: end;
102: degenerate := degen;
103: end Build_Elimination_Matrix;
104:
105: procedure Eliminate ( k,n : in natural; p : in Poly_Sys; m : in Matrix;
106: pivots : in Standard_Natural_Vectors.Vector;
107: q : in out Poly_Sys ) is
108: -- DESCRIPTION :
109: -- eliminates k unknowns of the polynomial system p.
110:
111: -- ON ENTRY :
112: -- k the number of unknowns to be eliminated;
113: -- n the number of polynomials and unknowns;
114: -- p a polynomial system with random coefficients;
115: -- m the diagonalized matrix, to be used for elimination,
116: -- for i in m'range(1) : m(i,pivots(i)) /= 0;
117: -- pivots a vector indicating nonzero entries in m.
118:
119: -- ON RETURN :
120: -- q the reduced system.
121:
122: begin
123: Clear(q); Copy(p,q);
124: for i in pivots'range loop
125: declare
126: h : Vector(0..n-i+1);
127: piv : natural := pivots(i)-i+1;
128: begin
129: h(0) := m(i,n+1);
130: h(1..n-i+1) := (1..n-i+1 => Create(0.0));
131: h(piv) := m(i,pivots(i));
132: for j in piv+1..n-i+1 loop
133: h(j) := m(i,j+i-1);
134: end loop;
135: Substitute(piv,h,q);
136: end;
137: end loop;
138: end Eliminate;
139:
140: procedure Update ( sols : in out Solution_List; k,n : in natural;
141: m : in Matrix;
142: pivots : in Standard_Natural_Vectors.Vector ) is
143:
144: -- DESCRIPTION :
145: -- Based on the elimination matrix m, solution components for x1..xk
146: -- can be computed.
147:
148: -- ON ENTRY :
149: -- sols a list of solutions, with values for x(n-k+1)..x(n);
150: -- k the number of unknowns that have been eliminated;
151: -- n total number of unknowns;
152: -- m the elimination matrix;
153: -- pivots for i in pivots'range : m(i,pivots(i)) /= 0.
154:
155: -- ON RETURN :
156: -- sols the solution list with n-dimensional vectors.
157:
158: tmp,res,res_last : Solution_List;
159:
160: begin
161: tmp := sols;
162: while not Is_Null(tmp) loop
163: declare
164: sol : Solution(n-k) := Head_Of(tmp).all;
165: soln : Solution(n);
166: ind : natural;
167: begin
168: soln.t := Create(0.0);
169: soln.m := sol.m;
170: soln.v := (1..n => Create(0.0));
171: for i in 1..k loop
172: soln.v(pivots(i)) := -m(i,n+1);
173: end loop;
174: ind := 0;
175: for i in 1..n loop
176: if soln.v(i) = Create(0.0)
177: then ind := ind + 1;
178: soln.v(i) := sol.v(ind);
179: end if;
180: end loop;
181: for i in 1..k loop
182: ind := pivots(i);
183: for j in ind+1..n loop
184: soln.v(ind) := soln.v(ind) - m(i,j)*soln.v(j);
185: end loop;
186: soln.v(ind) := soln.v(ind) / m(i,ind);
187: end loop;
188: Append(res,res_last,soln);
189: end;
190: tmp := Tail_Of(tmp);
191: end loop;
192: Clear(sols); Copy(res,sols);
193: end Update;
194:
195: procedure Build_Indices ( l,k,n : in natural; p : in Poly_Sys;
196: acc : in out Standard_Natural_Vectors.Vector;
197: bb : in out natural ) is
198:
199: -- DESCRIPTION :
200: -- Builds the indices in the vector acc.
201:
202: -- ON ENTRY :
203: -- l index in the array acc;
204: -- k number of unknowns to eliminate;
205: -- n number of polynomials and unknowns;
206: -- p a polynomial system;
207: -- acc(1..l-1) contains entries indicating hyperplanes
208: -- in Random_Product_System.
209:
210: -- ON RETURN :
211: -- acc(1..l) contains entries indicating hyperplanes
212: -- in Random_Product_System;
213: -- bb the mixed Bezout BKK bound.
214:
215: begin
216: if l > k
217: then declare
218: degenerate : boolean;
219: m : Matrix(1..k,1..n+1);
220: pivots : Standard_Natural_Vectors.Vector(1..k);
221: begin
222: pivots := (1..k => 0);
223: Build_Elimination_Matrix(k,n,acc,m,pivots,degenerate);
224: if not degenerate
225: then declare
226: q : Poly_Sys(p'range);
227: begin
228: Eliminate(k,n,p,m,pivots,q);
229: bb := bb + Bezout_BKK(0,n-k,q);
230: Clear(q);
231: end;
232: end if;
233: end;
234: else for j in 1..Random_Product_System.Number_of_Hyperplanes(l) loop
235: acc(l) := j;
236: Build_Indices(l+1,k,n,p,acc,bb);
237: end loop;
238: end if;
239: end Build_Indices;
240:
241: procedure Build_Indices ( l,k,n : in natural; p : in Poly_Sys;
242: tv : in Tree_of_Vectors;
243: acc : in out Standard_Natural_Vectors.Vector;
244: bb : in out natural ) is
245:
246: -- DESCRIPTION :
247: -- Builds the indices in the vector acc.
248:
249: -- ON ENTRY :
250: -- l index in the array acc;
251: -- k number of unknowns to eliminate;
252: -- n number of polynomials and unknowns;
253: -- p a polynomial system;
254: -- tv the tree of degrees used to compute the mixed volume;
255: -- acc(1..l-1) contains entries indicating hyperplanes
256: -- in Random_Product_System.
257:
258: -- ON RETURN :
259: -- acc(1..l) contains entries indicating hyperplanes
260: -- in Random_Product_System;
261: -- bb the mixed Bezout BKK bound.
262:
263: begin
264: if l > k
265: then declare
266: degenerate : boolean;
267: m : Matrix(1..k,1..n+1);
268: pivots : Standard_Natural_Vectors.Vector(1..k);
269: begin
270: pivots := (1..k => 0);
271: Build_Elimination_Matrix(k,n,acc,m,pivots,degenerate);
272: if not degenerate
273: then declare
274: q : Poly_Sys(p'range);
275: begin
276: Eliminate(k,n,p,m,pivots,q);
277: bb := bb + Bezout_BKK(0,n-k,q,tv);
278: Clear(q);
279: end;
280: end if;
281: end;
282: else for j in 1..Random_Product_System.Number_of_Hyperplanes(l) loop
283: acc(l) := j;
284: Build_Indices(l+1,k,n,p,tv,acc,bb);
285: end loop;
286: end if;
287: end Build_Indices;
288:
289: procedure Build_Indices2 ( l,k,n : in natural; p : in Poly_Sys;
290: tv : in out Tree_of_Vectors;
291: acc : in out Standard_Natural_Vectors.Vector;
292: bb : in out natural ) is
293:
294: -- DESCRIPTION :
295: -- Builds the indices in the vector acc.
296:
297: -- ON ENTRY :
298: -- l index in the array acc;
299: -- k number of unknowns to eliminate;
300: -- n number of polynomials and unknowns;
301: -- p a polynomial system;
302: -- acc(1..l-1) contains entries indicating hyperplanes
303: -- in Random_Product_System.
304:
305: -- ON RETURN :
306: -- tv the tree of degrees used to compute the mixed volume;
307: -- acc(1..l) contains entries indicating hyperplanes
308: -- in Random_Product_System;
309: -- bb the mixed Bezout BKK bound.
310:
311: begin
312: if l > k
313: then declare
314: degenerate : boolean;
315: m : Matrix(1..k,1..n+1);
316: pivots : Standard_Natural_Vectors.Vector(1..k);
317: begin
318: pivots := (1..k => 0);
319: Build_Elimination_Matrix(k,n,acc,m,pivots,degenerate);
320: if not degenerate
321: then declare
322: q : Poly_Sys(p'range);
323: qtv,tmp : Tree_of_Vectors;
324: mv : natural;
325: begin
326: Eliminate(k,n,p,m,pivots,q);
327: Bezout_BKK(0,n-k,q,qtv,mv);
328: bb := bb + mv;
329: tmp := qtv;
330: while not Is_Null(tmp) loop
331: declare
332: nd : node := Head_Of(tmp);
333: begin
334: if Is_In(tv,nd.d)
335: then Clear(nd);
336: else Construct(nd,tv);
337: end if;
338: end;
339: tmp := Tail_Of(tmp);
340: end loop;
341: Clear(q);
342: end;
343: end if;
344: end;
345: else for j in 1..Random_Product_System.Number_of_Hyperplanes(l) loop
346: acc(l) := j;
347: Build_Indices2(l+1,k,n,p,tv,acc,bb);
348: end loop;
349: end if;
350: end Build_Indices2;
351:
352: procedure Build_Indices
353: ( file : in file_type; p : in Poly_Sys; l,k,n : in natural;
354: acc : in out Standard_Natural_Vectors.Vector;
355: bb : in out natural;
356: sols,sols_last : in out Solution_List ) is
357:
358: -- DESCRIPTION :
359: -- Builds the indices in the vector acc.
360:
361: -- ON ENTRY :
362: -- file file to write intermediate results on;
363: -- p a polynomial system;
364: -- l index in the array acc;
365: -- k number of unknowns to eliminate;
366: -- n number of polynomials and unknowns;
367: -- acc(1..l-1) contains entries indicating hyperplanes
368: -- in Random_Product_System.
369:
370: -- ON RETURN :
371: -- acc(1..l) contains entries indicating hyperplanes
372: -- in Random_Product_System;
373: -- bb the mixed Bezout BKK bound;
374: -- sols the solutions of pi;
375: -- sols_last points to the last element of the list sols.
376:
377: begin
378: if l > k
379: then declare
380: degenerate : boolean;
381: m : Matrix(1..k,1..n+1);
382: pivots : Standard_Natural_Vectors.Vector(1..k);
383: begin
384: pivots := (1..k => 0);
385: Build_Elimination_Matrix(k,n,acc,m,pivots,degenerate);
386: if not degenerate
387: then declare
388: q : Poly_Sys(p'range);
389: begin
390: Eliminate(k,n,p,m,pivots,q);
391: declare
392: las : Laur_Sys(q'range);
393: qsols : Solution_List;
394: bkk : natural;
395: begin
396: las := Polynomial_to_Laurent_System(q);
397: Solve(file,las,bkk,qsols);
398: bb := bb + bkk;
399: Update(qsols,k,n,m,pivots);
400: Concat(sols,sols_last,qsols);
401: Clear(las); Shallow_Clear(qsols);
402: end;
403: Clear(q);
404: end;
405: end if;
406: end;
407: else for j in 1..Random_Product_System.Number_of_Hyperplanes(l) loop
408: acc(l) := j;
409: Build_Indices(file,p,l+1,k,n,acc,bb,sols,sols_last);
410: end loop;
411: end if;
412: end Build_Indices;
413:
414: procedure Build_Indices
415: ( file : in file_type; p : in Poly_Sys; l,k,n : in natural;
416: acc : in out Standard_Natural_Vectors.Vector;
417: tv : in Tree_of_Vectors; bb : in out natural;
418: sols,sols_last : in out Solution_List ) is
419:
420: -- DESCRIPTION :
421: -- Builds the indices in the vector acc.
422:
423: -- ON ENTRY :
424: -- file file to write intermediate results on;
425: -- p a polynomial system;
426: -- l index in the array acc;
427: -- k number of unknowns to eliminate;
428: -- n number of polynomials and unknowns;
429: -- acc(1..l-1) contains entries indicating hyperplanes
430: -- in Random_Product_System;
431: -- tv the tree containing useful directions;
432:
433: -- ON RETURN :
434: -- acc(1..l) contains entries indicating hyperplanes
435: -- in Random_Product_System;
436: -- bb the mixed Bezout BKK bound;
437: -- sols the solutions of pi;
438: -- sols_last points to the last element of the list sols.
439:
440: begin
441: if l > k
442: then declare
443: degenerate : boolean;
444: m : Matrix(1..k,1..n+1);
445: pivots : Standard_Natural_Vectors.Vector(1..k);
446: begin
447: pivots := (1..k => 0);
448: Build_Elimination_Matrix(k,n,acc,m,pivots,degenerate);
449: if not degenerate
450: then declare
451: q : Poly_Sys(p'range);
452: begin
453: Eliminate(k,n,p,m,pivots,q);
454: declare
455: las : Laur_Sys(q'range);
456: qsols : Solution_List;
457: bkk : natural;
458: begin
459: las := Polynomial_to_Laurent_System(q);
460: Solve(file,las,tv,bkk,qsols);
461: bb := bb + bkk;
462: Update(qsols,k,n,m,pivots);
463: Concat(sols,sols_last,qsols);
464: Clear(las); Shallow_Clear(qsols);
465: end;
466: Clear(q);
467: end;
468: end if;
469: end;
470: else for j in 1..Random_Product_System.Number_of_Hyperplanes(l) loop
471: acc(l) := j;
472: Build_Indices(file,p,l+1,k,n,acc,tv,bb,sols,sols_last);
473: end loop;
474: end if;
475: end Build_Indices;
476:
477: -- TARGET ROUTINES :
478:
479: function Bezout_BKK ( k,n : natural; p : Poly_Sys ) return natural is
480:
481: res : natural;
482:
483: begin
484: if k = 0
485: then declare
486: adl : Array_of_Lists(p'range) := Create(p);
487: begin
488: res := Mixed_Volume(n,adl);
489: Deep_Clear(adl);
490: end;
491: elsif k = n
492: then if Set_Structure.Empty
493: then Random_Product_Start_Systems.Build_Set_Structure(p);
494: res := Set_Structure.B;
495: Set_Structure.Clear;
496: else res := Set_Structure.B;
497: end if;
498: else Build_RPS(k,n,p);
499: declare
500: acc : Standard_Natural_Vectors.Vector(1..k);
501: q : Poly_Sys(1..n-k);
502: begin
503: acc := (1..k => 0); res := 0;
504: for i in q'range loop
505: q(i) := Complex_Randomize1(p(i+k));
506: end loop;
507: Build_Indices(1,k,n,q,acc,res);
508: Clear(q);
509: end;
510: Random_Product_System.Clear;
511: end if;
512: return res;
513: end Bezout_BKK;
514:
515: function Bezout_BKK ( k,n : natural; p : Poly_Sys; tv : Tree_of_Vectors )
516: return natural is
517: res : natural;
518:
519: begin
520: if k = 0
521: then declare
522: adl : Array_of_Lists(p'range) := Create(p);
523: begin
524: res := Mixed_Volume(n,adl,tv);
525: Deep_Clear(adl);
526: end;
527: elsif k = n
528: then if Set_Structure.Empty
529: then Random_Product_Start_Systems.Build_Set_Structure(p);
530: res := Set_Structure.B;
531: Set_Structure.Clear;
532: else res := Set_Structure.B;
533: end if;
534: else Build_RPS(k,n,p);
535: declare
536: acc : Standard_Natural_Vectors.Vector(1..k);
537: q : Poly_Sys(1..n-k);
538: begin
539: acc := (1..k => 0); res := 0;
540: for i in q'range loop
541: q(i) := Complex_Randomize1(p(i+k));
542: end loop;
543: Build_Indices(1,k,n,q,tv,acc,res);
544: Clear(q);
545: end;
546: Random_Product_System.Clear;
547: end if;
548: return res;
549: end Bezout_BKK;
550:
551: procedure Bezout_BKK ( k,n : in natural; p : in Poly_Sys;
552: tv : in out Tree_of_Vectors; bb : out natural ) is
553: res : natural;
554:
555: begin
556: if k = 0
557: then declare
558: adl : Array_of_Lists(p'range) := Create(p);
559: begin
560: Mixed_Volume(n,adl,tv,res);
561: Deep_Clear(adl);
562: end;
563: elsif k = n
564: then if Set_Structure.Empty
565: then Random_Product_Start_Systems.Build_Set_Structure(p);
566: res := Set_Structure.B;
567: Set_Structure.Clear;
568: else res := Set_Structure.B;
569: end if;
570: else Build_RPS(k,n,p);
571: declare
572: acc : Standard_Natural_Vectors.Vector(1..k);
573: q : Poly_Sys(1..n-k);
574: begin
575: acc := (1..k => 0); res := 0;
576: for i in q'range loop
577: q(i) := Complex_Randomize1(p(i+k));
578: end loop;
579: Build_Indices2(1,k,n,q,tv,acc,res);
580: Clear(q);
581: end;
582: Random_Product_System.Clear;
583: end if;
584: bb := res;
585: end Bezout_BKK;
586:
587: procedure Mixed_Solve
588: ( file : in file_type; k,n : in natural; p : in Poly_Sys;
589: bb : out natural;
590: g : in out Poly_Sys; sols : in out Solution_List ) is
591: begin
592: if k = 0
593: then declare
594: l : Laur_Sys(g'range);
595: tmp : Solution_List;
596: begin
597: l := Polynomial_to_Laurent_System(g);
598: Solve(file,l,bb,sols);
599: tmp := sols;
600: while not Is_Null(tmp) loop
601: Head_Of(tmp).t := Create(0.0);
602: tmp := Tail_Of(tmp);
603: end loop;
604: Clear(l);
605: end;
606: elsif k = n
607: then Random_Product_Start_Systems.Construct(p,g,sols);
608: bb := Length_Of(sols);
609: else Build_RPS(k,n,p);
610: declare
611: acc : Standard_Natural_Vectors.Vector(1..k);
612: q : Poly_Sys(1..n-k);
613: rq : Poly_Sys(p'range);
614: sols_last : Solution_List;
615: res : natural;
616: begin
617: acc := (1..k => 0); res := 0;
618: for i in q'range loop
619: q(i) := g(i+k);
620: end loop;
621: Build_Indices(file,q,1,k,n,acc,res,sols,sols_last);
622: bb := res;
623: rq := Random_Product_System.Polynomial_System;
624: g(1..k) := rq(1..k);
625: end;
626: Random_Product_System.Clear;
627: end if;
628: end Mixed_Solve;
629:
630: procedure Mixed_Solve
631: ( file : in file_type; k,n : in natural; p : in Poly_Sys;
632: tv : in Tree_of_Vectors; bb : out natural;
633: g : in out Poly_Sys; sols : in out Solution_List ) is
634: begin
635: if k = 0
636: then declare
637: l : Laur_Sys(g'range);
638: tmp : Solution_List;
639: begin
640: l := Polynomial_to_Laurent_System(g);
641: Solve(file,l,tv,bb,sols);
642: tmp := sols;
643: while not Is_Null(tmp) loop
644: Head_Of(tmp).t := Create(0.0);
645: tmp := Tail_Of(tmp);
646: end loop;
647: Clear(l);
648: end;
649: elsif k = n
650: then Random_Product_Start_Systems.Construct(p,g,sols);
651: bb := Length_Of(sols);
652: else Build_RPS(k,n,p);
653: declare
654: acc : Standard_Natural_Vectors.Vector(1..k);
655: q : Poly_Sys(1..n-k);
656: rq : Poly_Sys(p'range);
657: sols_last : Solution_List;
658: res : natural;
659: begin
660: acc := (1..k => 0); res := 0;
661: for i in q'range loop
662: q(i) := g(i+k);
663: end loop;
664: Build_Indices(file,q,1,k,n,acc,tv,res,sols,sols_last);
665: bb := res;
666: rq := Random_Product_System.Polynomial_System;
667: g(1..k) := rq(1..k);
668: end;
669: Random_Product_System.Clear;
670: end if;
671: end Mixed_Solve;
672:
673: end Set_Structures_and_Volumes;