Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Stalift/integer_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_Complex_Numbers; use Standard_Complex_Numbers;
4: with Standard_Complex_Numbers_io; use Standard_Complex_Numbers_io;
5: with Standard_Mathematical_Functions; use Standard_Mathematical_Functions;
6: with Standard_Integer_Vectors_io; use Standard_Integer_Vectors_io;
7: with Standard_Integer_VecVecs;
8: with Standard_Integer_VecVecs_io; use Standard_Integer_VecVecs_io;
9: with Standard_Floating_Vectors;
10: with Standard_Floating_VecVecs;
11: with Standard_Complex_Vectors;
12: with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
13: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
14: with Standard_Complex_Polynomials;
15: with Standard_Complex_Laur_Polys;
16: with Standard_Complex_Laur_Functions; use Standard_Complex_Laur_Functions;
17: with Standard_Complex_Poly_Systems; use Standard_Complex_Poly_Systems;
18: with Standard_Laur_Poly_Convertors; use Standard_Laur_Poly_Convertors;
19: with Standard_Complex_Substitutors; use Standard_Complex_Substitutors;
20: with Lists_of_Integer_Vectors; use Lists_of_Integer_Vectors;
21: with Power_Lists; use Power_Lists;
22: with Integer_Lifting_Utilities; use Integer_Lifting_Utilities;
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 Integer_Polyhedral_Continuation is
31:
32: procedure Write ( file : in file_type;
33: c : in Standard_Complex_VecVecs.VecVec ) is
34: begin
35: for i in c'range loop
36: put(file,i,1); put(file," : ");
37: for j in c(i)'range loop
38: if REAL_PART(c(i)(j)) = 0.0
39: then put(file," 0");
40: else put(file,c(i)(j),2,3,0);
41: end if;
42: end loop;
43: new_line(file);
44: end loop;
45: end Write;
46:
47: -- UTILITIES FOR POLYHEDRAL COEFFICIENT HOMOTOPY :
48:
49: function Power ( x,m : double_float ) return double_float is
50:
51: -- DESCRIPTION :
52: -- Computes x**m for high powers of m to avoid RANGE_ERROR.
53:
54: -- intm : integer := integer(m);
55: -- fltm : double_float;
56:
57: begin
58: -- if m < 10.0
59: -- then
60: return (x**m); -- FOR GNAT 3.07
61: -- else if double_float(intm) > m
62: -- then intm := intm-1;
63: -- end if;
64: -- fltm := m - double_float(intm);
65: -- return ((x**intm)*(x**fltm));
66: -- end if;
67: end Power;
68:
69: function Minimum ( v : Standard_Integer_Vectors.Vector ) return integer is
70:
71: -- DESCRIPTION :
72: -- Returns the minimal element in the vector v, different from zero.
73:
74: tmp,min : integer := 0;
75:
76: begin
77: for i in v'range loop
78: if v(i) /= 0
79: then if v(i) < 0
80: then tmp := -v(i);
81: end if;
82: if min = 0
83: then min := tmp;
84: elsif tmp < min
85: then min := tmp;
86: end if;
87: end if;
88: end loop;
89: return min;
90: end Minimum;
91:
92: function Scale ( v : Standard_Integer_Vectors.Vector )
93: return Standard_Floating_Vectors.Vector is
94:
95: -- DESCRIPTION :
96: -- Returns the vector v divided by its minimal element different from zero,
97: -- such that the smallest positive element in the scaled vector equals one.
98:
99: res : Standard_Floating_Vectors.Vector(v'range);
100: min : constant integer := Minimum(v);
101:
102: begin
103: if (min = 0) or (min = 1)
104: then for i in res'range loop
105: res(i) := double_float(v(i));
106: end loop;
107: else for i in res'range loop
108: res(i) := double_float(v(i))/double_float(min);
109: end loop;
110: end if;
111: return res;
112: end Scale;
113:
114: function Scale ( v : Standard_Integer_VecVecs.VecVec )
115: return Standard_Floating_VecVecs.VecVec is
116:
117: -- DESCRIPTION :
118: -- Returns an array of scaled vectors.
119:
120: res : Standard_Floating_VecVecs.VecVec(v'range);
121:
122: begin
123: for i in v'range loop
124: declare
125: sv : constant Standard_Floating_Vectors.Vector := Scale(v(i).all);
126: begin
127: res(i) := new Standard_Floating_Vectors.Vector(sv'range);
128: for j in sv'range loop
129: res(i)(j) := sv(j);
130: end loop;
131: end; -- detour set up for GNAT 3.07
132: -- res(i) := new Standard_Floating_Vectors.Vector'(Scale(v(i).all));
133: end loop;
134: return res;
135: end Scale;
136:
137: procedure Shift ( v : in out Standard_Integer_Vectors.Vector ) is
138:
139: -- DESCRIPTION :
140: -- Shifts the elements in v, such that the minimal element equals zero.
141:
142: min : integer := v(v'first);
143:
144: begin
145: for i in v'first+1..v'last loop
146: if v(i) < min
147: then min := v(i);
148: end if;
149: end loop;
150: if min /= 0
151: then for i in v'range loop
152: v(i) := v(i) - min;
153: end loop;
154: end if;
155: end Shift;
156:
157: function Create ( e : Standard_Integer_VecVecs.VecVec;
158: l : List; normal : Vector )
159: return Standard_Integer_Vectors.Vector is
160:
161: -- DESCRIPTION :
162: -- Returns a vector with all inner products of the normal with
163: -- the exponents in the list, such that minimal value equals zero.
164:
165: res : Standard_Integer_Vectors.Vector(e'range);
166: lei : Standard_Integer_Vectors.Vector(normal'first..normal'last-1);
167: found : boolean;
168: lif : integer;
169:
170: begin
171: for i in e'range loop
172: lei := e(i).all;
173: Search_Lifting(l,lei,found,lif);
174: if not found
175: then res(i) := 1000;
176: else res(i) := lei*normal(lei'range) + normal(normal'last)*lif;
177: end if;
178: end loop;
179: Shift(res);
180: return res;
181: end Create;
182:
183: function Create ( e : Exponent_Vectors_Array;
184: l : Array_of_Lists; mix,normal : Vector )
185: return Standard_Integer_VecVecs.VecVec is
186:
187: res : Standard_Integer_VecVecs.VecVec(e'range);
188: cnt : natural := res'first;
189:
190: begin
191: for i in mix'range loop
192: declare
193: rescnt : constant Standard_Integer_Vectors.Vector
194: := Create(e(cnt).all,l(i),normal);
195: begin
196: res(cnt) := new Standard_Integer_Vectors.Vector(rescnt'range);
197: for j in rescnt'range loop
198: res(cnt)(j) := rescnt(j);
199: end loop;
200: end;
201: Shift(res(cnt).all);
202: for k in 1..(mix(i)-1) loop
203: res(cnt+k) := new Standard_Integer_Vectors.Vector(res(cnt)'range);
204: for j in res(cnt)'range loop
205: res(cnt+k)(j) := res(cnt)(j);
206: end loop;
207: end loop;
208: cnt := cnt + mix(i);
209: end loop;
210: return res;
211: end Create;
212:
213: function Chebychev_Poly
214: ( k : integer; t : double_float ) return double_float is
215:
216: -- DESCRIPTION : returns T_k(t), T_k = kth Chebychev polynomial.
217: -- note that k can also be negative...
218:
219: begin
220: return COS(double_float(k)*ARCCOS(t));
221: end Chebychev_Poly;
222:
223: function Chebychev_Interpolate
224: ( k : natural; t : double_float ) return double_float is
225:
226: -- DESCRIPTION :
227: -- Uses Chebychev polynomials to return a polynomial f of degree k,
228: -- evaluated at t, such that f(1) = 1 and f(0) = 0 (except if k=0).
229:
230: arg : double_float;
231:
232: begin
233: if k = 0
234: then return 1.0;
235: else arg := (1.0-t)*COS(PI/(2.0*double_float(k))) + t;
236: return Chebychev_Poly(k,arg);
237: end if;
238: end Chebychev_Interpolate;
239:
240: function Interpolate
241: ( k : natural; t : double_float ) return Complex_Number is
242:
243: -- DESCRIPTION :
244: -- Evaluates t^k on the complex unit circle.
245:
246: arg : double_float := 2.0*double_float(k)*PI*t;
247:
248: begin
249: return ((Create(1.0) - Create(SIN(arg),COS(arg)))/Create(2.0));
250: end Interpolate;
251:
252: procedure Eval ( c : in Standard_Complex_Vectors.Vector;
253: t : in double_float; m : in Standard_Integer_Vectors.Vector;
254: ctm : in out Standard_Complex_Vectors.Vector ) is
255:
256: -- DESCRIPTION : ctm = c*t**m.
257:
258: begin
259: for i in ctm'range loop
260: ctm(i) := c(i)*Create((t**m(i)));
261: -- ctm(i) := c(i)*Create(Chebychev_Interpolate(m(i),t));
262: -- ctm(i) := c(i)*Interpolate(m(i),t);
263: end loop;
264: end Eval;
265:
266: procedure Eval ( c : in Standard_Complex_Vectors.Vector;
267: t : in double_float; m : in Standard_Floating_Vectors.Vector;
268: ctm : in out Standard_Complex_Vectors.Vector ) is
269:
270: -- DESCRIPTION : ctm = c*t**m.
271:
272: begin
273: for i in ctm'range loop
274: -- ctm(i) := c(i)*Create((t**m(i)));
275: if (REAL_PART(c(i)) = 0.0) and (IMAG_PART(c(i)) = 0.0)
276: then ctm(i) := Create(0.0);
277: else ctm(i) := c(i)*Create(Power(t,m(i)));
278: end if;
279: end loop;
280: end Eval;
281:
282: procedure Eval ( c : in Standard_Complex_VecVecs.VecVec;
283: t : in double_float;
284: m : in Standard_Integer_VecVecs.VecVec;
285: ctm : in out Standard_Complex_VecVecs.VecVec ) is
286:
287: -- DESCRIPTION : ctm = c*t**m.
288:
289: begin
290: for i in ctm'range loop
291: Eval(c(i).all,t,m(i).all,ctm(i).all);
292: end loop;
293: end Eval;
294:
295: procedure Eval ( c : in Standard_Complex_VecVecs.VecVec;
296: t : in double_float; m : in Standard_Floating_VecVecs.VecVec;
297: ctm : in out Standard_Complex_VecVecs.VecVec ) is
298:
299: -- DESCRIPTION : ctm = c*t**m.
300:
301: begin
302: for i in ctm'range loop
303: Eval(c(i).all,t,m(i).all,ctm(i).all);
304: end loop;
305: end Eval;
306:
307: -- pragma inline(Eval);
308:
309: -- HOMOTOPY CONSTRUCTOR :
310:
311: function Select_Subsystem ( p : Laur_Sys; mix : Vector; mic : Mixed_Cell )
312: return Laur_Sys is
313:
314: res : Laur_Sys(p'range);
315: cnt : natural := 0;
316:
317: begin
318: for k in mix'range loop
319: for l in 1..mix(k) loop
320: cnt := cnt + 1;
321: res(cnt) := Select_Terms(p(cnt),mic.pts(k));
322: end loop;
323: end loop;
324: return res;
325: end Select_Subsystem;
326:
327: function Construct_Homotopy ( p : Laur_Sys; normal : Vector )
328: return Laur_Sys is
329:
330: -- DESCRIPTION :
331: -- Given a Laurent polynomial system of dimension n*(n+1) and a
332: -- normal, a homotopy will be constructed, with t = x(n+1)
333: -- and so that the support of the start system corresponds with
334: -- all points which give the minimal product with the normal.
335:
336: res : Laur_Sys(p'range);
337: n : constant natural := p'length;
338: use Standard_Complex_Laur_Polys;
339:
340: function Construct_Polynomial ( p : Poly; v : Vector ) return Poly is
341:
342: res : Poly := Null_Poly;
343:
344: procedure Construct_Term ( t : in Term; cont : out boolean ) is
345:
346: rt : term;
347:
348: begin
349: rt.cf := t.cf;
350: rt.dg := new Standard_Integer_Vectors.Vector'(t.dg.all);
351: rt.dg(n+1) := t.dg.all*v;
352: Add(res,rt);
353: Clear(rt);
354: cont := true;
355: end Construct_Term;
356: procedure Construct_Terms is
357: new Standard_Complex_Laur_Polys.Visiting_Iterator(Construct_Term);
358:
359: begin
360: Construct_Terms(p);
361: return res;
362: end Construct_Polynomial;
363:
364: begin
365: -- SUBSTITUTIONS :
366: for k in p'range loop
367: res(k) := Construct_Polynomial(p(k),normal);
368: end loop;
369: -- SHIFT :
370: for k in res'range loop
371: declare
372: d : integer := Minimal_Degree(res(k),n+1);
373: t : Term;
374: begin
375: t.cf := Create(1.0);
376: t.dg := new Standard_Integer_Vectors.Vector'(1..n+1 => 0);
377: t.dg(n+1) := -d;
378: Mul(res(k),t);
379: Clear(t);
380: end;
381: end loop;
382: return res;
383: end Construct_Homotopy;
384:
385: function Determine_Power ( n : natural; h : Laur_Sys ) return positive is
386:
387: -- DESCRIPTION :
388: -- Returns the smallest power of the last unknown,
389: -- over all polynomials in h.
390:
391: res : positive := 1;
392: first : boolean := true;
393: d : integer;
394: use Standard_Complex_Laur_Polys;
395:
396: procedure Scan ( t : in Term; cont : out boolean ) is
397: begin
398: if (t.dg(n+1) > 0) and then (t.dg(n+1) < d)
399: then d := t.dg(n+1);
400: end if;
401: cont := true;
402: end Scan;
403: procedure Search_Positive_Minimum is new Visiting_Iterator(Scan);
404:
405: begin
406: for i in h'range loop
407: d := Maximal_Degree(h(i),n+1);
408: if d > 0
409: then Search_Positive_Minimum(h(i));
410: if first
411: then res := d;
412: first := false;
413: elsif d < res
414: then res := d;
415: end if;
416: end if;
417: exit when (d=1);
418: end loop;
419: if res = 1
420: then return res;
421: else return 2;
422: end if;
423: end Determine_Power;
424:
425: procedure Extract_Regular ( sols : in out Solution_List ) is
426:
427: function To_Be_Removed ( flag : in natural ) return boolean is
428: begin
429: return ( flag /= 1 );
430: end To_Be_Removed;
431: procedure Extract_Regular_Solutions is
432: new Standard_Complex_Solutions.Delete(To_Be_Removed);
433:
434: begin
435: Extract_Regular_Solutions(sols);
436: end Extract_Regular;
437:
438: procedure Refine ( file : in file_type; p : in Laur_Sys;
439: sols : in out Solution_List ) is
440:
441: -- DESCRIPTION :
442: -- Given a polyhedral homotopy p and a list of solution for t=1,
443: -- this list of solutions will be refined.
444:
445: pp : Poly_Sys(p'range) := Laurent_to_Polynomial_System(p);
446: n : constant natural := p'length;
447: eps : constant double_float := 10.0**(-12);
448: tolsing : constant double_float := 10.0**(-8);
449: max : constant natural := 3;
450: numb : natural := 0;
451:
452: begin
453: pp := Laurent_to_Polynomial_System(p);
454: Substitute(n+1,Create(1.0),pp);
455: -- Reporting_Root_Refiner(file,pp,sols,eps,eps,tolsing,numb,max,false);
456: Clear(pp); Extract_Regular(sols);
457: end Refine;
458:
459: -- FIRST LAYER OF CONTINUATION ROUTINES :
460:
461: procedure Mixed_Continuation
462: ( p : in Laur_Sys; normal : in Vector;
463: sols : in out Solution_List ) is
464:
465: n : constant natural := p'length;
466:
467: h : Laur_Sys(p'range) := Construct_Homotopy(p,normal);
468: hpe : Eval_Laur_Sys(h'range) := Create(h);
469: j : Jaco_Mat(h'range,h'first..h'last+1) := Create(h);
470: je : Eval_Jaco_Mat(j'range(1),j'range(2)) := Create(j);
471:
472: use Standard_Complex_Laur_Polys;
473:
474: function Eval ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
475: return Standard_Complex_Vectors.Vector is
476:
477: xt : Standard_Complex_Vectors.Vector(x'first..x'last+1);
478:
479: begin
480: xt(x'range) := x;
481: xt(xt'last) := t;
482: return Eval(hpe,xt);
483: end Eval;
484:
485: function dHt ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
486: return Standard_Complex_Vectors.Vector is
487:
488: res : Standard_Complex_Vectors.Vector(p'range);
489: xt : Standard_Complex_Vectors.Vector(x'first..x'last+1);
490:
491: begin
492: xt(x'range) := x;
493: xt(xt'last) := t;
494: for i in res'range loop
495: res(i) := Eval(je(i,xt'last),xt);
496: end loop;
497: return res;
498: end dHt;
499:
500: function dHx ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
501: return matrix is
502:
503: m : Matrix(x'range,x'range);
504: xt : Standard_Complex_Vectors.Vector(x'first..x'last+1);
505:
506: begin
507: xt(x'range) := x;
508: xt(xt'last) := t;
509: for i in m'range(1) loop
510: for j in m'range(1) loop
511: m(i,j) := Eval(je(i,j),xt);
512: end loop;
513: end loop;
514: return m;
515: end dHx;
516:
517: -- pragma inline(Eval,dHt,dHx);
518:
519: procedure Laur_Cont is new Silent_Continue(Max_Norm,Eval,dHt,dHx);
520:
521: begin
522: -- Continuation_Parameters.power_of_t := Determine_Power(h'length,h);
523: Laur_Cont(sols,false);
524: Clear(h); Clear(hpe); Clear(j); Clear(je);
525: Extract_Regular(sols);
526: end Mixed_Continuation;
527:
528: procedure Mixed_Continuation
529: ( file : in file_type; p : in Laur_Sys;
530: normal : in Vector; sols : in out Solution_List ) is
531:
532: h : Laur_Sys(p'range) := Construct_Homotopy(p,normal);
533: hpe : Eval_Laur_Sys(h'range) := Create(h);
534: j : Jaco_Mat(h'range,h'first..h'last+1) := Create(h);
535: je : Eval_Jaco_Mat(j'range(1),j'range(2)) := Create(j);
536:
537: use Standard_Complex_Laur_Polys;
538:
539: function Eval ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
540: return Standard_Complex_Vectors.Vector is
541:
542: xt : Standard_Complex_Vectors.Vector(x'first..x'last+1);
543:
544: begin
545: xt(x'range) := x;
546: xt(xt'last) := t;
547: return Eval(hpe,xt);
548: end Eval;
549:
550: function dHt ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
551: return Standard_Complex_Vectors.Vector is
552:
553: res : Standard_Complex_Vectors.Vector(p'range);
554: xt : Standard_Complex_Vectors.Vector(x'first..x'last+1);
555:
556: begin
557: xt(x'range) := x;
558: xt(xt'last) := t;
559: for i in res'range loop
560: res(i) := Eval(je(i,xt'last),xt);
561: end loop;
562: return res;
563: end dHt;
564:
565: function dHx ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
566: return Matrix is
567:
568: m : Matrix(x'range,x'range);
569: xt : Standard_Complex_Vectors.Vector(x'first..x'last+1);
570:
571: begin
572: xt(x'range) := x;
573: xt(xt'last) := t;
574: for i in m'range(1) loop
575: for j in m'range(1) loop
576: m(i,j) := Eval(je(i,j),xt);
577: end loop;
578: end loop;
579: return m;
580: end dHx;
581:
582: -- pragma inline(Eval,dHt,dHx);
583:
584: procedure Laur_Cont is new Reporting_Continue(Max_Norm,Eval,dHt,dHx);
585:
586: begin
587: -- Continuation_Parameters.power_of_t := Determine_Power(h'length,h);
588: Laur_Cont(file,sols,false);
589: Clear(h); Clear(hpe); Clear(j); Clear(je);
590: Extract_Regular(sols);
591: end Mixed_Continuation;
592:
593: procedure Mixed_Continuation
594: ( mix : in Standard_Integer_Vectors.Vector;
595: lifted : in Array_of_Lists; h : in Eval_Coeff_Laur_Sys;
596: c : in Standard_Complex_VecVecs.VecVec;
597: e : in Exponent_Vectors_Array;
598: j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
599: normal : in Vector; sols : in out Solution_List ) is
600:
601: pow : Standard_Integer_VecVecs.VecVec(c'range)
602: := Create(e,lifted,mix,normal);
603: -- scapow : Standard_Floating_VecVecs.VecVec(c'range) := Scale(pow);
604: ctm : Standard_Complex_VecVecs.VecVec(c'range);
605:
606: use Standard_Complex_Laur_Polys;
607:
608: function Eval ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
609: return Standard_Complex_Vectors.Vector is
610: begin
611: -- Eval(c,REAL_PART(t),scapow,ctm);
612: Eval(c,REAL_PART(t),pow,ctm);
613: return Eval(h,ctm,x);
614: end Eval;
615:
616: function dHt ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
617: return Standard_Complex_Vectors.Vector is
618:
619: res : Standard_Complex_Vectors.Vector(h'range);
620: xtl : constant integer := x'last+1;
621:
622: begin
623: -- Eval(c,REAL_PART(t),scapow,ctm);
624: Eval(c,REAL_PART(t),pow,ctm);
625: for i in res'range loop
626: res(i) := Eval(j(i,xtl),m(i,xtl).all,ctm(i).all,x);
627: end loop;
628: return res;
629: end dHt;
630:
631: function dHx ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
632: return matrix is
633:
634: mt : Matrix(x'range,x'range);
635:
636: begin
637: -- Eval(c,REAL_PART(t),scapow,ctm);
638: Eval(c,REAL_PART(t),pow,ctm);
639: for k in mt'range(1) loop
640: for l in mt'range(2) loop
641: mt(k,l) := Eval(j(k,l),m(k,l).all,ctm(k).all,x);
642: end loop;
643: end loop;
644: return mt;
645: end dHx;
646:
647: -- pragma inline(Eval,dHt,dHx);
648:
649: procedure Laur_Cont is new Silent_Continue(Max_Norm,Eval,dHt,dHx);
650:
651: begin
652: for i in c'range loop
653: ctm(i) := new Standard_Complex_Vectors.Vector(c(i)'range);
654: end loop;
655: Laur_Cont(sols,false);
656: Standard_Integer_VecVecs.Clear(pow);
657: -- Standard_Floating_VecVecs.Clear(scapow);
658: Standard_Complex_VecVecs.Clear(ctm);
659: Extract_Regular(sols);
660: end Mixed_Continuation;
661:
662: procedure Mixed_Continuation
663: ( file : in file_type; mix : in Standard_Integer_Vectors.Vector;
664: lifted : in Array_of_Lists; h : in Eval_Coeff_Laur_Sys;
665: c : in Standard_Complex_VecVecs.VecVec;
666: e : in Exponent_Vectors_Array;
667: j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
668: normal : in Vector; sols : in out Solution_List ) is
669:
670: pow : Standard_Integer_VecVecs.VecVec(c'range)
671: := Create(e,lifted,mix,normal);
672: -- scapow : Standard_Floating_VecVecs.VecVec(c'range) := Scale(pow);
673: ctm : Standard_Complex_VecVecs.VecVec(c'range);
674:
675: use Standard_Complex_Laur_Polys;
676:
677: function Eval ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
678: return Standard_Complex_Vectors.Vector is
679: begin
680: -- Eval(c,REAL_PART(t),scapow,ctm);
681: Eval(c,REAL_PART(t),pow,ctm);
682: return Eval(h,ctm,x);
683: end Eval;
684:
685: function dHt ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
686: return Standard_Complex_Vectors.Vector is
687:
688: res : Standard_Complex_Vectors.Vector(h'range);
689: xtl : constant integer := x'last+1;
690:
691: begin
692: -- Eval(c,REAL_PART(t),scapow,ctm);
693: Eval(c,REAL_PART(t),pow,ctm);
694: for i in res'range loop
695: res(i) := Eval(j(i,xtl),m(i,xtl).all,ctm(i).all,x);
696: end loop;
697: return res;
698: end dHt;
699:
700: function dHx ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
701: return Matrix is
702:
703: mt : Matrix(x'range,x'range);
704:
705: begin
706: -- Eval(c,REAL_PART(t),scapow,ctm);
707: Eval(c,REAL_PART(t),pow,ctm);
708: for k in m'range(1) loop
709: for l in m'range(1) loop
710: mt(k,l) := Eval(j(k,l),m(k,l).all,ctm(k).all,x);
711: end loop;
712: end loop;
713: return mt;
714: end dHx;
715:
716: -- pragma inline(Eval,dHt,dHx);
717:
718: procedure Laur_Cont is new Reporting_Continue(Max_Norm,Eval,dHt,dHx);
719:
720: begin
721: -- put(file,"The normal : "); put(file,normal); new_line(file);
722: -- put_line(file,"The exponent vector : ");
723: -- for i in pow'range loop
724: -- put(file,pow(i)); new_line(file);
725: -- end loop;
726: -- put_line(file,"The coefficient vector : "); Write(file,c);
727: for i in c'range loop
728: ctm(i) := new Standard_Complex_Vectors.Vector(c(i)'range);
729: end loop;
730: Laur_Cont(file,sols,false);
731: Standard_Integer_VecVecs.Clear(pow);
732: -- Standard_Floating_VecVecs.Clear(scapow);
733: Standard_Complex_VecVecs.Clear(ctm);
734: Extract_Regular(sols);
735: end Mixed_Continuation;
736:
737: -- UTILITIES FOR SECOND LAYER :
738:
739: function Sub_Lifting ( q : Laur_Sys; mix : Vector; mic : Mixed_Cell )
740: return Array_of_Lists is
741:
742: -- DESCRIPTION :
743: -- Returns the lifting used to subdivide the cell.
744:
745: res : Array_of_Lists(mix'range);
746: sup : Array_of_Lists(q'range);
747: n : constant natural := q'last;
748: cnt : natural := sup'first;
749:
750: begin
751: for i in mic.pts'range loop
752: sup(cnt) := Reduce(mic.pts(i),q'last+1);
753: for j in 1..(mix(i)-1) loop
754: Copy(sup(cnt),sup(cnt+j));
755: end loop;
756: cnt := cnt + mix(i);
757: end loop;
758: res := Induced_Lifting(n,mix,sup,mic.sub.all);
759: Deep_Clear(sup);
760: return res;
761: end Sub_Lifting;
762:
763: procedure Refined_Mixed_Solve
764: ( q : in Laur_Sys; mix : in Vector; mic : in Mixed_Cell;
765: qsols : in out Solution_List ) is
766:
767: -- DESCRIPTION :
768: -- Solves a subsystem q using the refinement of the cell mic.
769:
770: -- REQUIRED : mic.sub /= null.
771:
772: lif : Array_of_Lists(mix'range) := Sub_Lifting(q,mix,mic);
773: lifq : Laur_Sys(q'range) := Perform_Lifting(q'last,mix,lif,q);
774:
775: begin
776: Mixed_Solve(lifq,mix,mic.sub.all,qsols);
777: Deep_Clear(lif); Clear(lifq);
778: end Refined_Mixed_Solve;
779:
780: procedure Refined_Mixed_Solve
781: ( file : in file_type;
782: q : in Laur_Sys; mix : in Vector; mic : in Mixed_Cell;
783: qsols : in out Solution_List ) is
784:
785: -- DESCRIPTION :
786: -- Solves a subsystem q using the refinement of the cell mic.
787:
788: -- REQUIRED : mic.sub /= null.
789:
790: lif : Array_of_Lists(mix'range) := Sub_Lifting(q,mix,mic);
791: lifq : Laur_Sys(q'range) := Perform_Lifting(q'last,mix,lif,q);
792:
793: begin
794: Mixed_Solve(file,lifq,mix,mic.sub.all,qsols);
795: Deep_Clear(lif); Clear(lifq);
796: end Refined_Mixed_Solve;
797:
798: function Sub_Polyhedral_Homotopy
799: ( l : List; e : Standard_Integer_VecVecs.VecVec;
800: c : Standard_Complex_Vectors.Vector )
801: return Standard_Complex_Vectors.Vector is
802:
803: -- DESCRIPTION :
804: -- For every vector in e that does not belong to l, the corresponding
805: -- index in c will be set to zero, otherwise it is copied to the result.
806:
807: res : Standard_Complex_Vectors.Vector(c'range);
808: found : boolean;
809: lif : integer;
810:
811: begin
812: for i in e'range loop
813: Search_Lifting(l,e(i).all,found,lif);
814: if not found
815: then res(i) := Create(0.0);
816: else res(i) := c(i);
817: end if;
818: end loop;
819: return res;
820: end Sub_Polyhedral_Homotopy;
821:
822: function Sub_Polyhedral_Homotopy
823: ( mix : Vector; mic : Mixed_Cell;
824: e : Exponent_Vectors_Array;
825: c : Standard_Complex_VecVecs.VecVec )
826: return Standard_Complex_VecVecs.VecVec is
827:
828: -- DESCRIPTION :
829: -- Given a subsystem q of p and the coefficient vector of p, the
830: -- vector on return will have only nonzero entries for coefficients
831: -- that belong to q.
832:
833: res : Standard_Complex_VecVecs.VecVec(c'range);
834: ind : natural := 0;
835:
836: begin
837: for i in mix'range loop
838: ind := ind+1;
839: declare
840: cri : constant Standard_Complex_Vectors.Vector
841: := Sub_Polyhedral_Homotopy(mic.pts(i),e(ind).all,c(ind).all);
842: begin
843: res(ind) := new Standard_Complex_Vectors.Vector'(cri);
844: for j in 1..(mix(i)-1) loop
845: declare
846: crj : constant Standard_Complex_Vectors.Vector
847: := Sub_Polyhedral_Homotopy(mic.pts(i),e(i+j).all,c(i+j).all);
848: begin
849: ind := ind+1;
850: res(ind) := new Standard_Complex_Vectors.Vector'(crj);
851: end;
852: end loop;
853: end;
854: end loop;
855: return res;
856: end Sub_Polyhedral_Homotopy;
857:
858: procedure Refined_Mixed_Solve
859: ( q : in Laur_Sys; mix : in Vector; mic : in Mixed_Cell;
860: h : in Eval_Coeff_Laur_Sys;
861: c : in Standard_Complex_VecVecs.VecVec;
862: e : in Exponent_Vectors_Array;
863: j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
864: qsols : in out Solution_List ) is
865:
866: -- DESCRIPTION :
867: -- Polyhedral coeffient-homotopy for subsystem q.
868:
869: -- REQUIRED : mic.sub /= null.
870:
871: lif : Array_of_Lists(mix'range) := Sub_Lifting(q,mix,mic);
872: lq : Laur_Sys(q'range) := Perform_Lifting(q'last,mix,lif,q);
873: cq : Standard_Complex_VecVecs.VecVec(c'range)
874: := Sub_Polyhedral_Homotopy(mix,mic,e,c);
875:
876: begin
877: Mixed_Solve(lq,lif,h,cq,e,j,m,mix,mic.sub.all,qsols);
878: Standard_Complex_VecVecs.Clear(cq); Deep_Clear(lif); Clear(lq);
879: end Refined_Mixed_Solve;
880:
881: procedure Refined_Mixed_Solve
882: ( file : in file_type; q : in Laur_Sys; mix : in Vector;
883: mic : in Mixed_Cell; h : in Eval_Coeff_Laur_Sys;
884: c : in Standard_Complex_VecVecs.VecVec;
885: e : in Exponent_Vectors_Array;
886: j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
887: qsols : in out Solution_List ) is
888:
889: -- DESCRIPTION :
890: -- Polyhedral coeffient-homotopy for subsystem q.
891:
892: -- REQUIRED : mic.sub /= null.
893:
894: lif : Array_of_Lists(mix'range) := Sub_Lifting(q,mix,mic);
895: lq : Laur_Sys(q'range) := Perform_Lifting(q'last,mix,lif,q);
896: cq : Standard_Complex_VecVecs.VecVec(c'range)
897: := Sub_Polyhedral_Homotopy(mix,mic,e,c);
898:
899: begin
900: Mixed_Solve(file,lq,lif,h,cq,e,j,m,mix,mic.sub.all,qsols);
901: Standard_Complex_VecVecs.Clear(cq); Deep_Clear(lif); Clear(lq);
902: end Refined_Mixed_Solve;
903:
904: -- TARGET ROUTINES FOR SECOND LAYER :
905:
906: procedure Mixed_Solve
907: ( p : in Laur_Sys; mix : in Vector; mic : in Mixed_Cell;
908: sols,sols_last : in out Solution_List ) is
909:
910: q : Laur_Sys(p'range) := Select_Subsystem(p,mix,mic);
911: sq : Laur_Sys(q'range);
912: pq : Poly_Sys(q'range);
913: qsols : Solution_List;
914: len : natural := 0;
915: fail : boolean;
916:
917: begin
918: Reduce(q'last+1,q); sq := Shift(q);
919: Solve(sq,qsols,fail);
920: if fail
921: then if mic.sub = null
922: then pq := Laurent_to_Polynomial_System(sq);
923: qsols := Solve_by_Static_Lifting(pq);
924: Clear(pq);
925: else Refined_Mixed_Solve(q,mix,mic,qsols);
926: end if;
927: Set_Continuation_Parameter(qsols,Create(0.0));
928: end if;
929: len := Length_Of(qsols);
930: if len > 0
931: then Mixed_Continuation(p,mic.nor.all,qsols);
932: Concat(sols,sols_last,qsols);
933: end if;
934: Clear(sq); Clear(q); Shallow_Clear(qsols);
935: end Mixed_Solve;
936:
937: procedure Mixed_Solve
938: ( file : in file_type; p : in Laur_Sys;
939: mix : in Vector; mic : in Mixed_Cell;
940: sols,sols_last : in out Solution_List ) is
941:
942: q : Laur_Sys(p'range) := Select_Subsystem(p,mix,mic);
943: sq : Laur_Sys(q'range);
944: pq : Poly_Sys(q'range);
945: qsols : Solution_List;
946: len : natural := 0;
947: fail : boolean;
948:
949: begin
950: Reduce(q'last+1,q); sq := Shift(q);
951: Solve(sq,qsols,fail);
952: if not fail
953: then put_line(file,"It is a fewnomial system.");
954: else put_line(file,"No fewnomial system.");
955: if mic.sub = null
956: then put_line(file,"Calling the black box solver.");
957: pq := Laurent_to_Polynomial_System(sq);
958: qsols := Solve_by_Static_Lifting(file,pq);
959: Clear(pq);
960: else put_line(file,"Using the refinement of the cell.");
961: Refined_Mixed_Solve(file,q,mix,mic,qsols);
962: end if;
963: Set_Continuation_Parameter(qsols,Create(0.0));
964: end if;
965: len := Length_Of(qsols);
966: put(file,len,1); put_line(file," solutions found.");
967: if len > 0
968: then Mixed_Continuation(file,p,mic.nor.all,qsols);
969: Concat(sols,sols_last,qsols);
970: end if;
971: Clear(sq); Clear(q); Shallow_Clear(qsols);
972: end Mixed_Solve;
973:
974: procedure Mixed_Solve
975: ( p : in Laur_Sys; lifted : in Array_of_Lists;
976: h : in Eval_Coeff_Laur_Sys;
977: c : in Standard_Complex_VecVecs.VecVec;
978: e : in Exponent_Vectors_Array;
979: j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
980: mix : in Vector; mic : in Mixed_Cell;
981: sols,sols_last : in out Solution_List ) is
982:
983: q : Laur_Sys(p'range) := Select_Subsystem(p,mix,mic);
984: sq : Laur_Sys(q'range);
985: pq : Poly_Sys(q'range);
986: qsols : Solution_List;
987: len : natural := 0;
988: fail : boolean;
989:
990: begin
991: Reduce(q'last+1,q); sq := Shift(q);
992: if mic.sub = null
993: then Solve(sq,qsols,fail);
994: else fail := true;
995: end if;
996: if fail
997: then if mic.sub = null
998: then pq := Laurent_to_Polynomial_System(sq);
999: qsols := Solve_by_Static_Lifting(pq);
1000: Clear(pq);
1001: else -- Refined_Mixed_Solve(q,mix,mic,qsols);
1002: Refined_Mixed_Solve(q,mix,mic,h,c,e,j,m,qsols);
1003: end if;
1004: Set_Continuation_Parameter(qsols,Create(0.0));
1005: end if;
1006: len := Length_Of(qsols);
1007: if len > 0
1008: then Mixed_Continuation(mix,lifted,h,c,e,j,m,mic.nor.all,qsols);
1009: Concat(sols,sols_last,qsols);
1010: end if;
1011: Clear(sq); Clear(q); Shallow_Clear(qsols);
1012: end Mixed_Solve;
1013:
1014: procedure Mixed_Solve
1015: ( file : in file_type;
1016: p : in Laur_Sys; lifted : in Array_of_Lists;
1017: h : in Eval_Coeff_Laur_Sys;
1018: c : in Standard_Complex_VecVecs.VecVec;
1019: e : in Exponent_Vectors_Array;
1020: j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
1021: mix : in Vector; mic : in Mixed_Cell;
1022: sols,sols_last : in out Solution_List ) is
1023:
1024: q : Laur_Sys(p'range) := Select_Subsystem(p,mix,mic);
1025: sq : Laur_Sys(q'range);
1026: pq : Poly_Sys(q'range);
1027: qsols : Solution_List;
1028: len : natural := 0;
1029: fail : boolean;
1030:
1031: begin
1032: Reduce(q'last+1,q); sq := Shift(q);
1033: if mic.sub = null
1034: then Solve(sq,qsols,fail);
1035: else fail := true;
1036: end if;
1037: if not fail
1038: then put_line(file,"It is a fewnomial system.");
1039: else put_line(file,"No fewnomial system.");
1040: if mic.sub = null
1041: then put_line(file,"Calling the black box solver.");
1042: pq := Laurent_to_Polynomial_System(sq);
1043: qsols := Solve_by_Static_Lifting(file,pq);
1044: Clear(pq);
1045: else put_line(file,"Using the refinement of the cell.");
1046: -- Refined_Mixed_Solve(file,q,mix,mic,qsols);
1047: Refined_Mixed_Solve(file,q,mix,mic,h,c,e,j,m,qsols);
1048: end if;
1049: Set_Continuation_Parameter(qsols,Create(0.0));
1050: end if;
1051: len := Length_Of(qsols);
1052: put(file,len,1); put_line(file," solutions found.");
1053: if len > 0
1054: then Mixed_Continuation(file,mix,lifted,h,c,e,j,m,mic.nor.all,qsols);
1055: Concat(sols,sols_last,qsols);
1056: end if;
1057: Clear(sq); Clear(q); Shallow_Clear(qsols);
1058: end Mixed_Solve;
1059:
1060: -- THIRD LAYER :
1061:
1062: procedure Mixed_Solve
1063: ( p : in Laur_Sys;
1064: mix : in Vector; mixsub : in Mixed_Subdivision;
1065: sols : in out Solution_List ) is
1066:
1067: tmp : Mixed_Subdivision := mixsub;
1068: mic : Mixed_Cell;
1069: sols_last : Solution_List;
1070:
1071: begin
1072: while not Is_Null(tmp) loop
1073: mic := Head_Of(tmp);
1074: Mixed_Solve(p,mix,mic,sols,sols_last);
1075: tmp := Tail_Of(tmp);
1076: end loop;
1077: end Mixed_Solve;
1078:
1079: procedure Mixed_Solve
1080: ( file : in file_type; p : in Laur_Sys;
1081: mix : in Vector; mixsub : in Mixed_Subdivision;
1082: sols : in out Solution_List ) is
1083:
1084: tmp : Mixed_Subdivision := mixsub;
1085: mic : Mixed_Cell;
1086: sols_last : Solution_List;
1087: cnt : natural := 0;
1088:
1089: begin
1090: while not Is_Null(tmp) loop
1091: mic := Head_Of(tmp);
1092: cnt := cnt + 1;
1093: new_line(file);
1094: put(file,"*** PROCESSING SUBSYSTEM "); put(file,cnt,1);
1095: put_line(file," ***");
1096: new_line(file);
1097: Mixed_Solve(file,p,mix,mic,sols,sols_last);
1098: tmp := Tail_Of(tmp);
1099: end loop;
1100: end Mixed_Solve;
1101:
1102: procedure Mixed_Solve
1103: ( p : in Laur_Sys; lifted : in Array_of_Lists;
1104: h : in Eval_Coeff_Laur_Sys;
1105: c : in Standard_Complex_VecVecs.VecVec;
1106: e : in Exponent_Vectors_Array;
1107: j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
1108: mix : in Vector; mixsub : in Mixed_Subdivision;
1109: sols : in out Solution_List ) is
1110:
1111: tmp : Mixed_Subdivision := mixsub;
1112: mic : Mixed_Cell;
1113: sols_last : Solution_List;
1114:
1115: begin
1116: while not Is_Null(tmp) loop
1117: mic := Head_Of(tmp);
1118: Mixed_Solve(p,lifted,h,c,e,j,m,mix,mic,sols,sols_last);
1119: tmp := Tail_Of(tmp);
1120: end loop;
1121: end Mixed_Solve;
1122:
1123: procedure Mixed_Solve
1124: ( file : in file_type;
1125: p : in Laur_Sys; lifted : in Array_of_Lists;
1126: h : in Eval_Coeff_Laur_Sys;
1127: c : in Standard_Complex_VecVecs.VecVec;
1128: e : in Exponent_Vectors_Array;
1129: j : in Eval_Coeff_Jaco_Mat; m : in Mult_Factors;
1130: mix : in Vector; mixsub : in Mixed_Subdivision;
1131: sols : in out Solution_List ) is
1132:
1133: tmp : Mixed_Subdivision := mixsub;
1134: mic : Mixed_Cell;
1135: sols_last : Solution_List;
1136: cnt : natural := 0;
1137:
1138: begin
1139: while not Is_Null(tmp) loop
1140: mic := Head_Of(tmp);
1141: cnt := cnt + 1;
1142: new_line(file);
1143: put(file,"*** PROCESSING SUBSYSTEM "); put(file,cnt,1);
1144: put_line(file," ***");
1145: new_line(file);
1146: Mixed_Solve(file,p,lifted,h,c,e,j,m,mix,mic,sols,sols_last);
1147: tmp := Tail_Of(tmp);
1148: end loop;
1149: end Mixed_Solve;
1150:
1151: end Integer_Polyhedral_Continuation;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to integer_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