Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Implift/mixed_homotopy_continuation.adb, Revision 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