Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Stalift/floating_polyhedral_continuation.adb, Revision 1.1.1.1
1.1 maekawa 1: with integer_io; use integer_io;
2: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
3: with Standard_Floating_Numbers_io; use Standard_Floating_Numbers_io;
4: with Standard_Mathematical_Functions; use Standard_Mathematical_Functions;
5: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
6: with Standard_Floating_Vectors_io; use Standard_Floating_Vectors_io;
7: with Standard_Complex_Vectors;
8: with Standard_Complex_Vectors_io; use Standard_Complex_Vectors_io;
9: with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
10: with Standard_Integer_VecVecs;
11: with Standard_Floating_VecVecs;
12: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
13: with Standard_Complex_Polynomials;
14: with Standard_Complex_Laur_Polys;
15: with Standard_Complex_Laur_Functions; use Standard_Complex_Laur_Functions;
16: with Standard_Complex_Poly_Systems; use Standard_Complex_Poly_Systems;
17: with Standard_Laur_Poly_Convertors; use Standard_Laur_Poly_Convertors;
18: with Standard_Complex_Substitutors; use Standard_Complex_Substitutors;
19: with Power_Lists; use Power_Lists;
20: with Lists_of_Floating_Vectors; use Lists_of_Floating_Vectors;
21: with Floating_Lifting_Utilities; use Floating_Lifting_Utilities;
22: with Floating_Integer_Convertors; use Floating_Integer_Convertors;
23: with Transforming_Integer_Vector_Lists; use Transforming_Integer_Vector_Lists;
24: with Transforming_Laurent_Systems; use Transforming_Laurent_Systems;
25: with Continuation_Parameters;
26: with Increment_and_Fix_Continuation; use Increment_and_Fix_Continuation;
27: with Fewnomial_System_Solvers; use Fewnomial_System_Solvers;
28: with BKK_Bound_Computations; use BKK_Bound_Computations;
29:
30: package body Floating_Polyhedral_Continuation is
31:
32: -- UTILITIES FOR POLYHEDRAL COEFFICIENT HOMOTOPY :
33:
34: procedure put ( file : in file_type;
35: v : in Standard_Floating_Vectors.Vector;
36: fore,aft,exp : in natural ) is
37: begin
38: for i in v'range loop
39: put(file,v(i),fore,aft,exp);
40: end loop;
41: end put;
42:
43: procedure put ( file : in file_type;
44: v : in Standard_Complex_Vectors.Vector;
45: fore,aft,exp : in natural ) is
46: begin
47: for i in v'range loop
48: put(file,REAL_PART(v(i)),fore,aft,exp);
49: put(file,IMAG_PART(v(i)),fore,aft,exp);
50: end loop;
51: end put;
52:
53: function Power ( x,m : double_float ) return double_float is
54:
55: -- DESCRIPTION :
56: -- Computes x**m for high powers of m to avoid RANGE_ERROR.
57:
58: -- intm : integer := integer(m);
59: -- fltm : double_float;
60:
61: begin
62: -- if m < 10.0
63: -- then
64: return (x**m); -- GNAT is better at this
65: -- else if double_float(intm) > m
66: -- then intm := intm-1;
67: -- end if;
68: -- fltm := m - double_float(intm);
69: -- return ((x**intm)*(x**fltm));
70: -- end if;
71: end Power;
72:
73: function Minimum ( v : Standard_Floating_Vectors.Vector )
74: return double_float is
75:
76: -- DESCRIPTION :
77: -- Returns the minimal element (/= 0) in the vector v.
78:
79: tol : constant double_float := 10.0**(-7);
80: min : double_float := 0.0;
81: tmp : double_float;
82:
83: begin
84: for i in v'range loop
85: tmp := abs(v(i));
86: if tmp > tol
87: then if (min = 0.0) or else (tmp < min)
88: then min := tmp;
89: end if;
90: end if;
91: end loop;
92: return min;
93: end Minimum;
94:
95: function Minimum ( v : Standard_Floating_VecVecs.VecVec )
96: return double_float is
97:
98: -- DESCRIPTION :
99: -- Returns the minimal nonzero element in the vectors of v.
100:
101: tol : constant double_float := 10.0**(-7);
102: min : double_float := 0.0;
103: tmp : double_float;
104:
105: begin
106: for i in v'range loop
107: tmp := Minimum(v(i).all);
108: if abs(tmp) > tol
109: then if (min = 0.0) or else (tmp < min)
110: then min := tmp;
111: end if;
112: end if;
113: end loop;
114: return min;
115: end Minimum;
116:
117: function Scale ( v : Standard_Floating_Vectors.Vector; s : double_float )
118: return Standard_Floating_Vectors.Vector is
119:
120: -- DESCRIPTION :
121: -- Returns the scaled vector, by dividing every element by s.
122:
123: res : Standard_Floating_Vectors.Vector(v'range);
124:
125: begin
126: for i in res'range loop
127: res(i) := v(i)/s;
128: end loop;
129: return res;
130: end Scale;
131:
132: function Scale ( v : Standard_Floating_VecVecs.VecVec )
133: return Standard_Floating_VecVecs.VecVec is
134:
135: -- DESCRIPTION :
136: -- Returns an array of scaled vectors.
137:
138: res : Standard_Floating_VecVecs.VecVec(v'range);
139: min : constant double_float := Minimum(v);
140: tol : constant double_float := 10.0**(-8);
141:
142: begin
143: if (abs(min) > tol) and (min /= 1.0)
144: then
145: for i in v'range loop
146: res(i) := new Standard_Floating_Vectors.Vector'(Scale(v(i).all,min));
147: end loop;
148: else
149: for i in v'range loop
150: res(i) := new Standard_Floating_Vectors.Vector'(v(i).all);
151: end loop;
152: end if;
153: return res;
154: end Scale;
155:
156: procedure Shift ( v : in out Standard_Floating_Vectors.Vector ) is
157:
158: -- DESCRIPTION :
159: -- Shifts the elements in v, such that the minimal element equals zero.
160:
161: min : double_float := v(v'first);
162:
163: begin
164: for i in v'first+1..v'last loop
165: if v(i) < min
166: then min := v(i);
167: end if;
168: end loop;
169: if min /= 0.0
170: then for i in v'range loop
171: v(i) := v(i) - min;
172: end loop;
173: end if;
174: end Shift;
175:
176: function Create ( e : Standard_Integer_VecVecs.VecVec;
177: l : List; normal : Standard_Floating_Vectors.Vector )
178: return Standard_Floating_Vectors.Vector is
179:
180: -- DESCRIPTION :
181: -- Returns a vector with all inner products of the normal with
182: -- the exponents in the list, such that minimal value equals zero.
183:
184: res : Standard_Floating_Vectors.Vector(e'range);
185: use Standard_Floating_Vectors;
186: found : boolean;
187: lif : double_float;
188:
189: begin
190: for i in e'range loop
191: declare
192: fei : constant Standard_Floating_Vectors.Vector := Convert(e(i).all);
193: begin
194: Search_Lifting(l,fei,found,lif);
195: if not found
196: then res(i) := 1000.0;
197: else res(i) := fei*normal(fei'range) + lif*normal(normal'last);
198: end if;
199: end;
200: end loop;
201: Shift(res);
202: return res;
203: end Create;
204:
205: function Create ( e : Exponent_Vectors_Array;
206: l : Array_of_Lists; mix : Standard_Integer_Vectors.Vector;
207: normal : Standard_Floating_Vectors.Vector )
208: return Standard_Floating_VecVecs.VecVec is
209:
210: res : Standard_Floating_VecVecs.VecVec(normal'first..normal'last-1);
211: cnt : natural := res'first;
212:
213: begin
214: for i in mix'range loop
215: declare
216: rescnt : constant Standard_Floating_Vectors.Vector
217: := Create(e(cnt).all,l(i),normal);
218: begin
219: res(cnt) := new Standard_Floating_Vectors.Vector(rescnt'range);
220: for j in rescnt'range loop
221: res(cnt)(j) := rescnt(j);
222: end loop;
223: end;
224: Shift(res(cnt).all);
225: for k in 1..(mix(i)-1) loop
226: res(cnt+k) := new Standard_Floating_Vectors.Vector(res(cnt)'range);
227: for j in res(cnt)'range loop
228: res(cnt+k)(j) := res(cnt)(j);
229: end loop;
230: end loop;
231: cnt := cnt + mix(i);
232: end loop;
233: return res;
234: end Create;
235:
236: procedure Eval ( c : in Standard_Complex_Vectors.Vector;
237: t : in double_float; m : in Standard_Floating_Vectors.Vector;
238: ctm : in out Standard_Complex_Vectors.Vector ) is
239:
240: -- DESCRIPTION : ctm = c*t**m.
241:
242: begin
243: for i in ctm'range loop
244: -- ctm(i) := c(i)*Create((t**m(i)));
245: ctm(i) := c(i)*Create(Power(t,m(i)));
246: end loop;
247: end Eval;
248:
249: procedure Eval ( c : in Standard_Complex_VecVecs.VecVec;
250: t : in double_float; m : in Standard_Floating_VecVecs.VecVec;
251: ctm : in out Standard_Complex_VecVecs.VecVec ) is
252:
253: -- DESCRIPTION : ctm = c*t**m.
254:
255: begin
256: for i in ctm'range loop
257: Eval(c(i).all,t,m(i).all,ctm(i).all);
258: end loop;
259: end Eval;
260:
261: -- USEFUL PRIMITIVES :
262:
263: function Is_In ( l : List;
264: d : Standard_Complex_Laur_Polys.Degrees )
265: return boolean is
266:
267: -- DESCRIPTION :
268: -- Returns true if the degrees belong to the list l.
269:
270: tmp : List := l;
271: pt : Standard_Floating_Vectors.Link_to_Vector;
272: fld : Standard_Floating_Vectors.Vector(d'range);
273: equ : boolean;
274:
275: begin
276: for i in fld'range loop
277: fld(i) := double_float(d(i));
278: end loop;
279: while not Is_Null(tmp) loop
280: pt := Head_Of(tmp);
281: equ := true;
282: for i in fld'range loop
283: if pt(i) /= fld(i)
284: then equ := false;
285: end if;
286: exit when not equ;
287: end loop;
288: if equ
289: then return true;
290: else tmp := Tail_Of(tmp);
291: end if;
292: end loop;
293: return false;
294: end Is_In;
295:
296: function Select_Terms ( p : Standard_Complex_Laur_Polys.Poly; l : List )
297: return Standard_Complex_Laur_Polys.Poly is
298:
299: use Standard_Complex_Laur_Polys;
300: res : Poly := Null_Poly;
301:
302: procedure Visit_Term ( t : in Term; cont : out boolean ) is
303: begin
304: if Is_In(l,t.dg)
305: then Add(res,t);
306: end if;
307: cont := true;
308: end Visit_Term;
309: procedure Visit_Terms is new Visiting_Iterator(Visit_Term);
310:
311: begin
312: Visit_Terms(p);
313: return res;
314: end Select_Terms;
315:
316: function Select_Subsystem
317: ( p : Laur_Sys; mix : Standard_Integer_Vectors.Vector;
318: mic : Mixed_Cell ) return Laur_Sys is
319:
320: -- DESCRIPTION :
321: -- Given a Laurent polynomial system and a mixed cell,
322: -- the corresponding subsystem will be returned.
323:
324: -- ON ENTRY :
325: -- p a Laurent polynomial system;
326: -- mix type of mixture: occurencies of the supports;
327: -- mic a mixed cell.
328:
329: -- REQUIRED :
330: -- The polynomials in p must be ordered according to the type of mixture.
331:
332: res : Laur_Sys(p'range);
333: cnt : natural := 0;
334:
335: begin
336: for k in mix'range loop
337: for l in 1..mix(k) loop
338: cnt := cnt + 1;
339: res(cnt) := Select_Terms(p(cnt),mic.pts(k));
340: end loop;
341: end loop;
342: return res;
343: end Select_Subsystem;
344:
345: procedure Extract_Regular ( sols : in out Solution_List ) is
346:
347: function To_Be_Removed ( flag : in natural ) return boolean is
348: begin
349: return ( flag /= 1 );
350: end To_Be_Removed;
351: procedure Extract_Regular_Solutions is new
352: Standard_Complex_Solutions.Delete(To_Be_Removed);
353:
354: begin
355: Extract_Regular_Solutions(sols);
356: end Extract_Regular;
357:
358: procedure Refine ( file : in file_type; p : in Laur_Sys;
359: sols : in out Solution_List ) is
360:
361: -- DESCRIPTION :
362: -- Given a polyhedral homotopy p and a list of solution for t=1,
363: -- this list of solutions will be refined.
364:
365: pp : Poly_Sys(p'range) := Laurent_to_Polynomial_System(p);
366: n : constant natural := p'length;
367: eps : constant double_float := 10.0**(-12);
368: tolsing : constant double_float := 10.0**(-8);
369: max : constant natural := 3;
370: numb : natural := 0;
371:
372: begin
373: pp := Laurent_to_Polynomial_System(p);
374: Substitute(n+1,Create(1.0),pp);
375: -- Reporting_Root_Refiner(file,pp,sols,eps,eps,tolsing,numb,max,false);
376: Clear(pp); Extract_Regular(sols);
377: end Refine;
378:
379: -- FIRST LAYER OF CONTINUATION ROUTINES :
380:
381: procedure Mixed_Continuation
382: ( mix : in Standard_Integer_Vectors.Vector;
383: lifted : in Array_of_Lists; h : in Eval_Coeff_Laur_Sys;
384: c : in Standard_Complex_VecVecs.VecVec;
385: e : in Exponent_Vectors_Array;
386: j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
387: normal : in Standard_Floating_Vectors.Vector;
388: sols : in out Solution_List ) is
389:
390: pow : Standard_Floating_VecVecs.VecVec(c'range)
391: := Create(e,lifted,mix,normal);
392: scapow : Standard_Floating_VecVecs.VecVec(c'range) := Scale(pow);
393: ctm : Standard_Complex_VecVecs.VecVec(c'range);
394:
395: use Standard_Floating_Vectors;
396: use Standard_Complex_Laur_Polys;
397:
398: function Eval ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
399: return Standard_Complex_Vectors.Vector is
400: begin
401: Eval(c,REAL_PART(t),scapow,ctm);
402: return Eval(h,ctm,x);
403: end Eval;
404:
405: function dHt ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
406: return Standard_Complex_Vectors.Vector is
407:
408: res : Standard_Complex_Vectors.Vector(h'range);
409: xtl : constant integer := x'last+1;
410:
411: begin
412: Eval(c,REAL_PART(t),scapow,ctm);
413: for i in res'range loop
414: res(i) := Eval(j(i,xtl),m(i,xtl).all,ctm(i).all,x);
415: end loop;
416: return res;
417: end dHt;
418:
419: function dHx ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
420: return matrix is
421:
422: mt : Matrix(x'range,x'range);
423:
424: begin
425: Eval(c,REAL_PART(t),scapow,ctm);
426: for k in mt'range(1) loop
427: for l in mt'range(2) loop
428: mt(k,l) := Eval(j(k,l),m(k,l).all,ctm(k).all,x);
429: end loop;
430: end loop;
431: return mt;
432: end dHx;
433:
434: procedure Laur_Cont is new Silent_Continue(Max_Norm,Eval,dHt,dHx);
435:
436: begin
437: for i in c'range loop
438: ctm(i)
439: := new Standard_Complex_Vectors.Vector'(c(i).all'range => Create(0.0));
440: end loop;
441: Laur_Cont(sols,false);
442: Standard_Floating_VecVecs.Clear(pow);
443: Standard_Floating_VecVecs.Clear(scapow);
444: Standard_Complex_VecVecs.Clear(ctm);
445: Extract_Regular(sols);
446: end Mixed_Continuation;
447:
448: procedure Mixed_Continuation
449: ( file : in file_type;
450: mix : in Standard_Integer_Vectors.Vector;
451: lifted : in Array_of_Lists; h : in Eval_Coeff_Laur_Sys;
452: c : in Standard_Complex_VecVecs.VecVec;
453: e : in Exponent_Vectors_Array;
454: j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
455: normal : in Standard_Floating_Vectors.Vector;
456: sols : in out Solution_List ) is
457:
458: pow : Standard_Floating_VecVecs.VecVec(c'range)
459: := Create(e,lifted,mix,normal);
460: scapow : Standard_Floating_VecVecs.VecVec(c'range) := Scale(pow);
461: ctm : Standard_Complex_VecVecs.VecVec(c'range);
462:
463: use Standard_Floating_Vectors;
464: use Standard_Complex_Laur_Polys;
465:
466: function Eval ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
467: return Standard_Complex_Vectors.Vector is
468: begin
469: Eval(c,REAL_PART(t),scapow,ctm);
470: return Eval(h,ctm,x);
471: end Eval;
472:
473: function dHt ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
474: return Standard_Complex_Vectors.Vector is
475:
476: res : Standard_Complex_Vectors.Vector(h'range);
477: xtl : constant integer := x'last+1;
478:
479: begin
480: Eval(c,REAL_PART(t),scapow,ctm);
481: for i in res'range loop
482: res(i) := Eval(j(i,xtl),m(i,xtl).all,ctm(i).all,x);
483: end loop;
484: return res;
485: end dHt;
486:
487: function dHx ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
488: return Matrix is
489:
490: mt : Matrix(x'range,x'range);
491:
492: begin
493: Eval(c,REAL_PART(t),scapow,ctm);
494: for k in m'range(1) loop
495: for l in m'range(1) loop
496: mt(k,l) := Eval(j(k,l),m(k,l).all,ctm(k).all,x);
497: end loop;
498: end loop;
499: return mt;
500: end dHx;
501:
502: procedure Laur_Cont is new Reporting_Continue(Max_Norm,Eval,dHt,dHx);
503:
504: begin
505: -- put_line(file,"The coefficient vectors :" );
506: -- for i in c'range loop
507: -- put(file,c(i).all,3,3,3); new_line(file);
508: -- end loop;
509: for i in c'range loop
510: ctm(i)
511: := new Standard_Complex_Vectors.Vector'(c(i).all'range => Create(0.0));
512: end loop;
513: -- put(file,"The normal : "); put(file,normal,3,3,3); new_line(file);
514: -- put_line(file,"The exponent vector : ");
515: -- for i in pow'range loop
516: -- put(file,pow(i).all,3,3,3); new_line(file);
517: -- end loop;
518: -- put_line(file,"The scaled exponent vector : ");
519: -- for i in pow'range loop
520: -- put(file,scapow(i).all,3,3,3); new_line(file);
521: -- end loop;
522: Laur_Cont(file,sols,false);
523: Standard_Floating_VecVecs.Clear(pow);
524: Standard_Floating_VecVecs.Clear(scapow);
525: Standard_Complex_VecVecs.Clear(ctm);
526: Extract_Regular(sols);
527: end Mixed_Continuation;
528:
529: -- UTILITIES FOR SECOND LAYER :
530:
531: function Remove_Lifting ( l : List ) return List is
532:
533: -- DESCRIPTION :
534: -- Removes the lifting value from the vectors.
535:
536: tmp,res,res_last : List;
537:
538: begin
539: tmp := l;
540: while not Is_Null(tmp) loop
541: declare
542: d1 : constant Standard_Floating_Vectors.Vector := Head_Of(tmp).all;
543: d2 : constant Standard_Floating_Vectors.Vector
544: := d1(d1'first..d1'last-1);
545: begin
546: Append(res,res_last,d2);
547: end;
548: tmp := Tail_Of(tmp);
549: end loop;
550: return res;
551: end Remove_Lifting;
552:
553: function Sub_Lifting ( q : Laur_Sys; mix : Standard_Integer_Vectors.Vector;
554: mic : Mixed_Cell ) return Array_of_Lists is
555:
556: -- DESCRIPTION :
557: -- Returns the lifting used to subdivide the cell.
558:
559: res : Array_of_Lists(mix'range);
560: sup : Array_of_Lists(q'range);
561: n : constant natural := q'last;
562: cnt : natural := sup'first;
563:
564: begin
565: for i in mic.pts'range loop
566: sup(cnt) := Remove_Lifting(mic.pts(i));
567: for j in 1..(mix(i)-1) loop
568: Copy(sup(cnt),sup(cnt+j));
569: end loop;
570: cnt := cnt + mix(i);
571: end loop;
572: res := Induced_Lifting(n,mix,sup,mic.sub.all);
573: Deep_Clear(sup);
574: return res;
575: end Sub_Lifting;
576:
577: function Sub_Polyhedral_Homotopy
578: ( l : List; e : Standard_Integer_VecVecs.VecVec;
579: c : Standard_Complex_Vectors.Vector )
580: return Standard_Complex_Vectors.Vector is
581:
582: -- DESCRIPTION :
583: -- For every vector in e that does not belong to l, the corresponding
584: -- index in c will be set to zero, otherwise it is copied to the result.
585:
586: res : Standard_Complex_Vectors.Vector(c'range);
587: found : boolean;
588: lif : double_float;
589:
590: begin
591: for i in e'range loop
592: declare
593: fei : constant Standard_Floating_Vectors.Vector := Convert(e(i).all);
594: begin
595: Search_Lifting(l,fei,found,lif);
596: if not found
597: then res(i) := Create(0.0);
598: else res(i) := c(i);
599: end if;
600: end;
601: end loop;
602: return res;
603: end Sub_Polyhedral_Homotopy;
604:
605: function Sub_Polyhedral_Homotopy
606: ( mix : Standard_Integer_Vectors.Vector; mic : Mixed_Cell;
607: e : Exponent_Vectors_Array;
608: c : Standard_Complex_VecVecs.VecVec )
609: return Standard_Complex_VecVecs.VecVec is
610:
611: -- DESCRIPTION :
612: -- Given a subsystem q of p and the coefficient vector of p, the
613: -- vector on return will have only nonzero entries for coefficients
614: -- that belong to q.
615:
616: res : Standard_Complex_VecVecs.VecVec(c'range);
617:
618: begin
619: for i in mix'range loop
620: declare
621: cri : constant Standard_Complex_Vectors.Vector
622: := Sub_Polyhedral_Homotopy(mic.pts(i),e(i).all,c(i).all);
623: begin
624: res(i) := new Standard_Complex_Vectors.Vector'(cri);
625: for j in 1..(mix(i)-1) loop
626: declare
627: crj : constant Standard_Complex_Vectors.Vector
628: := Sub_Polyhedral_Homotopy(mic.pts(i),e(i+j).all,c(i+j).all);
629: begin
630: res(i+j) := new Standard_Complex_Vectors.Vector'(crj);
631: end;
632: end loop;
633: end;
634: end loop;
635: return res;
636: end Sub_Polyhedral_Homotopy;
637:
638: procedure Refined_Mixed_Solve
639: ( q : in Laur_Sys; mix : in Standard_Integer_Vectors.Vector;
640: mic : in Mixed_Cell; h : in Eval_Coeff_Laur_Sys;
641: c : in Standard_Complex_VecVecs.VecVec;
642: e : in Exponent_Vectors_Array;
643: j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
644: qsols : in out Solution_List ) is
645:
646: -- DESCRIPTION :
647: -- Polyhedral coeffient-homotopy for subsystem q.
648:
649: -- REQUIRED : mic.sub /= null.
650:
651: lif : Array_of_Lists(mix'range) := Sub_Lifting(q,mix,mic);
652: cq : Standard_Complex_VecVecs.VecVec(c'range)
653: := Sub_Polyhedral_Homotopy(mix,mic,e,c);
654:
655: begin
656: Mixed_Solve(q,lif,h,cq,e,j,m,mix,mic.sub.all,qsols);
657: Standard_Complex_VecVecs.Clear(cq); Deep_Clear(lif);
658: end Refined_Mixed_Solve;
659:
660: procedure Refined_Mixed_Solve
661: ( file : in file_type; q : in Laur_Sys;
662: mix : in Standard_Integer_Vectors.Vector;
663: mic : in Mixed_Cell; h : in Eval_Coeff_Laur_Sys;
664: c : in Standard_Complex_VecVecs.VecVec;
665: e : in Exponent_Vectors_Array;
666: j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
667: qsols : in out Solution_List ) is
668:
669: -- DESCRIPTION :
670: -- Polyhedral coeffient-homotopy for subsystem q.
671:
672: -- REQUIRED : mic.sub /= null.
673:
674: lif : Array_of_Lists(mix'range) := Sub_Lifting(q,mix,mic);
675: cq : Standard_Complex_VecVecs.VecVec(c'range)
676: := Sub_Polyhedral_Homotopy(mix,mic,e,c);
677:
678: begin
679: Mixed_Solve(file,q,lif,h,cq,e,j,m,mix,mic.sub.all,qsols);
680: Standard_Complex_VecVecs.Clear(cq); Deep_Clear(lif);
681: end Refined_Mixed_Solve;
682:
683: -- SECOND LAYER :
684:
685: procedure Mixed_Solve
686: ( p : in Laur_Sys; lifted : in Array_of_Lists;
687: h : in Eval_Coeff_Laur_Sys;
688: c : in Standard_Complex_VecVecs.VecVec;
689: e : in Exponent_Vectors_Array;
690: j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
691: mix : in Standard_Integer_Vectors.Vector;
692: mic : in Mixed_Cell;
693: sols,sols_last : in out Solution_List ) is
694:
695: q : Laur_Sys(p'range) := Select_Subsystem(p,mix,mic);
696: sq : Laur_Sys(q'range) := Shift(q);
697: pq : Poly_Sys(q'range);
698: qsols : Solution_List;
699: len : natural := 0;
700: fail : boolean;
701:
702: begin
703: Solve(sq,qsols,fail);
704: if fail
705: then if mic.sub = null
706: then pq := Laurent_to_Polynomial_System(sq);
707: qsols := Solve_by_Static_Lifting(pq);
708: Clear(pq);
709: else Refined_Mixed_Solve(q,mix,mic,h,c,e,j,m,qsols);
710: end if;
711: Set_Continuation_Parameter(qsols,Create(0.0));
712: end if;
713: len := Length_Of(qsols);
714: if len > 0
715: then Mixed_Continuation(mix,lifted,h,c,e,j,m,mic.nor.all,qsols);
716: Concat(sols,sols_last,qsols);
717: end if;
718: Clear(sq); Clear(q); Clear(qsols);
719: end Mixed_Solve;
720:
721: procedure Mixed_Solve
722: ( file : in file_type;
723: p : in Laur_Sys; lifted : in Array_of_Lists;
724: h : in Eval_Coeff_Laur_Sys;
725: c : in Standard_Complex_VecVecs.VecVec;
726: e : in Exponent_Vectors_Array;
727: j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
728: mix : in Standard_Integer_Vectors.Vector;
729: mic : in Mixed_Cell;
730: sols,sols_last : in out Solution_List ) is
731:
732: q : Laur_Sys(p'range) := Select_Subsystem(p,mix,mic);
733: sq : Laur_Sys(q'range) := Shift(q);
734: pq : Poly_Sys(q'range);
735: qsols : Solution_List;
736: len : natural := 0;
737: fail : boolean;
738:
739: begin
740: Solve(sq,qsols,fail);
741: if not fail
742: then put_line(file,"It is a fewnomial system.");
743: else put_line(file,"No fewnomial system.");
744: if mic.sub = null
745: then put_line(file,"Calling the black box solver.");
746: pq := Laurent_to_Polynomial_System(sq);
747: qsols := Solve_by_Static_Lifting(file,pq);
748: Clear(pq);
749: else put_line(file,"Using the refinement of the cell.");
750: Refined_Mixed_Solve(file,q,mix,mic,h,c,e,j,m,qsols);
751: end if;
752: Set_Continuation_Parameter(qsols,Create(0.0));
753: end if;
754: len := Length_Of(qsols);
755: put(file,len,1); put_line(file," solutions found.");
756: if len > 0
757: then
758: Mixed_Continuation(file,mix,lifted,h,c,e,j,m,mic.nor.all,qsols);
759: Concat(sols,sols_last,qsols);
760: end if;
761: Clear(sq); Clear(q); Clear(qsols);
762: end Mixed_Solve;
763:
764: -- THIRD LAYER :
765:
766: procedure Mixed_Solve
767: ( p : in Laur_Sys; lifted : in Array_of_Lists;
768: h : in Eval_Coeff_Laur_Sys;
769: c : in Standard_Complex_VecVecs.VecVec;
770: e : in Exponent_Vectors_Array;
771: j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
772: mix : in Standard_Integer_Vectors.Vector;
773: mixsub : in Mixed_Subdivision;
774: sols : in out Solution_List ) is
775:
776: tmp : Mixed_Subdivision := mixsub;
777: mic : Mixed_Cell;
778: sols_last : Solution_List;
779:
780: begin
781: while not Is_Null(tmp) loop
782: mic := Head_Of(tmp);
783: Mixed_Solve(p,lifted,h,c,e,j,m,mix,mic,sols,sols_last);
784: tmp := Tail_Of(tmp);
785: end loop;
786: end Mixed_Solve;
787:
788: procedure Mixed_Solve
789: ( file : in file_type;
790: p : in Laur_Sys; lifted : in Array_of_Lists;
791: h : in Eval_Coeff_Laur_Sys;
792: c : in Standard_Complex_VecVecs.VecVec;
793: e : in Exponent_Vectors_Array;
794: j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
795: mix : in Standard_Integer_Vectors.Vector;
796: mixsub : in Mixed_Subdivision;
797: sols : in out Solution_List ) is
798:
799: tmp : Mixed_Subdivision := mixsub;
800: mic : Mixed_Cell;
801: sols_last : Solution_List;
802: cnt : natural := 0;
803:
804: begin
805: while not Is_Null(tmp) loop
806: mic := Head_Of(tmp);
807: cnt := cnt + 1;
808: new_line(file);
809: put(file,"*** PROCESSING SUBSYSTEM "); put(file,cnt,1);
810: put_line(file," ***");
811: new_line(file);
812: Mixed_Solve(file,p,lifted,h,c,e,j,m,mix,mic,sols,sols_last);
813: tmp := Tail_Of(tmp);
814: end loop;
815: end Mixed_Solve;
816:
817: end Floating_Polyhedral_Continuation;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to floating_polyhedral_continuation.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 / Stalift