Annotation of OpenXM_contrib/PHC/Ada/Schubert/bracket_expansions.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
! 2: with Standard_Natural_Vectors; use Standard_Natural_Vectors;
! 3:
! 4: package body Bracket_Expansions is
! 5:
! 6: -- AUXILIARIES :
! 7:
! 8: function Create_Term ( n,d,i,j : natural ) return Term is
! 9:
! 10: -- DESCRIPTION :
! 11: -- Returns the term that corresponds with xij.
! 12:
! 13: res : Term;
! 14:
! 15: begin
! 16: res.cf := Create(1.0);
! 17: res.dg := new Standard_Natural_Vectors.Vector'(1..d*n => 0);
! 18: res.dg((i-1)*d + j) := 1;
! 19: return res;
! 20: end Create_Term;
! 21:
! 22: function Create_Localized_Term ( n,d,i,j,k : natural ) return Term is
! 23:
! 24: -- DESCRIPTION :
! 25: -- Creates a term in localized variables:
! 26: -- if i >= k, then xij equals either 1 if i=j, or 0 if i/=j.
! 27: -- If i >= k, then the `dg' of the term on return is null.
! 28:
! 29: res : Term;
! 30:
! 31: begin
! 32: if i < k
! 33: then res.cf := Create(1.0);
! 34: res.dg := new Standard_Natural_Vectors.Vector'(1..d*(n-d) => 0);
! 35: res.dg((i-1)*d + j) := 1;
! 36: else if i = j+k-1
! 37: then res.cf := Create(1.0);
! 38: else res.cf := Create(0.0);
! 39: end if;
! 40: return res;
! 41: end if;
! 42: return res;
! 43: end Create_Localized_Term;
! 44:
! 45: function Create_Term ( loc,n,d,i,j : natural ) return Term is
! 46:
! 47: -- DESCRIPTION :
! 48: -- Returns the term that corresponds with xij, with respect to
! 49: -- the localization information in loc, which is 0,1, or 2.
! 50:
! 51: res : Term;
! 52:
! 53: begin
! 54: if loc = 0
! 55: then res.cf := Create(0.0);
! 56: else res.cf := Create(1.0);
! 57: end if;
! 58: res.dg := new Standard_Natural_Vectors.Vector'(1..d*n => 0);
! 59: if loc = 2
! 60: then res.dg((i-1)*d + j) := 1;
! 61: end if;
! 62: return res;
! 63: end Create_Term;
! 64:
! 65: -- TARGET ROUTINES :
! 66:
! 67: function Expand2 ( n,d : natural; b : Bracket ) return Poly is
! 68:
! 69: -- DESCRIPTION :
! 70: -- Does the expansion in the two-dimensional case.
! 71:
! 72: res,subtmp : Poly;
! 73: b11 : Term := Create_Term(n,d,b(1),1);
! 74: b12 : Term := Create_Term(n,d,b(1),2);
! 75: b21 : Term := Create_Term(n,d,b(2),1);
! 76: b22 : Term := Create_Term(n,d,b(2),2);
! 77: d1 : Poly := Create(b11);
! 78: d2 : Poly := Create(b21);
! 79:
! 80: begin -- res := b11*b22 - b21*b12;
! 81: res := b22*d1;
! 82: subtmp := b12*d2;
! 83: Sub(res,subtmp);
! 84: Clear(subtmp);
! 85: Clear(b22); Clear(b12);
! 86: Clear(b11); Clear(b21);
! 87: Clear(d1); Clear(d2);
! 88: return res;
! 89: end Expand2;
! 90:
! 91: function Expand ( n,d : natural; b : Bracket ) return Poly is
! 92:
! 93: res : Poly;
! 94:
! 95: begin
! 96: if b'last = 2
! 97: then res := Expand2(n,d,b);
! 98: else res := Null_Poly;
! 99: declare
! 100: sig : integer := +1;
! 101: b1 : Bracket(1..b'last-1);
! 102: begin
! 103: b1(1..b'last-1) := b(1..b'last-1);
! 104: if b'last mod 2 = 0
! 105: then sig := -1;
! 106: end if;
! 107: for i in b'range loop
! 108: declare
! 109: xt : Term := Create_Term(n,d,b(i),b'last);
! 110: expb1,xtexpb1 : Poly;
! 111: begin
! 112: b1(1..i-1) := b(1..i-1);
! 113: b1(i..b1'last) := b(i+1..b'last);
! 114: expb1 := Expand(n,d,b1);
! 115: xtexpb1 := xt*expb1;
! 116: if sig > 0
! 117: then Add(res,xtexpb1);
! 118: else Sub(res,xtexpb1);
! 119: end if;
! 120: sig := -sig;
! 121: Clear(expb1); Clear(xt); Clear(xtexpb1);
! 122: end;
! 123: end loop;
! 124: end;
! 125: end if;
! 126: return res;
! 127: end Expand;
! 128:
! 129: function Expand ( n,d : natural; b : Bracket_Monomial ) return Poly is
! 130:
! 131: res : Poly := Null_Poly;
! 132: fst : boolean := true;
! 133:
! 134: procedure Expand_Bracket ( bb : in Bracket; continue : out boolean ) is
! 135:
! 136: expbb : Poly;
! 137:
! 138: begin
! 139: if fst
! 140: then res := Expand(n,d,bb);
! 141: fst := false;
! 142: else expbb := Expand(n,d,bb);
! 143: Mul(res,expbb);
! 144: Clear(expbb);
! 145: end if;
! 146: continue := true;
! 147: end Expand_Bracket;
! 148: procedure Expand_Brackets is new Enumerate_Brackets(Expand_Bracket);
! 149:
! 150: begin
! 151: Expand_Brackets(b);
! 152: return res;
! 153: end Expand;
! 154:
! 155: function Expand ( n,d : natural; b : Bracket_Term ) return Poly is
! 156:
! 157: res : Poly := Expand(n,d,b.monom);
! 158:
! 159: begin
! 160: Mul(res,b.coeff);
! 161: return res;
! 162: end Expand;
! 163:
! 164: function Expand ( n,d : natural; b : Bracket_Polynomial ) return Poly is
! 165:
! 166: res : Poly := Null_Poly;
! 167:
! 168: procedure Expand_Term ( t : in Bracket_Term; continue : out boolean ) is
! 169:
! 170: expt : Poly := Expand(n,d,t);
! 171:
! 172: begin
! 173: Add(res,expt);
! 174: Clear(expt);
! 175: continue := true;
! 176: end Expand_Term;
! 177: procedure Expand_Terms is new Enumerate_Terms(Expand_Term);
! 178:
! 179: begin
! 180: Expand_Terms(b);
! 181: return res;
! 182: end Expand;
! 183:
! 184: function Localized_Expand2
! 185: ( n,d : natural; b : Bracket; k : natural ) return Poly is
! 186:
! 187: -- DESCRIPTION :
! 188: -- Computes the localized expansion for two dimensions.
! 189: -- Exploits the fact that b(1) < b(2).
! 190:
! 191: res : Poly;
! 192: b11 : Term := Create_Localized_Term(n,d,b(1),1,k);
! 193: b12 : Term := Create_Localized_Term(n,d,b(1),2,k);
! 194: b21 : Term := Create_Localized_Term(n,d,b(2),1,k);
! 195: b22 : Term := Create_Localized_Term(n,d,b(2),2,k);
! 196:
! 197: begin
! 198: if b(2) < k
! 199: then declare
! 200: subtmp : Poly;
! 201: d1 : Poly := Create(b11);
! 202: d2 : Poly := Create(b21);
! 203: begin -- res := b11*b22 - b21*b12;
! 204: res := b22*d1;
! 205: subtmp := b12*d2;
! 206: Sub(res,subtmp);
! 207: Clear(subtmp);
! 208: Clear(b22); Clear(b12);
! 209: Clear(b11); Clear(b21);
! 210: Clear(d1); Clear(d2);
! 211: end;
! 212: else if b(1) < k
! 213: then b11.cf := b11.cf*b22.cf;
! 214: b12.cf := b12.cf*b21.cf;
! 215: res := Create(b11);
! 216: Sub(res,b12);
! 217: Clear(b11); Clear(b12);
! 218: else declare
! 219: rt : Term;
! 220: dd : constant natural := d*(n-d);
! 221: begin
! 222: rt.cf := b11.cf*b22.cf - b21.cf*b12.cf;
! 223: rt.dg := new Standard_Natural_Vectors.Vector'(1..dd => 0);
! 224: res := Create(rt);
! 225: end;
! 226: end if;
! 227: end if;
! 228: return res;
! 229: end Localized_Expand2;
! 230:
! 231: function Localized_Expand ( n,d : natural; b : Bracket ) return Poly is
! 232:
! 233: res : Poly;
! 234: k : constant natural := n-d+1;
! 235:
! 236: begin
! 237: if b'last = 2
! 238: then res := Localized_Expand2(n,d,b,k);
! 239: else res := Null_Poly;
! 240: declare
! 241: sig : integer := +1;
! 242: b1 : Bracket(1..b'last-1);
! 243: begin
! 244: b1(1..b'last-1) := b(1..b'last-1);
! 245: if b'last mod 2 = 0
! 246: then sig := -1;
! 247: end if;
! 248: for i in b'range loop
! 249: declare
! 250: xt : Term := Create_Localized_Term(n,d,b(i),b'last,k);
! 251: expb1,xtexpb1 : Poly;
! 252: begin
! 253: if xt.cf /= Create(0.0)
! 254: then b1(1..i-1) := b(1..i-1);
! 255: b1(i..b1'last) := b(i+1..b'last);
! 256: expb1 := Localized_Expand(n,d,b1);
! 257: xtexpb1 := xt*expb1;
! 258: if sig > 0
! 259: then Add(res,xtexpb1);
! 260: else Sub(res,xtexpb1);
! 261: end if;
! 262: Clear(expb1); Clear(xt); Clear(xtexpb1);
! 263: end if;
! 264: sig := -sig;
! 265: end;
! 266: end loop;
! 267: end;
! 268: end if;
! 269: return res;
! 270: end Localized_Expand;
! 271:
! 272: function Localization_Map ( n,d : natural ) return Matrix is
! 273:
! 274: res : Matrix(1..n,1..d);
! 275: row,col : natural;
! 276:
! 277: begin
! 278: for i in 1..n loop -- initialization
! 279: for j in 1..d loop
! 280: res(i,j) := -1; -- means empty space
! 281: end loop;
! 282: end loop;
! 283: col := 0;
! 284: for i in 1..n loop -- one free element in every row
! 285: col := col+1;
! 286: if col > d
! 287: then col := 1;
! 288: end if;
! 289: res(i,col) := 2;
! 290: end loop;
! 291: row := 0;
! 292: for j in 1..d loop -- one 1 in every column
! 293: loop
! 294: row := row+1;
! 295: if row > n
! 296: then row := 1;
! 297: end if;
! 298: exit when (res(row,j) = -1); -- empty space found
! 299: end loop;
! 300: res(row,j) := 1;
! 301: end loop;
! 302: for j in 1..d loop -- fill in d-1 zeros in every column
! 303: for i in 1..(d-1) loop
! 304: row := 0;
! 305: loop
! 306: row := row+1;
! 307: if row > n
! 308: then row := 1;
! 309: end if;
! 310: exit when (res(row,j) = -1); -- empty space found
! 311: end loop;
! 312: res(row,j) := 0;
! 313: end loop;
! 314: end loop;
! 315: for i in 1..n loop -- fill rest with free elements
! 316: for j in 1..d loop
! 317: if res(i,j) = -1
! 318: then res(i,j) := 2;
! 319: end if;
! 320: end loop;
! 321: end loop;
! 322: return res;
! 323: end Localization_Map;
! 324:
! 325: function Expand2 ( locmap : Matrix; b : Bracket ) return Poly is
! 326:
! 327: -- DESCRIPTION :
! 328: -- Expands a 2-by-2 minor of the matrix, selecting the rows with
! 329: -- entries in the 2-bracket b, respecting the localization map.
! 330:
! 331: res,subtmp : Poly;
! 332: n : constant natural := locmap'length(1);
! 333: d : constant natural := locmap'length(2);
! 334: b11 : Term := Create_Term(locmap(b(1),1),n,d,b(1),1);
! 335: b12 : Term := Create_Term(locmap(b(1),2),n,d,b(1),2);
! 336: b21 : Term := Create_Term(locmap(b(2),1),n,d,b(2),1);
! 337: b22 : Term := Create_Term(locmap(b(2),2),n,d,b(2),2);
! 338: d1 : Poly := Create(b11);
! 339: d2 : Poly := Create(b21);
! 340:
! 341: begin -- res := b11*b22 - b21*b12;
! 342: res := b22*d1;
! 343: subtmp := b12*d2;
! 344: Sub(res,subtmp);
! 345: Clear(subtmp);
! 346: Clear(b22); Clear(b12);
! 347: Clear(b11); Clear(b21);
! 348: Clear(d1); Clear(d2);
! 349: return res;
! 350: end Expand2;
! 351:
! 352: function Expand ( locmap : Matrix; b : Bracket ) return Poly is
! 353:
! 354: -- DESCRIPTION :
! 355: -- Expands a d-by-d minor of the matrix, selecting the rows with
! 356: -- entries in b, respecting the localization map.
! 357: -- The expansion starts at the last column and proceeds recursively.
! 358:
! 359: res : Poly;
! 360: n : constant natural := locmap'length(1);
! 361: d : constant natural := locmap'length(2);
! 362:
! 363: begin
! 364: if b'last = 2
! 365: then res := Expand2(locmap,b);
! 366: else res := Null_Poly;
! 367: declare
! 368: sig : integer := +1;
! 369: b1 : Bracket(1..b'last-1);
! 370: begin
! 371: b1(1..b'last-1) := b(1..b'last-1);
! 372: if b'last mod 2 = 0
! 373: then sig := -1;
! 374: end if;
! 375: for i in b'range loop
! 376: declare
! 377: xt : Term := Create_Term(locmap(b(i),b'last),n,d,b(i),b'last);
! 378: expb1,xtexpb1 : Poly;
! 379: begin
! 380: if xt.cf /= Create(0.0)
! 381: then b1(1..i-1) := b(1..i-1);
! 382: b1(i..b1'last) := b(i+1..b'last);
! 383: expb1 := Expand(locmap,b1);
! 384: xtexpb1 := xt*expb1;
! 385: if sig > 0
! 386: then Add(res,xtexpb1);
! 387: else Sub(res,xtexpb1);
! 388: end if;
! 389: Clear(expb1); Clear(xtexpb1);
! 390: end if;
! 391: Clear(xt);
! 392: end;
! 393: sig := -sig;
! 394: end loop;
! 395: end;
! 396: end if;
! 397: return res;
! 398: end Expand;
! 399:
! 400: function Reduce_Variables
! 401: ( locmap : Matrix; dg : Standard_Natural_Vectors.Vector )
! 402: return Standard_Natural_Vectors.Vector is
! 403:
! 404: -- DESCRIPTION :
! 405: -- Removes the #variables in the degrees vectors, removing all variables
! 406: -- that correspond to zeros in the localization map.
! 407:
! 408: res : Standard_Natural_Vectors.Vector(dg'range) := dg;
! 409: cnt : natural := res'last;
! 410: d : constant natural := locmap'length(2);
! 411:
! 412: begin
! 413: for i in reverse locmap'range(1) loop
! 414: for j in reverse locmap'range(2) loop
! 415: if locmap(i,j) /= 2
! 416: then cnt := cnt - 1;
! 417: for k in ((i-1)*d+j)..cnt loop
! 418: res(k) := res(k+1);
! 419: end loop;
! 420: end if;
! 421: end loop;
! 422: end loop;
! 423: return res(res'first..cnt);
! 424: end Reduce_Variables;
! 425:
! 426: function Reduce_Variables ( locmap : Matrix; t : Term ) return Term is
! 427:
! 428: -- DESCRIPTION :
! 429: -- Reduces the #variables in the term, removing all variables that
! 430: -- correspond to zeros in the localization map.
! 431:
! 432: res : Term;
! 433:
! 434: begin
! 435: res.cf := t.cf;
! 436: res.dg := new Standard_Natural_Vectors.Vector'
! 437: (Reduce_Variables(locmap,t.dg.all));
! 438: return res;
! 439: end Reduce_Variables;
! 440:
! 441: procedure Reduce_Variables ( locmap : in Matrix; p : in out Poly ) is
! 442:
! 443: rp : Poly := Null_Poly;
! 444:
! 445: procedure Reduce_Term ( t : in Term; continue : out boolean ) is
! 446:
! 447: rt : Term := Reduce_Variables(locmap,t);
! 448:
! 449: begin
! 450: Add(rp,rt);
! 451: continue := true;
! 452: end Reduce_Term;
! 453: procedure Reduce_Terms is new Visiting_Iterator(Reduce_Term);
! 454:
! 455: begin
! 456: Reduce_Terms(p);
! 457: Clear(p);
! 458: p := rp;
! 459: end Reduce_Variables;
! 460:
! 461: end Bracket_Expansions;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to bracket_expansions.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