Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Product/interpolating_homotopies.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_Natural_Vectors;
! 4: with Standard_Complex_Vectors;
! 5: with Standard_Complex_Matrices;
! 6: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
! 7: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
! 8: with Standard_Complex_Poly_Randomizers; use Standard_Complex_Poly_Randomizers;
! 9: with Sets_of_Unknowns; use Sets_of_Unknowns;
! 10: with Degrees_in_Sets_of_Unknowns; use Degrees_in_Sets_of_Unknowns;
! 11:
! 12: package body Interpolating_Homotopies is
! 13:
! 14: -- AUXILIARY OPERATIONS :
! 15:
! 16: function Initial_Term ( p : Poly ) return Term is
! 17:
! 18: -- DESCRIPTION :
! 19: -- Returns the initial term of p.
! 20:
! 21: res : Term;
! 22:
! 23: procedure Init_Term ( t : Term; continue : out boolean ) is
! 24: begin
! 25: res := t;
! 26: continue := false;
! 27: end Init_Term;
! 28: procedure Init_Terms is new Visiting_Iterator(Init_Term);
! 29:
! 30: begin
! 31: Init_Terms(p);
! 32: return res;
! 33: end Initial_Term;
! 34:
! 35: procedure Eliminate_Term ( p : in out Poly; dg : in Degrees ) is
! 36:
! 37: -- DESCRIPTION :
! 38: -- Eliminates the monomial in p that has the given exponent vector.
! 39:
! 40: c : Complex_Number := Coeff(p,dg);
! 41:
! 42: begin
! 43: if c /= Create(0.0)
! 44: then declare
! 45: t : Term;
! 46: begin
! 47: t.cf := c; t.dg := dg;
! 48: Sub(p,t);
! 49: end;
! 50: end if;
! 51: end Eliminate_Term;
! 52:
! 53: function Admitted ( t : Term; i : natural; z : partition;
! 54: d : Standard_Integer_Matrices.Matrix ) return boolean is
! 55:
! 56: -- DESCRIPTION :
! 57: -- Returns true if the given term does not violate the m-homogeneous
! 58: -- structure, i.e., if Degree(t,z(j)) <= d(i,j), for j in z'range.
! 59:
! 60: begin
! 61: for j in z'range loop
! 62: if Degree(t,z(j)) > d(i,j)
! 63: then return false;
! 64: end if;
! 65: end loop;
! 66: return true;
! 67: end Admitted;
! 68:
! 69: function Dense_Representation ( p : Poly; i : natural; z : Partition;
! 70: d : Matrix ) return Poly is
! 71:
! 72: -- DESCRIPTION :
! 73: -- Returns the dense representation of the given polynomial p, known as
! 74: -- p(i) of some polynomial system, of the given m-homogeneous structure.
! 75: -- The returned polynomial has all its coefficients equal to one.
! 76:
! 77: res : Poly := Null_Poly;
! 78: maxdegs,accu : Standard_Natural_Vectors.Vector(d'range(1));
! 79:
! 80: procedure Generate_Monomials
! 81: ( k : in natural; max : in Standard_Natural_Vectors.Vector;
! 82: acc : in out Standard_Natural_Vectors.Vector) is
! 83:
! 84: -- DESCRIPTION :
! 85: -- Generates all monomials with exponents in a box, with upper bounds
! 86: -- given by the vector max. Every monomial that respects the
! 87: -- m-homogeneous structure will be added to the result res.
! 88: -- The current unknown is indicated by the parameter k.
! 89:
! 90: t : Term;
! 91:
! 92: begin
! 93: for j in 0..max(k) loop
! 94: acc(k) := j;
! 95: t.cf := Create(1.0);
! 96: t.dg := new Standard_Natural_Vectors.Vector'(acc);
! 97: if Admitted(t,i,z,d)
! 98: then Add(res,t);
! 99: end if;
! 100: Clear(t);
! 101: if k < max'last
! 102: then Generate_Monomials(k+1,max,acc);
! 103: end if;
! 104: end loop;
! 105: end Generate_Monomials;
! 106:
! 107: begin
! 108: for j in maxdegs'range loop
! 109: maxdegs(j) := Degree(p,j); accu(j) := 0;
! 110: end loop;
! 111: Generate_Monomials(1,maxdegs,accu);
! 112: return res;
! 113: end Dense_Representation;
! 114:
! 115: function Evaluate ( t : Term; x : Standard_Complex_Vectors.Vector )
! 116: return Complex_Number is
! 117:
! 118: res : Complex_Number := Create(1.0);
! 119:
! 120: begin
! 121: for i in x'range loop
! 122: for k in 1..t.dg(i) loop
! 123: res := res*x(i);
! 124: end loop;
! 125: end loop;
! 126: return res;
! 127: end Evaluate;
! 128:
! 129: procedure Interpolation_System
! 130: ( p : in Poly; b : in natural; sols : in Solution_List;
! 131: mat : out Standard_Complex_Matrices.Matrix;
! 132: rhs : out Standard_Complex_Vectors.Vector ) is
! 133:
! 134: -- DESCRIPTION :
! 135: -- Returns the matrix and right hand side vector of the linear system
! 136: -- that expresses the interpolationg conditions.
! 137:
! 138: m : Standard_Complex_Matrices.Matrix(1..b,1..b);
! 139: v : Standard_Complex_Vectors.Vector(1..b);
! 140: cnt : natural := 0;
! 141:
! 142: procedure Scan_Term ( t : in Term; continue : out boolean ) is
! 143:
! 144: tmp : Solution_List := sols;
! 145: s : Link_to_Solution;
! 146:
! 147: begin
! 148: for row in m'range(1) loop -- fill in column indicated by cnt
! 149: s := Head_Of(tmp);
! 150: if cnt = 0
! 151: then v(row) := -t.cf*Evaluate(t,s.v);
! 152: elsif cnt <= b
! 153: then m(row,cnt) := Evaluate(t,s.v);
! 154: else v(row) := v(row) - t.cf*Evaluate(t,s.v);
! 155: end if;
! 156: tmp := Tail_Of(tmp);
! 157: end loop;
! 158: cnt := cnt + 1;
! 159: continue := true;
! 160: end Scan_Term;
! 161: procedure Scan_Poly is new Visiting_Iterator(Scan_Term);
! 162:
! 163: begin
! 164: Scan_Poly(p);
! 165: mat := m; rhs := v;
! 166: end Interpolation_System;
! 167:
! 168: procedure Interpolation_System
! 169: ( p : in Poly; i,b : in natural; sols : in Solution_List;
! 170: mat : out Standard_Complex_Matrices.Matrix;
! 171: rhs : out Standard_Complex_Vectors.Vector ) is
! 172:
! 173: -- DESCRIPTION :
! 174: -- Returns the matrix and right hand side vector of the linear system
! 175: -- that expresses the interpolationg conditions. Monomials with degree
! 176: -- one in x_i will not appear in the interpolation matrix.
! 177:
! 178: m : Standard_Complex_Matrices.Matrix(1..b,1..b);
! 179: v : Standard_Complex_Vectors.Vector(1..b);
! 180: cnt : natural := 0;
! 181:
! 182: procedure Scan_Term ( t : in Term; continue : out boolean ) is
! 183:
! 184: tmp : Solution_List := sols;
! 185: s : Link_to_Solution;
! 186:
! 187: begin
! 188: for row in m'range(1) loop -- fill in column indicated by cnt
! 189: s := Head_Of(tmp);
! 190: if cnt = 0
! 191: then v(row) := -t.cf*Evaluate(t,s.v);
! 192: else if cnt <= b and t.dg(i) /= 1
! 193: then m(row,cnt) := Evaluate(t,s.v);
! 194: else v(row) := v(row) - t.cf*Evaluate(t,s.v);
! 195: end if;
! 196: end if;
! 197: tmp := Tail_Of(tmp);
! 198: end loop;
! 199: if cnt = 0 or t.dg(i) /= 1
! 200: then cnt := cnt + 1;
! 201: end if;
! 202: continue := true;
! 203: end Scan_Term;
! 204: procedure Scan_Poly is new Visiting_Iterator(Scan_Term);
! 205:
! 206: begin
! 207: Scan_Poly(p);
! 208: mat := m; rhs := v;
! 209: end Interpolation_System;
! 210:
! 211: procedure Construct_Interpolant ( p : in out Poly;
! 212: cv : in Standard_Complex_Vectors.Vector ) is
! 213:
! 214: -- DESCRIPTION :
! 215: -- With the coefficients of its monomials, the interpolant will be
! 216: -- constructed.
! 217:
! 218: cnt : natural := 0;
! 219:
! 220: procedure Scan_Term ( t : in out Term; continue : out boolean ) is
! 221: begin
! 222: if cnt > cv'last
! 223: then continue := false;
! 224: else if cnt >= cv'first
! 225: then t.cf := cv(cnt);
! 226: end if;
! 227: cnt := cnt + 1;
! 228: continue := true;
! 229: end if;
! 230: end Scan_Term;
! 231: procedure Scan_Poly is new Changing_Iterator(Scan_Term);
! 232:
! 233: begin
! 234: Scan_Poly(p);
! 235: end Construct_Interpolant;
! 236:
! 237: procedure Construct_Interpolant ( p : in out Poly; i : in natural;
! 238: cv : in Standard_Complex_Vectors.Vector ) is
! 239:
! 240: -- DESCRIPTION :
! 241: -- With the coefficients of its monomials, the interpolant will be
! 242: -- constructed. Monomials that have degree one in x_i will be ignored.
! 243:
! 244: cnt : natural := 0;
! 245:
! 246: procedure Scan_Term ( t : in out Term; continue : out boolean ) is
! 247: begin
! 248: if cnt > cv'last
! 249: then continue := false;
! 250: else if cnt = 0 or t.dg(i) /= 1
! 251: then if cnt >= cv'first
! 252: then t.cf := cv(cnt);
! 253: end if;
! 254: cnt := cnt + 1;
! 255: end if;
! 256: continue := true;
! 257: end if;
! 258: end Scan_Term;
! 259: procedure Scan_Poly is new Changing_Iterator(Scan_Term);
! 260:
! 261: begin
! 262: Scan_Poly(p);
! 263: end Construct_Interpolant;
! 264:
! 265: function Interpolate ( p : Poly; b : natural; sols : Solution_List )
! 266: return Poly is
! 267:
! 268: -- DESCRIPTION :
! 269: -- Constructs the interpolating polynomial with the same monomial
! 270: -- structure as the given polynomial.
! 271:
! 272: res : Poly := Complex_Randomize(p);
! 273: mat : Standard_Complex_Matrices.Matrix(1..b,1..b);
! 274: rhs : Standard_Complex_Vectors.Vector(1..b);
! 275: ipvt : Standard_Natural_Vectors.Vector(1..b);
! 276: info : integer;
! 277:
! 278: begin
! 279: Interpolation_System(res,b,sols,mat,rhs);
! 280: for i in ipvt'range loop
! 281: ipvt(i) := i;
! 282: end loop;
! 283: lufac(mat,b,ipvt,info);
! 284: if info = 0
! 285: then lusolve(mat,b,ipvt,rhs);
! 286: end if;
! 287: Construct_Interpolant(res,rhs);
! 288: return res;
! 289: end Interpolate;
! 290:
! 291: function Interpolate ( p : Poly; i,b : natural; sols : Solution_List )
! 292: return Poly is
! 293:
! 294: -- DESCRIPTION :
! 295: -- Constructs the interpolating polynomial with the same monomial
! 296: -- structure as the given polynomial. The monomials that have degree
! 297: -- one in x_i will not appear in the interpolation matrix.
! 298:
! 299: res : Poly := Complex_Randomize(p);
! 300: mat : Standard_Complex_Matrices.Matrix(1..b,1..b);
! 301: rhs : Standard_Complex_Vectors.Vector(1..b);
! 302: ipvt : Standard_Natural_Vectors.Vector(1..b);
! 303: info : integer;
! 304:
! 305: begin
! 306: Interpolation_System(res,i,b,sols,mat,rhs);
! 307: for j in ipvt'range loop
! 308: ipvt(j) := j;
! 309: end loop;
! 310: lufac(mat,b,ipvt,info);
! 311: if info = 0
! 312: then lusolve(mat,b,ipvt,rhs);
! 313: end if;
! 314: Construct_Interpolant(res,i,rhs);
! 315: return res;
! 316: end Interpolate;
! 317:
! 318: -- TARGET ROUTINES :
! 319:
! 320: function Dense_Representation
! 321: ( p : Poly_Sys; z : partition ) return Poly_Sys is
! 322:
! 323: d : constant Standard_Integer_Matrices.Matrix := Degree_Table(p,z);
! 324:
! 325: begin
! 326: return Dense_Representation(p,z,d);
! 327: end Dense_Representation;
! 328:
! 329: function Dense_Representation
! 330: ( p : Poly_Sys; z : partition; d : Matrix ) return Poly_Sys is
! 331:
! 332: res : Poly_Sys(d'range(1));
! 333:
! 334: begin
! 335: for i in res'range loop
! 336: if p(i) /= Null_Poly
! 337: then res(i) := Dense_Representation(p(i),i,z,d);
! 338: end if;
! 339: end loop;
! 340: return res;
! 341: end Dense_Representation;
! 342:
! 343: function Independent_Representation ( p : Poly_Sys ) return Poly_Sys is
! 344:
! 345: res : Poly_Sys(p'range);
! 346: it : Term;
! 347:
! 348: begin
! 349: Copy(p,res);
! 350: for i in res'range loop
! 351: if p(i) /= Null_Poly
! 352: then it := Initial_Term(res(i));
! 353: for j in res'range loop
! 354: if j /= i and then (p(j) /= Null_Poly)
! 355: then Eliminate_Term(res(j),it.dg);
! 356: end if;
! 357: end loop;
! 358: end if;
! 359: end loop;
! 360: return res;
! 361: end Independent_Representation;
! 362:
! 363: function Independent_Roots ( p : Poly_Sys ) return natural is
! 364:
! 365: ntp : natural := 0;
! 366: min : natural := ntp;
! 367:
! 368: begin
! 369: for i in p'first..p'last loop
! 370: if p(i) /= Null_Poly
! 371: then ntp := Number_of_Terms(p(i));
! 372: if min = 0
! 373: then min := ntp;
! 374: elsif ntp < min
! 375: then min := ntp;
! 376: end if;
! 377: end if;
! 378: end loop;
! 379: if min = 0
! 380: then return 0;
! 381: else return (min-1);
! 382: end if;
! 383: end Independent_Roots;
! 384:
! 385: function Number_of_Terms ( p : Poly; i : natural ) return natural is
! 386:
! 387: -- DESCRIPTION :
! 388: -- Returns the number of monomials of p, without those monomials that
! 389: -- have degree one in x_i.
! 390:
! 391: cnt : natural := 0;
! 392:
! 393: procedure Count_Term ( t : in Term; continue : out boolean ) is
! 394: begin
! 395: if t.dg(i) /= 1
! 396: then cnt := cnt+1;
! 397: end if;
! 398: continue := true;
! 399: end Count_Term;
! 400: procedure Count_Terms is new Visiting_Iterator(Count_Term);
! 401:
! 402: begin
! 403: Count_Terms(p);
! 404: return cnt;
! 405: end Number_of_Terms;
! 406:
! 407: function Independent_Roots ( p : Poly_Sys; i : natural ) return natural is
! 408:
! 409: ntp : natural := 0;
! 410: min : natural := ntp;
! 411:
! 412: begin
! 413: for j in p'first..p'last loop
! 414: if p(j) /= Null_Poly
! 415: then ntp := Number_of_Terms(p(j),i);
! 416: if min = 0
! 417: then min := ntp;
! 418: elsif ntp < min
! 419: then min := ntp;
! 420: end if;
! 421: end if;
! 422: end loop;
! 423: if min = 0
! 424: then return 0;
! 425: else return (min-1);
! 426: end if;
! 427: end Independent_Roots;
! 428:
! 429: function Interpolate ( p : Poly_Sys; b : natural; sols : Solution_List )
! 430: return Poly_Sys is
! 431:
! 432: res : Poly_Sys(p'range);
! 433:
! 434: begin
! 435: for i in p'range loop
! 436: if p(i) /= Null_Poly
! 437: then res(i) := Interpolate(p(i),b,sols);
! 438: end if;
! 439: end loop;
! 440: return res;
! 441: end Interpolate;
! 442:
! 443: function Interpolate ( p : Poly_Sys; i,b : natural; sols : Solution_List )
! 444: return Poly_Sys is
! 445:
! 446: res : Poly_Sys(p'range);
! 447:
! 448: begin
! 449: for j in p'range loop
! 450: if p(j) /= Null_Poly
! 451: then res(j) := Interpolate(p(j),i,b,sols);
! 452: end if;
! 453: end loop;
! 454: return res;
! 455: end Interpolate;
! 456:
! 457: end Interpolating_Homotopies;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to interpolating_homotopies.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 / Product