Annotation of OpenXM_contrib/PHC/Ada/Schubert/numeric_minor_equations.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
! 2: with Standard_Natural_Vectors;
! 3: with Standard_Complex_Poly_Functions; use Standard_Complex_Poly_Functions;
! 4: with Standard_Complex_Poly_SysFun; use Standard_Complex_Poly_SysFun;
! 5: with Brackets; use Brackets;
! 6: with Plane_Representations; use Plane_Representations;
! 7: with Evaluated_Minors; use Evaluated_Minors;
! 8: with Symbolic_Minor_Equations; use Symbolic_Minor_Equations;
! 9:
! 10: package body Numeric_Minor_Equations is
! 11:
! 12: tol : constant double_float := 10.0**(-10);
! 13:
! 14: -- EXPANDING ACCORDING A BRACKET MONOMIAL :
! 15:
! 16: function Expanded_Minors
! 17: ( cffmat : Standard_Floating_Matrices.Matrix;
! 18: polmat : Standard_Complex_Poly_Matrices.Matrix;
! 19: bm : Bracket_Monomial ) return Poly is
! 20:
! 21: res : Poly := Null_Poly;
! 22: first : boolean := true;
! 23:
! 24: procedure Visit_Bracket ( b : in Bracket; continue : out boolean ) is
! 25:
! 26: factor : double_float;
! 27: minor : Poly;
! 28:
! 29: begin
! 30: if first
! 31: then declare
! 32: bb : Bracket(b'first..b'last-1) := b(b'first+1..b'last);
! 33: begin
! 34: factor := Determinant(cffmat,bb);
! 35: end;
! 36: first := false;
! 37: else minor := Expanded_Minor(polmat,b);
! 38: if (minor /= Null_Poly) and (abs(factor) > tol)
! 39: then Mul(minor,Create(factor));
! 40: Add(res,minor);
! 41: end if;
! 42: Clear(minor);
! 43: end if;
! 44: continue := true;
! 45: end Visit_Bracket;
! 46: procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
! 47:
! 48: begin
! 49: Visit_Brackets(bm);
! 50: return res;
! 51: end Expanded_Minors;
! 52:
! 53: function Expanded_Minors
! 54: ( cffmat : Standard_Complex_Matrices.Matrix;
! 55: polmat : Standard_Complex_Poly_Matrices.Matrix;
! 56: bm : Bracket_Monomial ) return Poly is
! 57:
! 58: res : Poly := Null_Poly;
! 59: first : boolean := true;
! 60: factor : Complex_Number;
! 61:
! 62: procedure Visit_Bracket ( b : in Bracket; continue : out boolean ) is
! 63:
! 64: minor : Poly;
! 65:
! 66: begin
! 67: if first
! 68: then declare
! 69: bb : Bracket(b'first..b'last-1) := b(b'first+1..b'last);
! 70: begin
! 71: factor := Determinant(cffmat,bb);
! 72: end;
! 73: first := false;
! 74: else minor := Expanded_Minor(polmat,b);
! 75: if (minor /= Null_Poly) and (AbsVal(factor) > tol)
! 76: then Mul(minor,factor);
! 77: Add(res,minor);
! 78: end if;
! 79: Clear(minor);
! 80: end if;
! 81: continue := true;
! 82: end Visit_Bracket;
! 83: procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
! 84:
! 85: begin
! 86: Visit_Brackets(bm);
! 87: return res;
! 88: end Expanded_Minors;
! 89:
! 90: function Expanded_Minors
! 91: ( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
! 92: bm : Bracket_Monomial ) return Poly is
! 93:
! 94: res : Poly := Null_Poly;
! 95: first : boolean := true;
! 96: factor : Poly;
! 97:
! 98: procedure Visit_Bracket ( b : in Bracket; continue : out boolean ) is
! 99:
! 100: minor : Poly;
! 101:
! 102: begin
! 103: if first
! 104: then declare
! 105: bb : Bracket(b'first..b'last-1) := b(b'first+1..b'last);
! 106: begin
! 107: factor := Expanded_Minor(cntmat,bb);
! 108: end;
! 109: first := false;
! 110: else minor := Expanded_Minor(polmat,b);
! 111: if (minor /= Null_Poly) and (factor /= Null_Poly)
! 112: then Mul(minor,factor);
! 113: Add(res,minor);
! 114: end if;
! 115: Clear(factor); Clear(minor);
! 116: end if;
! 117: continue := true;
! 118: end Visit_Bracket;
! 119: procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
! 120:
! 121: begin
! 122: Visit_Brackets(bm);
! 123: return res;
! 124: end Expanded_Minors;
! 125:
! 126: function Lifted_Expanded_Minors
! 127: ( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
! 128: bm : Bracket_Monomial ) return Poly is
! 129:
! 130: res : Poly := Null_Poly;
! 131: first : boolean := true;
! 132: factor : Poly;
! 133:
! 134: procedure Visit_Bracket ( b : in Bracket; continue : out boolean ) is
! 135:
! 136: minor,extmin : Poly;
! 137:
! 138: begin
! 139: if first
! 140: then declare
! 141: bb : Bracket(b'first..b'last-1) := b(b'first+1..b'last);
! 142: begin
! 143: factor := Expanded_Minor(cntmat,bb);
! 144: end;
! 145: first := false;
! 146: else minor := Expanded_Minor(polmat,b);
! 147: if (minor /= Null_Poly) and (factor /= Null_Poly)
! 148: then extmin := Extend_Zero_Lifting(minor);
! 149: Mul(extmin,factor);
! 150: Add(res,extmin);
! 151: end if;
! 152: Clear(factor); Clear(minor); Clear(extmin);
! 153: end if;
! 154: continue := true;
! 155: end Visit_Bracket;
! 156: procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
! 157:
! 158: begin
! 159: Visit_Brackets(bm);
! 160: return res;
! 161: end Lifted_Expanded_Minors;
! 162:
! 163: -- EXPANDING ACCORDING A BRACKET POLYNOMIAL :
! 164:
! 165: function Expanded_Minors
! 166: ( cffmat : Standard_Floating_Matrices.Matrix;
! 167: polmat : Standard_Complex_Poly_Matrices.Matrix;
! 168: bp : Bracket_Polynomial ) return Poly is
! 169:
! 170: res : Poly := Null_Poly;
! 171:
! 172: procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
! 173:
! 174: minor : Poly := Expanded_Minors(cffmat,polmat,t.monom);
! 175:
! 176: begin
! 177: if REAL_PART(t.coeff) > 0.0
! 178: then Add(res,minor);
! 179: else Sub(res,minor);
! 180: end if;
! 181: Clear(minor);
! 182: continue := true;
! 183: end Visit_Term;
! 184: procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
! 185:
! 186: begin
! 187: Visit_Terms(bp);
! 188: return res;
! 189: end Expanded_Minors;
! 190:
! 191: function Expanded_Minors
! 192: ( cffmat : Standard_Complex_Matrices.Matrix;
! 193: polmat : Standard_Complex_Poly_Matrices.Matrix;
! 194: bp : Bracket_Polynomial ) return Poly is
! 195:
! 196: res : Poly := Null_Poly;
! 197:
! 198: procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
! 199:
! 200: minor : Poly := Expanded_Minors(cffmat,polmat,t.monom);
! 201:
! 202: begin
! 203: if REAL_PART(t.coeff) > 0.0
! 204: then Add(res,minor);
! 205: else Sub(res,minor);
! 206: end if;
! 207: Clear(minor);
! 208: continue := true;
! 209: end Visit_Term;
! 210: procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
! 211:
! 212: begin
! 213: Visit_Terms(bp);
! 214: return res;
! 215: end Expanded_Minors;
! 216:
! 217: function Expanded_Minors
! 218: ( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
! 219: bp : Bracket_Polynomial ) return Poly is
! 220:
! 221: res : Poly := Null_Poly;
! 222:
! 223: procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
! 224:
! 225: minor : Poly := Expanded_Minors(cntmat,polmat,t.monom);
! 226:
! 227: begin
! 228: if REAL_PART(t.coeff) > 0.0
! 229: then Add(res,minor);
! 230: else Sub(res,minor);
! 231: end if;
! 232: Clear(minor);
! 233: continue := true;
! 234: end Visit_Term;
! 235: procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
! 236:
! 237: begin
! 238: Visit_Terms(bp);
! 239: return res;
! 240: end Expanded_Minors;
! 241:
! 242: function Lifted_Expanded_Minors
! 243: ( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
! 244: bp : Bracket_Polynomial ) return Poly is
! 245:
! 246: res : Poly := Null_Poly;
! 247:
! 248: procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
! 249:
! 250: minor : Poly := Lifted_Expanded_Minors(cntmat,polmat,t.monom);
! 251:
! 252: begin
! 253: if REAL_PART(t.coeff) > 0.0
! 254: then Add(res,minor);
! 255: else Sub(res,minor);
! 256: end if;
! 257: Clear(minor);
! 258: continue := true;
! 259: end Visit_Term;
! 260: procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
! 261:
! 262: begin
! 263: Visit_Terms(bp);
! 264: return res;
! 265: end Lifted_Expanded_Minors;
! 266:
! 267: -- EXPANDING TO CONSTRUCT POLYNOMIAL SYSTEMS :
! 268:
! 269: function Expanded_Minors
! 270: ( cffmat : Standard_Floating_Matrices.Matrix;
! 271: polmat : Standard_Complex_Poly_Matrices.Matrix;
! 272: bs : Bracket_System ) return Poly_Sys is
! 273:
! 274: res : Poly_Sys(bs'first+1..bs'last);
! 275:
! 276: begin
! 277: for i in res'range loop
! 278: res(i) := Expanded_Minors(cffmat,polmat,bs(i));
! 279: end loop;
! 280: return res;
! 281: end Expanded_Minors;
! 282:
! 283: function Expanded_Minors
! 284: ( cffmat : Standard_Complex_Matrices.Matrix;
! 285: polmat : Standard_Complex_Poly_Matrices.Matrix;
! 286: bs : Bracket_System ) return Poly_Sys is
! 287:
! 288: res : Poly_Sys(bs'first+1..bs'last);
! 289:
! 290: begin
! 291: for i in res'range loop
! 292: res(i) := Expanded_Minors(cffmat,polmat,bs(i));
! 293: end loop;
! 294: return res;
! 295: end Expanded_Minors;
! 296:
! 297: function Expanded_Minors
! 298: ( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
! 299: bs : Bracket_System ) return Poly_Sys is
! 300:
! 301: res : Poly_Sys(bs'first+1..bs'last);
! 302:
! 303: begin
! 304: for i in res'range loop
! 305: res(i) := Expanded_Minors(cntmat,polmat,bs(i));
! 306: end loop;
! 307: return res;
! 308: end Expanded_Minors;
! 309:
! 310: function Lifted_Expanded_Minors
! 311: ( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
! 312: bs : Bracket_System ) return Poly_Sys is
! 313:
! 314: res : Poly_Sys(bs'first+1..bs'last);
! 315:
! 316: begin
! 317: for i in res'range loop
! 318: res(i) := Lifted_Expanded_Minors(cntmat,polmat,bs(i));
! 319: end loop;
! 320: return res;
! 321: end Lifted_Expanded_Minors;
! 322:
! 323: function Evaluate ( p : Poly; x : Standard_Complex_Matrices.Matrix )
! 324: return Complex_Number is
! 325:
! 326: xv : constant Standard_Complex_Vectors.Vector := Vector_Rep(x);
! 327:
! 328: begin
! 329: return Eval(p,xv);
! 330: end Evaluate;
! 331:
! 332: function Evaluate ( p : Poly_Sys; x : Standard_Complex_Matrices.Matrix )
! 333: return Standard_Complex_Vectors.Vector is
! 334:
! 335: xv : constant Standard_Complex_Vectors.Vector := Vector_Rep(x);
! 336:
! 337: begin
! 338: return Eval(p,xv);
! 339: end Evaluate;
! 340:
! 341: procedure Embed ( t : in out Term ) is
! 342:
! 343: dg : Degrees
! 344: := new Standard_Natural_Vectors.Vector(t.dg'first..t.dg'last+1);
! 345:
! 346: begin
! 347: dg(t.dg'range) := t.dg.all;
! 348: dg(dg'last) := 0;
! 349: Clear(t.dg);
! 350: t.dg := dg;
! 351: end Embed;
! 352:
! 353: function Embed ( t : Term ) return Term is
! 354:
! 355: res : Term;
! 356:
! 357: begin
! 358: res.dg := new Standard_Natural_Vectors.Vector(t.dg'first..t.dg'last+1);
! 359: res.dg(t.dg'range) := t.dg.all;
! 360: res.dg(res.dg'last) := 0;
! 361: res.cf := t.cf;
! 362: return res;
! 363: end Embed;
! 364:
! 365: procedure Embed ( p : in out Poly ) is
! 366:
! 367: res : Poly := Null_Poly;
! 368:
! 369: procedure Embed_Term ( t : in Term; continue : out boolean ) is
! 370:
! 371: et : Term := Embed(t);
! 372:
! 373: begin
! 374: Add(res,et);
! 375: Clear(et);
! 376: continue := true;
! 377: end Embed_Term;
! 378: procedure Embed_Terms is new Visiting_Iterator(Embed_Term);
! 379:
! 380: begin
! 381: Embed_Terms(p);
! 382: Clear(p);
! 383: p := res;
! 384: end Embed;
! 385:
! 386: procedure Embed ( p : in out Poly_Sys ) is
! 387: begin
! 388: for i in p'range loop
! 389: Embed(p(i));
! 390: end loop;
! 391: end Embed;
! 392:
! 393: procedure Embed ( m : in out Standard_Complex_Poly_Matrices.Matrix ) is
! 394: begin
! 395: for i in m'range(1) loop
! 396: for j in m'range(2) loop
! 397: if m(i,j) /= Null_Poly
! 398: then Embed(m(i,j));
! 399: end if;
! 400: end loop;
! 401: end loop;
! 402: end Embed;
! 403:
! 404: function Linear_Homotopy ( target,start : Poly ) return Poly is
! 405:
! 406: res : Poly := Null_Poly;
! 407:
! 408: procedure Embed_Target_Term ( t : in Term; continue : out boolean ) is
! 409:
! 410: et : Term := Embed(t);
! 411:
! 412: begin
! 413: et.dg(et.dg'last) := 1;
! 414: Add(res,et);
! 415: Clear(et);
! 416: continue := true;
! 417: end Embed_Target_Term;
! 418: procedure Embed_Target_Terms is new Visiting_Iterator(Embed_Target_Term);
! 419:
! 420: procedure Embed_Start_Term ( t : in Term; continue : out boolean ) is
! 421:
! 422: et : Term := Embed(t);
! 423:
! 424: begin
! 425: Add(res,et);
! 426: et.dg(et.dg'last) := 1;
! 427: Sub(res,et);
! 428: Clear(et);
! 429: continue := true;
! 430: end Embed_Start_Term;
! 431: procedure Embed_Start_Terms is new Visiting_Iterator(Embed_Start_Term);
! 432:
! 433: begin
! 434: Embed_Target_Terms(target);
! 435: Embed_Start_Terms(start);
! 436: return res;
! 437: end Linear_Homotopy;
! 438:
! 439: function Linear_Interpolation
! 440: ( target,start : Poly; k : natural ) return Poly is
! 441:
! 442: res : Poly := Null_Poly;
! 443:
! 444: procedure Embed_Target_Term ( t : in Term; continue : out boolean ) is
! 445:
! 446: et : Term;
! 447:
! 448: begin
! 449: Copy(t,et);
! 450: et.dg(k) := et.dg(k) + 1; -- multiply with t
! 451: Add(res,et);
! 452: Clear(et);
! 453: continue := true;
! 454: end Embed_Target_Term;
! 455: procedure Embed_Target_Terms is new Visiting_Iterator(Embed_Target_Term);
! 456:
! 457: procedure Embed_Start_Term ( t : in Term; continue : out boolean ) is
! 458:
! 459: et : Term;
! 460:
! 461: begin
! 462: Copy(t,et);
! 463: Add(res,et); -- res := res + et
! 464: et.dg(k) := et.dg(k) + 1; -- multiply with t
! 465: Sub(res,et); -- res := res + et - t*et
! 466: Clear(et);
! 467: continue := true;
! 468: end Embed_Start_Term;
! 469: procedure Embed_Start_Terms is new Visiting_Iterator(Embed_Start_Term);
! 470:
! 471: begin
! 472: Embed_Target_Terms(target);
! 473: Embed_Start_Terms(start);
! 474: return res;
! 475: end Linear_Interpolation;
! 476:
! 477: procedure Divide_Common_Factor ( p : in out Poly; k : in natural ) is
! 478:
! 479: first : boolean := true;
! 480: min : natural;
! 481:
! 482: procedure Min_Power ( t : in Term; continue : out boolean ) is
! 483: begin
! 484: if first
! 485: then first := false;
! 486: min := t.dg(k); -- initialize minimal power
! 487: else if t.dg(k) < min
! 488: then min := t.dg(k);
! 489: end if;
! 490: end if;
! 491: continue := true;
! 492: end Min_Power;
! 493: procedure Find_Min_Power is new Visiting_Iterator(Min_Power);
! 494:
! 495: procedure Divide ( t : in out Term; continue : out boolean ) is
! 496: begin
! 497: t.dg(k) := t.dg(k) - min;
! 498: continue := true;
! 499: end Divide;
! 500: procedure Divide_Min_Power is new Changing_Iterator(Divide);
! 501:
! 502: begin
! 503: Find_Min_Power(p);
! 504: if min > 0
! 505: then Divide_Min_Power(p);
! 506: end if;
! 507: end Divide_Common_Factor;
! 508:
! 509: end Numeric_Minor_Equations;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to numeric_minor_equations.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Schubert