Annotation of OpenXM_contrib/PHC/Ada/Schubert/straightening_syzygies.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 Straightening_Syzygies is
! 5:
! 6: -- AUXILIARIES :
! 7:
! 8: function Create ( v1,v2 : Vector ) return Bracket_Term is
! 9:
! 10: b1 : Bracket(v1'range);
! 11: b2 : Bracket(v2'range);
! 12: sig1,sig2 : integer;
! 13: bm : Bracket_Monomial;
! 14: bt : Bracket_Term;
! 15:
! 16: begin
! 17: Create(v1,b1,sig1);
! 18: if Is_Zero(b1)
! 19: then bt.coeff := Create(0.0);
! 20: else Create(v2,b2,sig2);
! 21: if Is_Zero(b2)
! 22: then bt.coeff := Create(0.0);
! 23: else bm := Create(b1);
! 24: Multiply(bm,b2);
! 25: if sig1*sig2 > 0
! 26: then bt.coeff := Create(1.0);
! 27: else bt.coeff := -Create(1.0);
! 28: end if;
! 29: end if;
! 30: end if;
! 31: bt.monom := bm;
! 32: return bt;
! 33: end Create;
! 34:
! 35: function Sort ( v : Vector ) return Vector is
! 36:
! 37: -- DESCRIPTION :
! 38: -- Returns the sorted vector v, in increasing order.
! 39:
! 40: res : Vector(v'range) := v;
! 41: ind,min : natural;
! 42:
! 43: begin
! 44: for i in v'first..v'last-1 loop
! 45: min := res(i);
! 46: ind := i;
! 47: for j in i+1..res'last loop
! 48: if res(j) < min
! 49: then min := res(j);
! 50: ind := j;
! 51: end if;
! 52: end loop;
! 53: if ind /= i
! 54: then res(ind) := res(i);
! 55: res(i) := min;
! 56: end if;
! 57: end loop;
! 58: return res;
! 59: end Sort;
! 60:
! 61: function Sign ( v1,v2 : Vector ) return integer is
! 62:
! 63: -- DESCRIPTION :
! 64: -- Returns the sign of the permutation (v1,v2).
! 65:
! 66: d : constant natural := v1'length+v2'length;
! 67: b : Bracket(1..d);
! 68: v : Vector(1..d);
! 69: sig : integer;
! 70:
! 71: begin
! 72: v(v1'range) := v1;
! 73: v(v1'last+1..v'last) := v2;
! 74: Create(v,b,sig);
! 75: return sig;
! 76: end Sign;
! 77:
! 78: function Complement ( n : natural; v : Vector ) return Vector is
! 79:
! 80: -- DESCRIPTION :
! 81: -- Returns the complement of the vector w.r.t. the set 1..n.
! 82:
! 83: res : Vector(1..n-v'length);
! 84: cnt : natural := 0;
! 85: found : boolean;
! 86:
! 87: begin
! 88: for i in 1..n loop
! 89: found := false;
! 90: for j in v'range loop
! 91: if v(j) = i
! 92: then found := true;
! 93: exit;
! 94: end if;
! 95: end loop;
! 96: if not found
! 97: then cnt := cnt + 1;
! 98: res(cnt) := i;
! 99: end if;
! 100: end loop;
! 101: return res;
! 102: end Complement;
! 103:
! 104: procedure Enumerate ( start,i,n : in natural; accu : in out Vector;
! 105: s : in natural; b1,b2 : in Bracket;
! 106: frame : in Vector; bp : in out Bracket_Polynomial ) is
! 107:
! 108: -- DESCRIPTION :
! 109: -- Enumerates all subsets of 1..n, of size accu'length, starting to
! 110: -- fill up accu(i) with entries in start..n.
! 111:
! 112: v1,v2 : Vector(1..b1'last);
! 113:
! 114: begin
! 115: if i > accu'last
! 116: then declare
! 117: comp : constant Vector := Complement(n,accu);
! 118: begin
! 119: v1(b1'first..s-1)
! 120: := Standard_Natural_Vectors.Vector(b1(b1'first..s-1));
! 121: for j in accu'range loop
! 122: v1(s+j-1) := frame(accu(j));
! 123: end loop;
! 124: v2(s+1..b2'last)
! 125: := Standard_Natural_Vectors.Vector(b2(s+1..b2'last));
! 126: for j in comp'range loop
! 127: v2(j) := frame(comp(j));
! 128: end loop;
! 129: declare
! 130: bt : Bracket_Term := Create(v1,v2);
! 131: begin
! 132: if bt.coeff /= Create(0.0)
! 133: then if Sign(accu,comp) > 0
! 134: then Frontal_Add(bp,bt);
! 135: else Frontal_Min(bp,bt);
! 136: end if;
! 137: end if;
! 138: Clear(bt);
! 139: end;
! 140: end;
! 141: else for l in start..n loop
! 142: accu(i) := l;
! 143: Enumerate(l+1,i+1,n,accu,s,b1,b2,frame,bp);
! 144: end loop;
! 145: end if;
! 146: end Enumerate;
! 147:
! 148: function Create ( b1,b2 : Bracket ) return Bracket_Term is
! 149:
! 150: bm : Bracket_Monomial := Create(b1);
! 151: bt : Bracket_Term;
! 152:
! 153: begin
! 154: Multiply(bm,b2);
! 155: bt.coeff := Create(1.0);
! 156: bt.monom := bm;
! 157: return bt;
! 158: end Create;
! 159:
! 160: procedure Enumerate2 ( n,d,k,start : in natural; b : in Bracket;
! 161: accu : in out Bracket;
! 162: cntstd,cntnonstd : in out natural;
! 163: nonstd : in out Bracket_Polynomial ) is
! 164:
! 165: -- DESCRIPTION :
! 166: -- Lexicographic enumeration of all brackets, with b < accu and with
! 167: -- a test whether the pair b*accu forms a Standard monomial.
! 168:
! 169: -- ON ENTRY :
! 170: -- n total number of indices to choose from;
! 171: -- d number of indices in the brackets;
! 172: -- k current entry in accu that is to be filled;
! 173: -- start indices will be taken in between start..n;
! 174: -- b previously enumerated bracket;
! 175: -- accu accumulating parameter, filled in upto (k-1)th entry;
! 176: -- cntstd current number of quadratic standard monomials;
! 177: -- cntnonstd current number of quadratic nonstandard monomials;
! 178: -- nonstd current polynomial of quadratic nonstandard monomials.
! 179:
! 180: -- ON RETURN :
! 181: -- accu updated accumulating parameter, accu(k) is filled in;
! 182: -- cnstd updated number of quadratic standard monomials;
! 183: -- cntnonstd updated number of quadratic nonstandard monomials;
! 184: -- nonstd updated polynomial of quadratic nonstandard monomials.
! 185:
! 186: s : natural;
! 187:
! 188: begin
! 189: if k > d
! 190: then if b < accu
! 191: then s := Brackets.Is_Standard(b,accu);
! 192: if s = 0
! 193: then cntstd := cntstd + 1;
! 194: else cntnonstd := cntnonstd + 1;
! 195: Add(nonstd,Create(b,accu));
! 196: end if;
! 197: end if;
! 198: else for l in start..n loop
! 199: accu(k) := l;
! 200: Enumerate2(n,d,k+1,l+1,b,accu,cntstd,cntnonstd,nonstd);
! 201: end loop;
! 202: end if;
! 203: end Enumerate2;
! 204:
! 205: procedure Enumerate1 ( n,d,k,start : natural; acc1,acc2 : in out Bracket;
! 206: cntstd,cntnonstd : in out natural;
! 207: nonstd : in out Bracket_Polynomial ) is
! 208:
! 209: -- DESCRIPTION :
! 210: -- Lexicographic enumeration of all brackets, with acc1 < acc2 and with
! 211: -- a test whether the pair acc1*acc2 forms a Standard monomial.
! 212: -- Counts the standard and nonstandard monomials and adds every
! 213: -- nonStandard monomial to the polynomial nonstd.
! 214:
! 215: begin
! 216: if k > d
! 217: then Enumerate2(n,d,1,acc1(1),acc1,acc2,cntstd,cntnonstd,nonstd);
! 218: else for l in start..n loop
! 219: acc1(k) := l;
! 220: Enumerate1(n,d,k+1,l+1,acc1,acc2,cntstd,cntnonstd,nonstd);
! 221: end loop;
! 222: end if;
! 223: end Enumerate1;
! 224:
! 225: procedure Enumerate3 ( n,d,k,start : in natural; accu : in out Vector;
! 226: bp : in out Bracket_Polynomial ) is
! 227:
! 228: -- DESCRIPTION :
! 229: -- Lexicographic enumerations of all vectors with d-entries, out
! 230: -- of a set of n natural numbers.
! 231:
! 232: -- ON ENTRY :
! 233: -- n total number of indices to choose from;
! 234: -- d number of indices in the brackets;
! 235: -- k current entry in accu that is to be filled;
! 236: -- start indices will be taken in between start..n;
! 237: -- accu accumulating parameter, filled in upto (k-1)th entry;
! 238:
! 239: -- ON RETURN :
! 240: -- accu filled in upto the kth entry.
! 241:
! 242: begin
! 243: if k > d
! 244: then declare
! 245: comp : constant Vector := Complement(n,accu);
! 246: acc0 : Vector(1..d+1);
! 247: t : Bracket_Term;
! 248: begin
! 249: -- put(" accu : "); put(accu); put(" complement : "); put(comp);
! 250: acc0(1) := 0;
! 251: acc0(2..d+1) := accu(1..d);
! 252: -- put(" acc0 : "); put(acc0); new_line;
! 253: t := Create(acc0,comp);
! 254: if Sign(accu,comp) > 0
! 255: then Frontal_Add(bp,t);
! 256: else Frontal_Min(bp,t);
! 257: end if;
! 258: Clear(t);
! 259: end;
! 260: else for l in start..n-d+k loop
! 261: accu(k) := l;
! 262: Enumerate3(n,d,k+1,l+1,accu,bp);
! 263: end loop;
! 264: end if;
! 265: end Enumerate3;
! 266:
! 267: -- TARGET OPERATIONS :
! 268:
! 269: function Laplace_Expansion ( n,d : natural ) return Bracket_Polynomial is
! 270:
! 271: res : Bracket_Polynomial;
! 272: acc : Vector(1..d);
! 273:
! 274: begin
! 275: Enumerate3(n,d,1,1,acc,res);
! 276: return res;
! 277: end Laplace_Expansion;
! 278:
! 279: function Straightening_Syzygy ( b1,b2 : Bracket )
! 280: return Bracket_Polynomial is
! 281:
! 282: s : constant natural := Is_Standard(b1,b2);
! 283: bm : Bracket_Monomial;
! 284: bp : Bracket_Polynomial;
! 285:
! 286: begin
! 287: if s = 0
! 288: then bm := Create(b1);
! 289: Multiply(bm,b2);
! 290: bp := Create(bm);
! 291: else declare
! 292: d1 : constant natural := b1'last+1;
! 293: frame : Vector(1..d1);
! 294: accu : Vector(1..d1-s);
! 295: begin
! 296: for i in s..b1'last loop
! 297: frame(i-s+1) := b1(i);
! 298: end loop;
! 299: for i in b2'first..s loop
! 300: frame(d1-s+i) := b2(i);
! 301: end loop;
! 302: frame := Sort(frame);
! 303: Enumerate(1,1,d1,accu,s,b1,b2,frame,bp);
! 304: end;
! 305: end if;
! 306: return bp;
! 307: end Straightening_Syzygy;
! 308:
! 309: function Straightening_Syzygy ( b : Bracket_Monomial )
! 310: return Bracket_Polynomial is
! 311:
! 312: b1,b2 : Link_to_Bracket;
! 313: bp : Bracket_Polynomial;
! 314:
! 315: procedure Get_Bracket ( bb : in Bracket; continue : out boolean ) is
! 316: begin
! 317: if b1 = null
! 318: then b1 := new Bracket'(bb);
! 319: else b2 := new Bracket'(bb);
! 320: end if;
! 321: continue := true;
! 322: end Get_Bracket;
! 323: procedure Get_Brackets is new Enumerate_Brackets(Get_Bracket);
! 324:
! 325: begin
! 326: Get_Brackets(b);
! 327: bp := Straightening_Syzygy(b1.all,b2.all);
! 328: Clear(b1); Clear(b2);
! 329: return bp;
! 330: end Straightening_Syzygy;
! 331:
! 332: function nonStandard_Monomials ( n,d : natural ) return Bracket_Polynomial is
! 333:
! 334: nonstd : Bracket_Polynomial;
! 335: b1,b2 : Bracket(1..d);
! 336: cntstd,cntnonstd : natural := 0;
! 337:
! 338: begin
! 339: Enumerate1(n,d,1,1,b1,b2,cntstd,cntnonstd,nonstd);
! 340: return nonstd;
! 341: end nonStandard_Monomials;
! 342:
! 343: procedure Enumerate_Syzygies ( p : in Bracket_Polynomial ) is
! 344:
! 345: procedure Process_Syzygy ( t : in Bracket_Term; continue : out boolean ) is
! 346: begin
! 347: Process(Straightening_Syzygy(t.monom),continue);
! 348: end Process_Syzygy;
! 349: procedure Process_Syzygies is new Enumerate_Terms(Process_Syzygy);
! 350:
! 351: begin
! 352: Process_Syzygies(p);
! 353: end Enumerate_Syzygies;
! 354:
! 355: function Straighten ( b1,b2 : Bracket ) return Bracket_Polynomial is
! 356:
! 357: bp : Bracket_Polynomial := Straightening_Syzygy(b1,b2);
! 358:
! 359: begin
! 360: return bp;
! 361: end Straighten;
! 362:
! 363: function Straighten ( b : Bracket_Monomial ) return Bracket_Polynomial is
! 364:
! 365: bp : Bracket_Polynomial;
! 366:
! 367: begin
! 368: return bp;
! 369: end Straighten;
! 370:
! 371: function Straighten ( b : Bracket_Term ) return Bracket_Polynomial is
! 372:
! 373: bp : Bracket_Polynomial;
! 374:
! 375: begin
! 376: return bp;
! 377: end Straighten;
! 378:
! 379: function Straighten ( b : Bracket_Polynomial ) return Bracket_Polynomial is
! 380:
! 381: bp : Bracket_Polynomial;
! 382:
! 383: begin
! 384: return bp;
! 385: end Straighten;
! 386:
! 387: end Straightening_Syzygies;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to straightening_syzygies.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