Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Implift/mixed_homotopy_continuation.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
2: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
3: with Standard_Complex_Numbers_Polar; use Standard_Complex_Numbers_Polar;
4: with Standard_Random_Numbers; use Standard_Random_Numbers;
5: with Standard_Integer_Vectors;
6: with Standard_Integer_Vectors_io; use Standard_Integer_Vectors_io;
7: with Standard_Complex_Vectors;
8: with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
9: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
10: with Standard_Complex_Polynomials;
11: with Standard_Complex_Laur_Polys; use Standard_Complex_Laur_Polys;
12: with Standard_Complex_Poly_Systems; use Standard_Complex_Poly_Systems;
13: with Standard_Complex_Poly_Functions; use Standard_Complex_Poly_Functions;
14: with Standard_Complex_Poly_SysFun; use Standard_Complex_Poly_SysFun;
15: with Standard_Complex_Jaco_Matrices; use Standard_Complex_Jaco_Matrices;
16: with Standard_Poly_Laur_Convertors; use Standard_Poly_Laur_Convertors;
17: with Standard_Laur_Poly_Convertors; use Standard_Laur_Poly_Convertors;
18: with Lists_of_Integer_Vectors; use Lists_of_Integer_Vectors;
19: with Integer_Support_Functions; use Integer_Support_Functions;
20: with Arrays_of_Integer_Vector_Lists; use Arrays_of_Integer_Vector_Lists;
21: with Increment_and_Fix_Continuation; use Increment_and_Fix_Continuation;
22: with Standard_Root_Refiners; use Standard_Root_Refiners;
23: with Integer_Vectors_Utilities; use Integer_Vectors_Utilities;
24: with Lists_of_Vectors_Utilities; use Lists_of_Vectors_Utilities;
25: with Arrays_of_Lists_Utilities; use Arrays_of_Lists_Utilities;
26: with Transformations; use Transformations;
27: with Transforming_Solutions; use Transforming_Solutions;
28: with Transforming_Laurent_Systems; use Transforming_Laurent_Systems;
29: with Transforming_Integer_Vector_Lists; use Transforming_Integer_Vector_Lists;
30: with Power_Lists; use Power_Lists;
31: with Volumes;
32: with Durand_Kerner;
33:
34: package body Mixed_Homotopy_Continuation is
35:
36: -- INVARIANT CONDITION :
37: -- The procedures and functions in this package `mirror' corresponding
38: -- routines in the package Volumes.
39:
40: -- AUXILIARIES :
41:
42: type bar is array ( integer range <> ) of boolean;
43:
44: function Interchange2 ( p : Laur_Sys; index : integer ) return Laur_Sys is
45:
46: -- DESCRIPTION :
47: -- Returns a polynomial system where the first equation is interchanged
48: -- with the equation given by the index.
49:
50: res : Laur_Sys(p'range);
51:
52: begin
53: if index = p'first
54: then res := p;
55: else res(res'first) := p(index);
56: res(index) := p(p'first);
57: res(res'first+1..index-1) := p(p'first+1..index-1);
58: res(index+1..res'last) := p(index+1..p'last);
59: end if;
60: return res;
61: end Interchange2;
62:
63: procedure Interchange2 ( adl : in out Array_of_Lists; index : integer ) is
64:
65: -- DESCRIPTION :
66: -- Interchanges the first list with the list given by the index.
67:
68: tmp : List;
69:
70: begin
71: if index /= adl'first
72: then tmp := adl(adl'first);
73: adl(adl'first) := adl(index);
74: adl(index) := tmp;
75: end if;
76: end Interchange2;
77:
78: function Permute ( perm : Standard_Integer_Vectors.Vector; p : Laur_Sys )
79: return laur_Sys is
80:
81: -- DESCRIPTION :
82: -- Returns a permuted polynomial system.
83:
84: res : Laur_Sys(p'range);
85:
86: begin
87: for i in p'range loop
88: res(i) := p(perm(i));
89: end loop;
90: return res;
91: end Permute;
92:
93: function Initial_Degrees ( p : Poly ) return Degrees is
94:
95: -- DESCRIPTION :
96: -- Returns the initial degrees of the polynomial p.
97:
98: init : Degrees;
99:
100: procedure Init_Term ( t : in Term; cont : out boolean ) is
101: begin
102: init := new Standard_Integer_Vectors.Vector'(t.dg.all);
103: cont := false;
104: end Init_Term;
105: procedure Initial_Term is new Visiting_Iterator (Init_Term);
106:
107: begin
108: Initial_Term(p);
109: return init;
110: end Initial_Degrees;
111:
112: procedure Binomial ( p : in Poly;
113: d : out Standard_Integer_Vectors.Link_to_Vector;
114: k : in out integer; c : in out Complex_Number ) is
115:
116: -- DESCRIPTION :
117: -- p consists of two terms, this procedure gets the degrees d and
118: -- the constant c of the binomial equation.
119:
120: first : boolean := true;
121: dd : Degrees;
122:
123: procedure Scan_Term ( t : in Term; cont : out boolean ) is
124: begin
125: if first
126: then dd := new Standard_Integer_Vectors.Vector'(t.dg.all);
127: c := t.cf;
128: first := false;
129: else k := dd'first - 1;
130: for i in dd'range loop
131: dd(i) := dd(i) - t.dg(i);
132: if k < dd'first and then dd(i) /= 0
133: then k := i;
134: end if;
135: end loop;
136: c := -t.cf/c;
137: end if;
138: cont := true;
139: end Scan_Term;
140: procedure Scan_Terms is new Visiting_Iterator (Scan_Term);
141:
142: begin
143: Scan_Terms(p);
144: d := new Standard_Integer_Vectors.Vector'(dd.all);
145: Standard_Integer_Vectors.Clear(Standard_Integer_Vectors.Link_to_Vector(dd));
146: end Binomial;
147:
148: procedure Normalize ( p : in Laur_Sys; dl : in out List;
149: wp : in out Laur_Sys; shifted : out boolean ) is
150:
151: -- DESCRIPTION :
152: -- Makes sure that the first element of dl contains all zeroes.
153:
154: -- REQUIRED :
155: -- The list dl is not empty.
156:
157: -- ON ENTRY :
158: -- p a Laurent polynomial system;
159: -- dl the power list of p(p'first).
160:
161: -- ON RETURN :
162: -- dl the first element of contains all zeroes,
163: -- therefore dl has been shifted;
164: -- wp a Laurent polynomial system where
165: -- dl is the power list of wp(wp'first);
166: -- shifted is true if dl has been shifted.
167:
168: use Standard_Integer_Vectors;
169:
170: first : Link_to_Vector := Head_Of(dl);
171: nullvec : Vector(first'range) := (first'range => 0);
172: shiftvec : Link_to_Vector;
173:
174: begin
175: if not Is_In(dl,nullvec)
176: then shiftvec := Graded_Max(dl);
177: Shift(dl,shiftvec);
178: Clear(shiftvec);
179: Copy(p,wp);
180: for i in p'range loop
181: Shift(wp(i));
182: end loop;
183: else wp := p;
184: shifted := false;
185: end if;
186: Move_to_Front(dl,nullvec);
187: end Normalize;
188:
189: function Evaluate ( p : Poly; x : Complex_Number; k : integer )
190: return Poly is
191:
192: -- DESCRIPTION :
193: -- This function returns a polynomial where the kth unknown has
194: -- been replaced by x.
195: -- It is important to use this function as the term order in p
196: -- must remain the same!
197:
198: res : Poly;
199:
200: procedure Evaluate_Term ( t : in out Term; cont : out boolean ) is
201:
202: fac : Complex_Number;
203: pow : natural;
204:
205: begin
206: if t.dg(k) < 0
207: then fac := Create(1.0)/x;
208: pow := -t.dg(k);
209: else fac := x;
210: pow := t.dg(k);
211: end if;
212: for i in 1..pow loop
213: t.cf := t.cf * fac;
214: end loop;
215: declare
216: tmp : constant Standard_Integer_Vectors.Vector := t.dg.all;
217: begin
218: Clear(t);
219: t.dg := new Standard_Integer_Vectors.Vector'(Reduce(tmp,k));
220: end;
221: cont := true;
222: end Evaluate_Term;
223: procedure Evaluate_Terms is new Changing_Iterator (Evaluate_Term);
224:
225: begin
226: Copy(p,res);
227: Evaluate_Terms(res);
228: return res;
229: end Evaluate;
230:
231: function Re_Arrange ( p : poly ) return Poly is
232:
233: -- DESCRIPTION :
234: -- Returns a polynomial whose terms are sorted
235: -- in graded lexicographical ordening.
236:
237: res : Poly := Null_Poly;
238:
239: procedure Re_Arrange_Term ( t : in Term; cont : out boolean ) is
240: begin
241: Add(res,t);
242: cont := true;
243: end Re_Arrange_Term;
244: procedure Re_Arrange_Terms is new Visiting_Iterator (Re_Arrange_Term);
245:
246: begin
247: Re_Arrange_Terms(p);
248: return res;
249: end Re_Arrange;
250:
251: function Substitute ( p : Poly; v : Standard_Complex_Vectors.Vector;
252: k : integer ) return Poly is
253:
254: -- DESCRIPTION :
255: -- Substitutes the values in v into the polynomial p,
256: -- starting at the last unknowns of p.
257:
258: init : Degrees := Initial_Degrees(p);
259: index : integer := init'last;
260: res,tmp : Poly;
261:
262: begin
263: if index = k
264: then index := index - 1;
265: res := Evaluate(p,v(v'last),index);
266: else res := Evaluate(p,v(v'last),index);
267: end if;
268: for i in reverse v'first..(v'last-1) loop
269: index := index - 1;
270: if index = k
271: then index := index - 1;
272: end if;
273: tmp := Evaluate(res,v(i),index);
274: Clear(res); Copy(tmp,res); Clear(tmp);
275: end loop;
276: Standard_Integer_Vectors.Clear
277: (Standard_Integer_Vectors.Link_to_Vector(init));
278: return res;
279: end Substitute;
280:
281: procedure Refine_and_Concat
282: ( p : in Laur_Sys;
283: newsols,sols,sols_last : in out Solution_List ) is
284:
285: -- DESCRIPTION :
286: -- This procedure refines the solutions of a given
287: -- polynomial system and adds them to the solution list.
288: -- This can be very useful for eliminating rounding errors
289: -- after transformating the solutions.
290:
291: pp : Poly_Sys(p'range) := Laurent_to_Polynomial_System(p);
292: numb : natural := 0;
293:
294: begin
295: Silent_Root_Refiner(pp,newsols,10.0**(-12),10.0**(-12),10.0**(-8),numb,5);
296: Concat(sols,sols_last,newsols);
297: Clear(pp); Shallow_Clear(newsols);
298: end Refine_and_Concat;
299:
300: procedure Refine_and_Concat
301: ( file : in file_type; p : in Laur_Sys;
302: newsols,sols,sols_last : in out Solution_List ) is
303:
304: -- DESCRIPTION :
305: -- This procedure refines the solutions of a given
306: -- polynomial system and adds them to the solution list.
307: -- This can be very useful for eliminating rounding errors
308: -- after transformating the solutions.
309:
310: pp : Poly_Sys(p'range) := Laurent_to_Polynomial_System(p);
311: numb : natural := 0;
312:
313: begin
314: Reporting_Root_Refiner
315: (file,pp,newsols,10.0**(-12),10.0**(-12),10.0**(-8),numb,5,false);
316: Concat(sols,sols_last,newsols);
317: Clear(pp); Shallow_Clear(newsols);
318: end Refine_and_Concat;
319:
320: procedure Write_Direction
321: ( file : in file_type;
322: v : in Standard_Integer_Vectors.Link_to_Vector ) is
323: begin
324: put(file,"++++ considering direction "); put(file,v);
325: put_line(file," ++++");
326: end Write_Direction;
327:
328: -- INTERMEDIATE LAYER :
329:
330: procedure Mixed_Continue
331: ( file : in file_type; p : in Laur_Sys;
332: k : in integer; m : in Standard_Integer_Vectors.Vector;
333: sols : in out Solution_List ) is
334:
335: -- DESCRIPTION :
336: -- This continuation routine computes a part of the solution list
337: -- of a Laurent polynomial system.
338:
339: -- ON ENTRY :
340: -- file a file where the intermediate results are written;
341: -- p the transformed Laurent polynomial system to be solved;
342: -- k the index;
343: -- m m(k) = p1(v), m(k/=l) = Maximal_Degree(p(l),k);
344: -- sols the start solutions.
345:
346: -- ON RETURN :
347: -- sols the computed solutions.
348:
349: h : Laur_Sys(p'range);
350: hp : Poly_Sys(h'range);
351: hpe : Eval_Poly_Sys(hp'range);
352: j : Jaco_Mat(p'range,p'first..p'last+1);
353: je : Eval_Jaco_Mat(j'range(1),j'range(2));
354:
355: function Construct_Homotopy
356: ( p : Laur_Sys; k : integer;
357: m : Standard_Integer_Vectors.Vector ) return Laur_Sys is
358:
359: res : Laur_Sys(p'range);
360: ran : Complex_Number;
361: re : boolean;
362: zeroes : Degrees := new Standard_Integer_Vectors.Vector'(p'range => 0);
363:
364: function Construct_First_Polynomial
365: ( pp : Poly; max : integer ) return Poly is
366:
367: r : Poly := Null_Poly;
368:
369: procedure Construct_Term ( t : in Term; cont : out boolean ) is
370:
371: rt : Term;
372:
373: begin
374: rt.cf := t.cf;
375: rt.dg := new Standard_Integer_Vectors.Vector(t.dg'first..t.dg'last+1);
376: rt.dg(t.dg'range) := t.dg.all;
377: rt.dg(k) := -t.dg(k) + max;
378: if Equal(t.dg,zeroes)
379: then rt.dg(rt.dg'last) := 0;
380: re := ( IMAG_PART(rt.cf) + 1.0 = 1.0 );
381: else rt.dg(rt.dg'last) := rt.dg(k);
382: end if;
383: Add(r,rt);
384: Clear(rt);
385: cont := true;
386: end Construct_Term;
387: procedure Construct_Terms is new Visiting_Iterator(Construct_Term);
388:
389: begin
390: Construct_Terms(pp);
391: Standard_Integer_Vectors.Clear
392: (Standard_Integer_Vectors.Link_to_Vector(zeroes));
393: return r;
394: end Construct_First_Polynomial;
395:
396: function Construct_Polynomial ( pp : Poly; max : integer ) return Poly is
397:
398: r : Poly := Null_Poly;
399:
400: procedure Construct_Term ( t : in Term; cont : out boolean ) is
401:
402: rt : Term;
403:
404: begin
405: rt.cf := t.cf;
406: rt.dg := new Standard_Integer_Vectors.Vector(t.dg'first..t.dg'last+1);
407: rt.dg(t.dg'range) := t.dg.all;
408: rt.dg(k) := -t.dg(k) + max;
409: rt.dg(rt.dg'last) := rt.dg(k);
410: Add(r,rt);
411: Clear(rt);
412: cont := true;
413: end Construct_Term;
414: procedure Construct_Terms is new Visiting_Iterator(Construct_Term);
415:
416: begin
417: Construct_Terms(pp);
418: return r;
419: end Construct_Polynomial;
420:
421: begin
422: res(res'first) := Construct_First_Polynomial(p(p'first),m(m'first));
423: for i in p'first+1..p'last loop
424: res(i) := Construct_Polynomial(p(i),m(i));
425: end loop;
426: if re
427: then for i in res'range loop
428: ran := Random;
429: Mul(res(i),ran);
430: end loop;
431: end if;
432: return res;
433: end Construct_Homotopy;
434:
435: function To_Be_Removed ( flag : in natural ) return boolean is
436: begin
437: return ( flag /= 1 );
438: end To_Be_Removed;
439: procedure Extract_Regular_Solutions is
440: new Standard_Complex_Solutions.Delete(To_Be_Removed);
441:
442: begin
443: h := Construct_Homotopy(p,k,m); -- CONSTRUCTION OF HOMOTOPY
444: hp := Laurent_to_Polynomial_System(h);
445: hpe := Create(hp);
446: j := Create(hp);
447: je := Create(j);
448: declare -- CONTINUATION
449:
450: use Standard_Complex_Vectors;
451:
452: function Eval ( x : Vector; t : Complex_Number ) return Vector is
453:
454: xt : Vector(x'first..x'last+1);
455:
456: begin
457: xt(x'range) := x;
458: xt(xt'last) := t;
459: return Eval(hpe,xt);
460: end Eval;
461:
462: function dHt ( x : Vector; t : Complex_Number ) return Vector is
463:
464: xt : Vector(x'first..x'last+1);
465: res : Vector(p'range);
466:
467: begin
468: xt(x'range) := x;
469: xt(xt'last) := t;
470: for i in res'range loop
471: res(i) := Eval(je(i,xt'last),xt);
472: end loop;
473: return res;
474: end dHt;
475:
476: function dHx ( x : Vector; t : Complex_Number ) return Matrix is
477:
478: xt : Vector(x'first..x'last+1);
479: m : Matrix(x'range,x'range);
480:
481: begin
482: xt(x'range) := x;
483: xt(xt'last) := t;
484: for i in m'range(1) loop
485: for j in m'range(1) loop
486: m(i,j) := Eval(je(i,j),xt);
487: end loop;
488: end loop;
489: return m;
490: end dHx;
491:
492: procedure Invert ( k : in integer; sols : in out Solution_List ) is
493:
494: -- DESCRIPTION :
495: -- For all solutions s in sols : s.v(k) := 1/s.v(k).
496:
497: tmp : Solution_List := sols;
498:
499: begin
500: while not Is_Null(tmp) loop
501: declare
502: l : Link_to_Solution := Head_Of(tmp);
503: begin
504: l.v(k) := Create(1.0)/l.v(k);
505: l.t := Create(0.0);
506: end;
507: tmp := Tail_Of(tmp);
508: end loop;
509: end Invert;
510:
511: procedure Cont is new Reporting_Continue(Max_Norm,Eval,dHt,dHx);
512:
513: begin
514: Invert(k,sols);
515: Cont(file,sols,false);
516: Invert(k,sols);
517: Extract_Regular_Solutions(sols);
518: end;
519: Clear(h); Clear(hp); Clear(hpe); Clear(j); Clear(je);
520: end Mixed_Continue;
521:
522: -- THE SOLVERS :
523:
524: function One_Unknown_Solve ( p : Poly )
525: return Standard_Complex_Vectors.Vector is
526:
527: -- DESCRIPTION :
528: -- Returns the solution vector of p, a polynomial in one unknown.
529: -- p will be solved by using the method of Durand-Kerner.
530:
531: p1 : Poly := Re_Arrange(p);
532: init : Degrees := Initial_Degrees(p1);
533: index : integer := init'first;
534: min : integer := -Minimal_Degree(p1,index);
535: pv : Standard_Complex_Vectors.Vector(0..Degree(p1)+min);
536: z,res : Standard_Complex_Vectors.Vector(1..pv'last);
537: maxsteps : constant natural := 10;
538: eps : constant double_float := 10.0**(-10);
539: nb : integer := pv'last + 1;
540:
541: procedure Store_Coeff ( t : in Term; cont : out boolean ) is
542: begin
543: nb := nb - 1;
544: if t.dg(index) = (nb - min)
545: then pv(nb) := t.cf;
546: else for i in reverse nb..(t.dg(index)+min+1) loop
547: pv(i) := Create(0.0);
548: end loop;
549: nb := t.dg(index) + min;
550: pv(nb) := t.cf;
551: end if;
552: cont := true;
553: end Store_Coeff;
554: procedure Polynomial_To_Vector is new Visiting_Iterator (Store_Coeff);
555:
556: procedure Write ( step : in natural;
557: z,res : in Standard_Complex_Vectors.Vector ) is
558: begin
559: null; -- no output desired during the iterations
560: end Write;
561: procedure DuKe is new Durand_Kerner (Write);
562:
563: begin
564: for i in pv'range loop -- initialize coefficient vector
565: pv(i) := Create(0.0);
566: end loop;
567: Standard_Integer_Vectors.Clear
568: (Standard_Integer_Vectors.Link_to_Vector(init));
569: Polynomial_To_Vector(p1);
570: Clear(p1);
571: for i in z'range loop
572: z(i) := Random;
573: res(i) := z(i);
574: end loop;
575: DuKe(pv,z,res,maxsteps,eps,nb);
576: return z;
577: end One_Unknown_Solve;
578:
579: procedure One_Unknown_Solve ( p : in Poly; sols : in out Solution_List ) is
580:
581: -- DESCRIPTION :
582: -- If p is a polynomial in one unknown,
583: -- p can be solved efficiently by the application of Durand-Kerner.
584:
585: init : Degrees := Initial_Degrees(p);
586:
587: function Make_Solutions ( x : in Standard_Complex_Vectors.Vector )
588: return Solution_List is
589:
590: res,res_last : Solution_List;
591: s : Solution(1);
592:
593: begin
594: s.m := 1;
595: s.t := Create(0.0);
596: for i in x'range loop
597: s.v(init'first) := x(i);
598: Append(res,res_last,s);
599: end loop;
600: return res;
601: end Make_Solutions;
602:
603: begin
604: sols := Make_Solutions(One_Unknown_Solve(p));
605: Standard_Integer_Vectors.Clear
606: (Standard_Integer_Vectors.Link_to_Vector(init));
607: end One_Unknown_Solve;
608:
609: procedure Two_Terms_Solve
610: ( file : in file_type; p : in Laur_Sys;
611: tv : in Tree_of_Vectors; bkk : out natural;
612: sols : in out Solution_List ) is
613:
614: -- DESCRIPTION :
615: -- The first polynomial of p consists of two terms.
616: -- A binomial system can be solved efficiently by
617: -- transforming and using de Moivre's rule.
618:
619: d : Standard_Integer_Vectors.Link_to_Vector;
620: c : Complex_Number := Create(0.0);
621: k : natural := 0;
622: sols_last : Solution_List;
623:
624: begin
625: -- put_line(file,"Applying Two_Terms_Solve on "); Write(file,p);
626: Binomial(p(p'first),d,k,c);
627: if k < d'first
628: then bkk := 0;
629: Standard_Integer_Vectors.Clear(d); return;
630: elsif ( c + Create(1.0) = Create(1.0) )
631: then bkk := 0;
632: Standard_Integer_Vectors.Clear(d); return;
633: elsif d(k) < 0
634: then Standard_Integer_Vectors.Min(d);
635: c := Create(1.0)/c;
636: end if;
637: declare
638: t : Transfo := Rotate(d,k);
639: tmp_bkk : natural := 0;
640: begin
641: Apply(t,d);
642: declare
643: v : Standard_Complex_Vectors.Vector(1..d(k));
644: tmp : Poly;
645: rtp : Laur_Sys(p'first..(p'last-1));
646: rtp_bkk : natural;
647: rtp_sols : Solution_List;
648: begin
649: for i in v'range loop
650: v(i) := Root(c,d(k),i);
651: for j in rtp'range loop
652: tmp := Transform(t,p(j+1));
653: rtp(j) := Evaluate(tmp,v(i),k);
654: Clear(tmp);
655: end loop;
656: Solve(file,rtp,tv,rtp_bkk,rtp_sols);
657: Clear(rtp);
658: tmp_bkk := tmp_bkk + rtp_bkk;
659: Insert(v(i),k,rtp_sols);
660: Transform(t,rtp_sols);
661: --Concat(sols,sols_last,rtp_sols);
662: Refine_and_Concat(file,p,rtp_sols,sols,sols_last);
663: end loop;
664: end;
665: Clear(t);
666: bkk := tmp_bkk;
667: end;
668: -- end if;
669: Standard_Integer_Vectors.Clear(d);
670: -- put_line(file,"The solutions found : "); put(file,sols);
671: end Two_Terms_Solve;
672:
673: function Project_on_First_and_Solve
674: ( p : Poly; k : integer; sols : Solution_List )
675: return Solution_List is
676:
677: -- ON ENTRY :
678: -- p a Laurent polynomial in n unknowns x1,..,xk,..,xn;
679: -- sols contains values for x1,..,xn, except xk.
680:
681: -- ON RETURN :
682: -- a solution list for p, obtained after substition of the values
683: -- for x1,..,xn into the polynomial p.
684:
685: tmp,res,res_last : Solution_List;
686: init : Degrees := Initial_Degrees(p);
687:
688: begin
689: -- put_line(file,"Calling Project_on_First_and_Solve");
690: -- put_line(file," with polynomial : ");
691: -- put(file,Laurent_Polynomial_to_Polynomial(p)); new_line(file);
692: tmp := sols;
693: while not Is_Null(tmp) loop
694: declare
695: p1 : Poly := Substitute(p,Head_Of(tmp).v,k);
696: sols1 : Solution_List;
697: begin
698: -- put(file,"k : "); put(file,k,1); new_line(file);
699: -- put(file,"v : "); put(file,Head_Of(tmp).v,3,3,3); new_line(file);
700: -- put_line(file,"After substitution : "); Write(file,p1);
701: if Number_of_Terms(p1) < 2
702: then null;
703: elsif Number_of_Terms(p1) = 2
704: then
705: declare
706: d : Standard_Integer_Vectors.Link_to_Vector;
707: l : natural := 0;
708: c : Complex_Number := Create(0.0);
709: begin
710: Binomial(p1,d,l,c);
711: if l < init'first
712: then null;
713: elsif ( c + Create(1.0) = Create(1.0) )
714: then null;
715: else
716: if d(l) < 0
717: then d(l) := -d(l); c := Create(1.0)/c;
718: end if;
719: declare
720: v : Standard_Complex_Vectors.Vector(1..d(l));
721: begin
722: for i in v'range loop
723: v(i) := Root(c,d(l),i);
724: end loop;
725: sols1 := Insert(v,k,Head_Of(tmp).all);
726: -- put_line(file,"Sols1 after Binomial :");
727: -- put(file,sols1);
728: Concat(res,res_last,sols1);
729: Shallow_Clear(sols1);
730: end;
731: end if;
732: Standard_Integer_Vectors.Clear(d);
733: end;
734: else
735: sols1 := Insert(One_Unknown_Solve(p1),k,Head_Of(tmp).all);
736: Concat(res,res_last,sols1);
737: -- put_line(file,"Sols1 :"); put(file,sols1);
738: Shallow_Clear(sols1);
739: end if;
740: Clear(p1);
741: end;
742: tmp := Tail_Of(tmp);
743: end loop;
744: Standard_Integer_Vectors.Clear
745: (Standard_Integer_Vectors.Link_to_Vector(init));
746: return res;
747: end Project_on_First_and_Solve;
748:
749: procedure Project_and_Solve
750: ( file : in file_type; p : in Laur_Sys; k : in integer;
751: m : in out Standard_Integer_Vectors.Vector;
752: nd : in node; bkk : out natural;
753: sols : in out Solution_List ) is
754:
755: -- DESCRIPTION :
756: -- Solves the projected start system of p along a direction v.
757:
758: -- ON ENTRY :
759: -- file a file to write intermediate results on;
760: -- p a Laurent polynomial system;
761: -- k entry in the degrees of p;
762: -- m m(m'first) equals Maximal_Support(p(p'first),v) > 0;
763: -- nd a node in the tree of useful directions.
764:
765: -- ON RETURN :
766: -- m m(m'first+1..m'last) contains the maximal degrees
767: -- of the last n-1 equations of p in xk;
768: -- bkk the BKK bound of the projected system;
769: -- sols the solutions of the projected system.
770:
771: g_v : Laur_Sys(p'first..(p'last-1));
772: bkk_g_v : natural;
773: sols_g_v : Solution_List;
774:
775: begin
776: -- put_line(file,"Applying Project_and_Solve on"); Write(file,p);
777: for i in g_v'range loop
778: m(i+1) := Maximal_Degree(p(i+1),k);
779: g_v(i) := Face(k,m(i+1),p(i+1));
780: Reduce(k,g_v(i));
781: end loop;
782: if (nd.ltv = null) or else Is_Null(nd.ltv.all)
783: then Solve(file,g_v,bkk_g_v,sols_g_v);
784: else Solve(file,g_v,nd.ltv.all,bkk_g_v,sols_g_v);
785: end if;
786: -- put(file,"After Solve (without tv) bkk_g_v = "); put(file,bkk_g_v,1);
787: -- new_line(file);
788: declare
789: p0 : Poly := Re_Arrange(p(p'first));
790: p1 : Poly := Face(k,m(m'first),p0);
791: cnst : Term;
792: begin
793: cnst.dg := new Standard_Integer_Vectors.Vector'(p'range => 0);
794: if Coeff(p1,cnst.dg) = Create(0.0)
795: then cnst.cf := Coeff(p0,cnst.dg);
796: Add(p1,cnst);
797: end if;
798: Standard_Integer_Vectors.Clear
799: (Standard_Integer_Vectors.Link_to_Vector(cnst.dg));
800: Clear(p0);
801: sols := Project_on_First_and_Solve(p1,k,sols_g_v);
802: -- sols := Project_on_First_and_Solve(file,p1,k,sols_g_v);
803: Set_Continuation_Parameter(sols,Create(0.0));
804: Clear(p1);
805: end;
806: bkk := m(m'first)*bkk_g_v;
807: Clear(sols_g_v);
808: Clear(g_v);
809: -- put_line(file,"The solutions found : "); put(file,sols);
810: end Project_and_Solve;
811:
812: procedure Unmixed_Solve
813: ( file : in file_type; p : in Laur_Sys; dl : in List;
814: tv : in Tree_of_Vectors;
815: bkk : out natural; sols : in out Solution_List ) is
816:
817: -- DESCRIPTION :
818: -- Solves a Laurent polynomial system where all polytopes are the same.
819:
820: -- ON ENTRY :
821: -- file a file to write intermediate results on;
822: -- p a Laurent polynomial system;
823: -- dl the list of powers for p;
824: -- tv the tree of degrees containing the useful directions.
825:
826: -- ON RETURN :
827: -- bkk the bkk bound;
828: -- sols the list of solutions.
829:
830: sols_last : Solution_List;
831: tmp_bkk : natural := 0;
832: tmp : Tree_of_Vectors := tv;
833:
834: begin
835: tmp_bkk := 0;
836: tmp := tv;
837: while not Is_Null(tmp) loop
838: declare
839: nd : node := Head_Of(tmp);
840: v : Standard_Integer_Vectors.Link_to_Vector := nd.d;
841: i : integer := Pivot(v);
842: pv : integer := Maximal_Support(dl,v.all);
843: t : Transfo := Build_Transfo(v,i);
844: tp : Laur_Sys(p'range) := Transform(t,p);
845: bkk_tp : natural;
846: sols_tp : Solution_List;
847: max : Standard_Integer_Vectors.Vector(p'range);
848: begin
849: Write_Direction(file,v);
850: max(max'first) := pv;
851: -- if (nd.ltv = null) or else Is_Null(nd.ltv.all)
852: -- then Projected_Solve(file,tp,i,max,bkk_tp,sols_tp);
853: -- else Projected_Solve
854: -- (file,tp,i,max,nd.ltv.all,bkk_tp,sols_tp);
855: -- end if;
856: Project_and_Solve(file,tp,i,max,nd,bkk_tp,sols_tp);
857: Mixed_Continue(file,tp,i,max,sols_tp);
858: tmp_bkk := tmp_bkk + bkk_tp;
859: Transform(t,sols_tp);
860: --Concat(sols,sols_last,sols_tp);
861: Refine_and_Concat(file,p,sols_tp,sols,sols_last);
862: Clear(t); Clear(tp);
863: end;
864: tmp := Tail_Of(tmp);
865: end loop;
866: bkk := tmp_bkk;
867: end Unmixed_Solve;
868:
869: procedure Unmixed_Solve
870: ( file : in file_type; p : in Laur_Sys; dl : in List;
871: bkk : out natural; sols : in out Solution_List ) is
872:
873: tv : Tree_of_Vectors;
874:
875: begin
876: Volumes.Volume(p'last,dl,tv,bkk);
877: Unmixed_Solve(file,p,dl,tv,bkk,sols);
878: Clear(tv);
879: end Unmixed_Solve;
880:
881: procedure Mixed_Solve
882: ( file : in file_type; p : in Laur_Sys;
883: adl : in out Array_of_Lists; tv : in Tree_of_Vectors;
884: bkk : out natural; sols : in out Solution_List ) is
885:
886: -- DESCRIPTION :
887: -- Computes the solutions of the Laurent polynomial system p,
888: -- where p has more than one equation.
889:
890: -- NOTE :
891: -- This procedure mirrors the procedure Minkowski_Sum in the body
892: -- of the package Volumes.
893:
894: tmp_bkk,len : natural;
895: tmp : Tree_of_Vectors;
896: index : integer := Index2(adl);
897: wp : Laur_Sys(p'range);
898: sols_last : Solution_List;
899: shifted : boolean;
900: perm,mix : Standard_Integer_Vectors.Link_to_Vector;
901:
902: begin
903: Interchange2(adl,index);
904: len := Length_Of(adl(adl'first));
905: -- put_line(file,"Applying Mixed_Solve on"); Write(file,p);
906: if len = 2
907: then wp := Interchange2(p,index);
908: Two_Terms_Solve(file,wp,tv,bkk,sols);
909: else
910: if len > 2
911: then -- INITIALIZATION :
912: Mixture(adl,perm,mix);
913: wp := Permute(perm.all,p);
914: declare
915: zeroes : Degrees
916: := new Standard_Integer_Vectors.Vector'(p'range => 0);
917: tmpwpi : Poly;
918: begin
919: if Coeff(wp(wp'first),zeroes) = Create(0.0)
920: then shifted := true;
921: -- wp(wp'first) := Shift(p(p'first));
922: Copy(p(index),tmpwpi); wp(wp'first) := tmpwpi;
923: Shift(wp(wp'first));
924: else shifted := false;
925: end if;
926: Standard_Integer_Vectors.Clear
927: (Standard_Integer_Vectors.Link_to_Vector(zeroes));
928: end;
929: -- MIXED HOMOTOPY CONTINUATION :
930: tmp_bkk := 0;
931: tmp := tv;
932: while not Is_Null(tmp) loop
933: declare
934: nd : node := Head_Of(tmp);
935: v : Standard_Integer_Vectors.Link_to_Vector := nd.d;
936: k : integer := Pivot(v);
937: pv : integer := Maximal_Support(wp(wp'first),v);
938: t : Transfo := Build_Transfo(v,k);
939: twp : Laur_Sys(wp'range) := Transform(t,wp);
940: bkk_twp : natural;
941: sols_twp : Solution_List;
942: m : Standard_Integer_Vectors.Vector(wp'range);
943: begin
944: Write_Direction(file,v);
945: m(m'first) := pv;
946: -- if (nd.ltv = null) or else Is_Null(nd.ltv.all)
947: -- then Projected_Solve(file,twp,k,m,bkk_twp,sols_twp);
948: -- else Projected_Solve
949: -- (file,twp,k,m,nd.ltv.all,bkk_twp,sols_twp);
950: -- end if;
951: Project_and_Solve(file,twp,k,m,nd,bkk_twp,sols_twp);
952: Mixed_Continue(file,twp,k,m,sols_twp);
953: tmp_bkk := tmp_bkk + bkk_twp;
954: Transform(t,sols_twp);
955: --Concat(sols,sols_last,sols_twp);
956: Refine_and_Concat(file,wp,sols_twp,sols,sols_last);
957: Clear(t); Clear(twp);
958: end;
959: tmp := Tail_Of(tmp);
960: end loop;
961: bkk := tmp_bkk;
962: Standard_Integer_Vectors.Clear(perm);
963: Standard_Integer_Vectors.Clear(mix);
964: if shifted
965: then Clear(wp(wp'first));
966: end if;
967: else -- len < 2
968: bkk := 0;
969: end if;
970: end if;
971: end Mixed_Solve;
972:
973: procedure Mixed_Solve
974: ( file : in file_type; p : in Laur_Sys;
975: adl : in out Array_of_Lists; bkk : out natural;
976: sols : in out Solution_List ) is
977:
978: tv : Tree_of_Vectors;
979:
980: begin
981: Volumes.Mixed_Volume(adl'last,adl,tv,bkk);
982: Mixed_Solve(file,p,adl,tv,bkk,sols);
983: Clear(tv);
984: end Mixed_Solve;
985:
986: procedure Solve ( file : in file_type; p : in Laur_Sys;
987: bkk : out natural; sols : in out Solution_List ) is
988:
989: al : Array_of_Lists(p'range) := Create(p);
990: tv : Tree_of_Vectors;
991:
992: begin
993: Volumes.Mixed_Volume(p'last,al,tv,bkk);
994: Solve(file,p,tv,bkk,sols);
995: Deep_Clear(al); Clear(tv);
996: end Solve;
997:
998: procedure Solve ( file : in file_type; p : in Laur_Sys;
999: tv : in Tree_of_Vectors; bkk : out natural;
1000: sols : in out Solution_List ) is
1001:
1002: -- NOTE :
1003: -- This procedure mirrors the procedure Volumes.Mixed_Volume,
1004: -- with a tree of useful directions on entry.
1005:
1006: begin
1007: if p'first = p'last
1008: then One_Unknown_Solve(p(p'first),sols);
1009: bkk := Length_Of(sols);
1010: else --if Is_Fewnomial_System(p)
1011: -- then
1012: -- declare
1013: -- fail : boolean;
1014: -- begin
1015: -- Fewnomials.Solve(p,sols,fail);
1016: -- if fail
1017: -- then bkk := 0; Clear(sols);
1018: -- else bkk := Length_Of(sols);
1019: -- end if;
1020: -- end;
1021: -- else
1022: declare
1023: adl : Array_of_Lists(p'range) := Create(p);
1024: begin
1025: if All_Equal(adl)
1026: then
1027: for i in (adl'first+1)..adl'last loop
1028: Deep_Clear(adl(i));
1029: end loop;
1030: declare
1031: wp : Laur_Sys(p'range);
1032: shifted : boolean;
1033: begin
1034: Normalize(p,adl(adl'first),wp,shifted);
1035: if Is_Null(tv)
1036: then Unmixed_Solve(file,wp,adl(adl'first),bkk,sols);
1037: else Unmixed_Solve(file,wp,adl(adl'first),tv,bkk,sols);
1038: end if;
1039: if shifted
1040: then Clear(wp);
1041: end if;
1042: end;
1043: elsif Is_Null(tv)
1044: then Mixed_Solve(file,p,adl,bkk,sols);
1045: else Mixed_Solve(file,p,adl,tv,bkk,sols);
1046: end if;
1047: end;
1048: end if;
1049: end Solve;
1050:
1051: end Mixed_Homotopy_Continuation;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to mixed_homotopy_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 / Implift