Annotation of OpenXM_contrib/PHC/Ada/Schubert/sagbi_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; use Standard_Natural_Vectors;
! 4: with Standard_Complex_Poly_Functions; use Standard_Complex_Poly_Functions;
! 5: with Brackets; use Brackets;
! 6: with Bracket_Monomials; use Bracket_Monomials;
! 7: with Bracket_Polynomials; use Bracket_Polynomials;
! 8: with Straightening_Syzygies; use Straightening_Syzygies;
! 9: with Bracket_Expansions; use Bracket_Expansions;
! 10: with Evaluated_Minors; use Evaluated_Minors;
! 11:
! 12: package body SAGBI_Homotopies is
! 13:
! 14: function Coordinatize_Hexadecimal ( b : Bracket ) return natural is
! 15:
! 16: -- DESCRIPTION :
! 17: -- Returns the hexadecimal expansion of the entries in the bracket.
! 18:
! 19: res : natural := 0;
! 20:
! 21: begin
! 22: for i in b'range loop
! 23: res := res*16 + b(i);
! 24: end loop;
! 25: return res;
! 26: end Coordinatize_Hexadecimal;
! 27:
! 28: function Unsigned ( i : integer ) return natural is
! 29:
! 30: -- DESCRIPTION :
! 31: -- Returns the unsigned integer.
! 32:
! 33: n : natural;
! 34:
! 35: begin
! 36: if i < 0
! 37: then n := -i;
! 38: else n := i;
! 39: end if;
! 40: return n;
! 41: end Unsigned;
! 42:
! 43: function Bracketize_Hexadecimal ( n,d : natural ) return Bracket is
! 44:
! 45: -- DESCRIPTION :
! 46: -- Returns the d-bracket from the hexadecimal expansion n.
! 47:
! 48: res : Bracket(1..d);
! 49: nn : natural := n;
! 50:
! 51: begin
! 52: for i in reverse 1..d loop
! 53: res(i) := nn mod 16;
! 54: nn := nn/16;
! 55: end loop;
! 56: return res;
! 57: end Bracketize_Hexadecimal;
! 58:
! 59: function Coordinatize ( p : Bracket_Polynomial ) return Bracket_Polynomial is
! 60:
! 61: -- DESCRIPTION :
! 62: -- Replaces the first bracket in every monomial by the decimal expansion.
! 63:
! 64: res : Bracket_Polynomial;
! 65:
! 66: procedure Coordinatize_Term ( t : in Bracket_Term; cont1 : out boolean ) is
! 67:
! 68: first,second : boolean;
! 69: bm : Bracket_Monomial;
! 70: bt : Bracket_Term;
! 71:
! 72: procedure Coordinatize_Bracket ( b : in Bracket; cont2 : out boolean ) is
! 73: begin
! 74: if first
! 75: then bt.coeff := Create(double_float(Coordinatize_Hexadecimal(b)));
! 76: first := false;
! 77: second := true;
! 78: elsif second
! 79: then bm := Create(b);
! 80: else Multiply(bm,b);
! 81: end if;
! 82: cont2 := true;
! 83: end Coordinatize_Bracket;
! 84: procedure Coordinatize_Brackets is
! 85: new Enumerate_Brackets(Coordinatize_Bracket);
! 86:
! 87: begin
! 88: first := true; second := false;
! 89: Coordinatize_Brackets(t.monom);
! 90: bt.monom := bm;
! 91: if REAL_PART(t.coeff) < 0.0
! 92: then Min(res,bt);
! 93: else Add(res,bt);
! 94: end if;
! 95: cont1 := true;
! 96: end Coordinatize_Term;
! 97: procedure Coordinatize_Terms is new Enumerate_Terms(Coordinatize_Term);
! 98:
! 99: begin
! 100: Coordinatize_Terms(p);
! 101: return res;
! 102: end Coordinatize;
! 103:
! 104: procedure Divide ( p : in out Poly; w : in natural ) is
! 105:
! 106: -- DESCRIPTION :
! 107: -- Divides the polynomial by t^w.
! 108:
! 109: procedure Divide_Term ( t : in out Term; continue : out boolean ) is
! 110: begin
! 111: t.dg(t.dg'last) := t.dg(t.dg'last) - w;
! 112: continue := true;
! 113: end Divide_Term;
! 114: procedure Divide_Terms is new Changing_Iterator(Divide_Term);
! 115:
! 116: begin
! 117: Divide_Terms(p);
! 118: end Divide;
! 119:
! 120: function Weight ( e : Standard_Natural_Vectors.Vector; n,d : natural )
! 121: return natural is
! 122:
! 123: -- DESCRIPTION :
! 124: -- Returns the weight of the exponent vector for the localization that
! 125: -- takes the d-by-d identitity matrix in the lower-right of the d-plane.
! 126: -- The lifting recipe is xij*t^((i-1)*(d-j)).
! 127:
! 128: res : natural := 0;
! 129: jmp,ind : natural;
! 130:
! 131: begin
! 132: for j in 1..d loop
! 133: jmp := d-j;
! 134: for i in 1..n-d loop
! 135: ind := (i-1)*d + j;
! 136: if e(ind) > 0
! 137: then res := res + (i-1)*jmp;
! 138: end if;
! 139: end loop;
! 140: end loop;
! 141: return res;
! 142: end Weight;
! 143:
! 144: function Weight ( locmap : Standard_Natural_Matrices.Matrix;
! 145: e : Standard_Natural_Vectors.Vector ) return natural is
! 146:
! 147: -- DESCRIPTION :
! 148: -- Returns the weight of the exponent vector as xij*t^((i-1)*(d-j))
! 149: -- for the localization pattern in locmap.
! 150:
! 151: res : natural := 0;
! 152: d : constant natural := locmap'length(2);
! 153: jmp,ind : natural;
! 154:
! 155: begin
! 156: ind := 0;
! 157: for i in locmap'range(1) loop
! 158: for j in locmap'range(2) loop
! 159: jmp := d-j;
! 160: if locmap(i,j) = 2
! 161: then ind := ind+1;
! 162: if e(ind) > 0
! 163: then res := res + (i-1)*jmp;
! 164: end if;
! 165: end if;
! 166: end loop;
! 167: end loop;
! 168: return res;
! 169: end Weight;
! 170:
! 171: function Lift ( p : Poly; n,d : natural ) return Poly is
! 172:
! 173: -- DESCRIPTION :
! 174: -- Returns the lifted polynomial, where the xij is lifted according
! 175: -- to xij*t^((i-1)*(d-j)). The lowest powers of t are divided out.
! 176: -- The d-by-d identity matrix is the lower-right of the d-plane.
! 177:
! 178: res : Poly := Null_Poly;
! 179: first : boolean := true;
! 180: minwei : natural;
! 181:
! 182: procedure Lift_Term ( t : in Term; continue : out boolean ) is
! 183:
! 184: tt : Term;
! 185: wei : natural;
! 186:
! 187: begin
! 188: tt.cf := t.cf;
! 189: tt.dg := new Standard_Natural_Vectors.Vector(1..t.dg'last+1);
! 190: tt.dg(t.dg'range) := t.dg.all;
! 191: wei := Weight(t.dg.all,n,d);
! 192: tt.dg(tt.dg'last) := wei;
! 193: Add(res,tt);
! 194: Clear(tt.dg);
! 195: if first
! 196: then minwei := wei;
! 197: first := false;
! 198: elsif wei < minwei
! 199: then minwei := wei;
! 200: end if;
! 201: continue := true;
! 202: end Lift_Term;
! 203: procedure Lift_Terms is new Visiting_Iterator(Lift_Term);
! 204:
! 205: begin
! 206: Lift_Terms(p);
! 207: if minwei /= 0
! 208: then Divide(res,minwei);
! 209: end if;
! 210: return res;
! 211: end Lift;
! 212:
! 213: function Lift ( locmap : Standard_Natural_Matrices.Matrix; p : Poly )
! 214: return Poly is
! 215:
! 216: -- DESCRIPTION :
! 217: -- Lifts p as to xij*t^((i-1)*(d-j)) and divides by the lowest powers
! 218: -- of t, respecting the localization pattern in locmap.
! 219:
! 220: res : Poly := Null_Poly;
! 221: first : boolean := true;
! 222: minwei : natural;
! 223:
! 224: procedure Lift_Term ( t : in Term; continue : out boolean ) is
! 225:
! 226: tt : Term;
! 227: wei : natural;
! 228:
! 229: begin
! 230: tt.cf := t.cf;
! 231: tt.dg := new Standard_Natural_Vectors.Vector(1..t.dg'last+1);
! 232: tt.dg(t.dg'range) := t.dg.all;
! 233: wei := Weight(locmap,t.dg.all);
! 234: tt.dg(tt.dg'last) := wei;
! 235: Add(res,tt);
! 236: Clear(tt.dg);
! 237: if first
! 238: then minwei := wei;
! 239: first := false;
! 240: elsif wei < minwei
! 241: then minwei := wei;
! 242: end if;
! 243: continue := true;
! 244: end Lift_Term;
! 245: procedure Lift_Terms is new Visiting_Iterator(Lift_Term);
! 246:
! 247: begin
! 248: Lift_Terms(p);
! 249: if minwei /= 0
! 250: then Divide(res,minwei);
! 251: end if;
! 252: return res;
! 253: end Lift;
! 254:
! 255: -- TARGET ROUTINES :
! 256:
! 257: function Lifted_Localized_Laplace_Expansion ( n,d : natural ) return Poly is
! 258:
! 259: res : Poly := Null_Poly;
! 260: p : Bracket_Polynomial := Laplace_Expansion(n,n-d);
! 261:
! 262: procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
! 263:
! 264: first : boolean := true;
! 265: cf : integer;
! 266:
! 267: procedure Visit_Bracket ( b : in Bracket; cont : out boolean ) is
! 268:
! 269: pb,lp : Poly;
! 270:
! 271: begin
! 272: if first
! 273: then cf := Coordinatize_Hexadecimal(b);
! 274: first := false;
! 275: else pb := Localized_Expand(n,d,b);
! 276: lp := Lift(pb,n,d); Clear(pb);
! 277: Mul(lp,Create(double_float(cf)));
! 278: Add(res,lp);
! 279: Clear(lp);
! 280: end if;
! 281: cont := true;
! 282: end Visit_Bracket;
! 283: procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
! 284:
! 285: begin
! 286: Visit_Brackets(t.monom);
! 287: continue := true;
! 288: end Visit_Term;
! 289: procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
! 290:
! 291: begin
! 292: Visit_Terms(p);
! 293: return res;
! 294: end Lifted_Localized_Laplace_Expansion;
! 295:
! 296: function Lifted_Localized_Laplace_Expansion
! 297: ( locmap : Standard_Natural_Matrices.Matrix ) return Poly is
! 298:
! 299: res : Poly := Null_Poly;
! 300: n : constant natural := locmap'length(1);
! 301: d : constant natural := locmap'length(2);
! 302: p : Bracket_Polynomial := Laplace_Expansion(n,n-d);
! 303:
! 304: procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
! 305:
! 306: first : boolean := true;
! 307: cf : integer;
! 308:
! 309: procedure Visit_Bracket ( b : in Bracket; cont : out boolean ) is
! 310:
! 311: pb,lp : Poly;
! 312:
! 313: begin
! 314: if first
! 315: then cf := Coordinatize_Hexadecimal(b);
! 316: first := false;
! 317: else pb := Expand(locmap,b);
! 318: Reduce_Variables(locmap,pb);
! 319: lp := Lift(locmap,pb); Clear(pb);
! 320: Mul(lp,Create(double_float(cf)));
! 321: Add(res,lp);
! 322: Clear(lp);
! 323: end if;
! 324: cont := true;
! 325: end Visit_Bracket;
! 326: procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
! 327:
! 328: begin
! 329: Visit_Brackets(t.monom);
! 330: continue := true;
! 331: end Visit_Term;
! 332: procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
! 333:
! 334: begin
! 335: Visit_Terms(p);
! 336: return res;
! 337: end Lifted_Localized_Laplace_Expansion;
! 338:
! 339: function Intersection_Coefficients
! 340: ( m : Standard_Floating_Matrices.Matrix;
! 341: c : Standard_Complex_Vectors.Vector )
! 342: return Standard_Complex_Vectors.Vector is
! 343:
! 344: res : Standard_Complex_Vectors.Vector(c'range);
! 345: nmd : constant natural := m'last(2);
! 346: ind : integer;
! 347: b : Bracket(1..nmd);
! 348:
! 349: begin
! 350: for i in c'range loop
! 351: ind := integer(REAL_PART(c(i)));
! 352: b := Bracketize_Hexadecimal(Unsigned(ind),nmd);
! 353: if ind > 0
! 354: then res(i) := Create(Determinant(m,b));
! 355: else res(i) := Create(-Determinant(m,b));
! 356: end if;
! 357: end loop;
! 358: return res;
! 359: end Intersection_Coefficients;
! 360:
! 361: function Intersection_Coefficients
! 362: ( m : Standard_Complex_Matrices.Matrix;
! 363: c : Standard_Complex_Vectors.Vector )
! 364: return Standard_Complex_Vectors.Vector is
! 365:
! 366: res : Standard_Complex_Vectors.Vector(c'range);
! 367: nmd : constant natural := m'last(2);
! 368: ind : integer;
! 369: b : Bracket(1..nmd);
! 370:
! 371: begin
! 372: for i in c'range loop
! 373: ind := integer(REAL_PART(c(i)));
! 374: b := Bracketize_Hexadecimal(Unsigned(ind),nmd);
! 375: if ind > 0
! 376: then res(i) := Determinant(m,b);
! 377: else res(i) := -Determinant(m,b);
! 378: end if;
! 379: end loop;
! 380: return res;
! 381: end Intersection_Coefficients;
! 382:
! 383: function Intersection_Condition
! 384: ( m : Standard_Floating_Matrices.Matrix; p : Poly ) return Poly is
! 385:
! 386: res : Poly := Null_Poly;
! 387: nmd : constant natural := m'last(2);
! 388:
! 389: procedure Visit_Term ( t : in Term; continue : out boolean ) is
! 390:
! 391: c : integer := integer(REAL_PART(t.cf));
! 392: b : Bracket(1..nmd) := Bracketize_Hexadecimal(Unsigned(c),nmd);
! 393: det : double_float := Determinant(m,b);
! 394: rt : Term;
! 395:
! 396: begin
! 397: if c > 0
! 398: then rt.cf := Create(det);
! 399: else rt.cf := Create(-det);
! 400: end if;
! 401: rt.dg := t.dg;
! 402: Add(res,rt);
! 403: continue := true;
! 404: end Visit_Term;
! 405: procedure Visit_Terms is new Visiting_Iterator(Visit_Term);
! 406:
! 407: begin
! 408: Visit_Terms(p);
! 409: return res;
! 410: end Intersection_Condition;
! 411:
! 412: function Intersection_Condition
! 413: ( m : Standard_Complex_Matrices.Matrix; p : Poly ) return Poly is
! 414:
! 415: res : Poly := Null_Poly;
! 416: nmd : constant natural := m'last(2);
! 417:
! 418: procedure Visit_Term ( t : in Term; continue : out boolean ) is
! 419:
! 420: c : integer := integer(REAL_PART(t.cf));
! 421: b : Bracket(1..nmd) := Bracketize_Hexadecimal(Unsigned(c),nmd);
! 422: det : Complex_Number := Determinant(m,b);
! 423: rt : Term;
! 424:
! 425: begin
! 426: if c > 0
! 427: then rt.cf := det;
! 428: else rt.cf := -det;
! 429: end if;
! 430: rt.dg := t.dg;
! 431: Add(res,rt);
! 432: continue := true;
! 433: end Visit_Term;
! 434: procedure Visit_Terms is new Visiting_Iterator(Visit_Term);
! 435:
! 436: begin
! 437: Visit_Terms(p);
! 438: return res;
! 439: end Intersection_Condition;
! 440:
! 441: end SAGBI_Homotopies;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to sagbi_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 / Schubert