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