Annotation of OpenXM_contrib/PHC/Ada/Schubert/ts_straighten.adb, Revision 1.1
1.1 ! maekawa 1: with text_io,integer_io; use text_io,integer_io;
! 2: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
! 3: with Brackets,Brackets_io; use Brackets,Brackets_io;
! 4: with Bracket_Monomials; use Bracket_Monomials;
! 5: with Bracket_Polynomials; use Bracket_Polynomials;
! 6: with Bracket_Polynomials_io; use Bracket_Polynomials_io;
! 7: with Straightening_Syzygies; use Straightening_Syzygies;
! 8:
! 9: procedure ts_straighten is
! 10:
! 11: -- DESCRIPTION :
! 12: -- Enumerates all brackets in the Pluecker embedding and generates
! 13: -- the van der Waerden syzygies.
! 14:
! 15: function Number_of_Brackets ( n,d : natural ) return natural is
! 16:
! 17: -- DESCIPTION :
! 18: -- Returns the number of brackets of d entries chosen from n numbers.
! 19:
! 20: a,b : natural;
! 21:
! 22: begin
! 23: a := 1;
! 24: for i in d+1..n loop
! 25: a := a*i;
! 26: end loop;
! 27: b := 1;
! 28: for i in 1..n-d loop
! 29: b := b*i;
! 30: end loop;
! 31: return a/b;
! 32: end Number_of_Brackets;
! 33:
! 34: function Number_of_Linear_Equations ( n,d : natural ) return natural is
! 35:
! 36: -- DESCRIPTION :
! 37: -- Returns the number of linear equations you need to determine a
! 38: -- d-plane in C^n completely.
! 39:
! 40: begin
! 41: return (n-d)*d;
! 42: end Number_of_Linear_Equations;
! 43:
! 44: function Number_of_Zero_Brackets ( n,d : natural ) return natural is
! 45:
! 46: -- DESCRIPTION :
! 47: -- Returns the maximal number of brackets that can be set to zero.
! 48:
! 49: begin
! 50: return (Number_of_Brackets(n,d) - Number_of_Linear_Equations(n,d) - 1);
! 51: end Number_of_Zero_Brackets;
! 52:
! 53: function Create ( b1,b2 : Bracket ) return Bracket_Term is
! 54:
! 55: bm : Bracket_Monomial := Create(b1);
! 56: bt : Bracket_Term;
! 57:
! 58: begin
! 59: Multiply(bm,b2);
! 60: bt.coeff := Create(1.0);
! 61: bt.monom := bm;
! 62: return bt;
! 63: end Create;
! 64:
! 65: procedure Enumerate ( n,d,k,start : in natural; accu : in out Bracket;
! 66: cnt : in out natural ) is
! 67:
! 68: -- DESCRIPTION :
! 69: -- Lexicographic enumeration of all brackets.
! 70:
! 71: begin
! 72: if k > d
! 73: then -- put(accu); new_line;
! 74: cnt := cnt + 1;
! 75: else for l in start..n loop
! 76: accu(k) := l;
! 77: Enumerate(n,d,k+1,l+1,accu,cnt);
! 78: end loop;
! 79: end if;
! 80: end Enumerate;
! 81:
! 82: procedure Enumerate2 ( n,d,k,start : in natural; b : in Bracket;
! 83: accu : in out Bracket;
! 84: cntstd,cntnonstd : in out natural;
! 85: nonstd : in out Bracket_Polynomial ) is
! 86:
! 87: -- DESCRIPTION :
! 88: -- Lexicographic enumeration of all brackets, with b < accu and with
! 89: -- a test whether the pair b*accu forms a Standard monomial.
! 90:
! 91: -- ON ENTRY :
! 92: -- n total number of indices to choose from;
! 93: -- d number of indices in the brackets;
! 94: -- k current entry in accu that is to be filled;
! 95: -- start indices will be taken in between start..n;
! 96: -- b previously enumerated bracket;
! 97: -- accu accumulating parameter, filled in upto (k-1)th entry;
! 98: -- cntstd current number of quadratic standard monomials;
! 99: -- cntnonstd current number of quadratic nonstandard monomials;
! 100: -- nonstd current polynomial of quadratic nonstandard monomials.
! 101:
! 102: -- ON RETURN :
! 103: -- accu updated accumulating parameter, accu(k) is filled in;
! 104: -- cnstd updated number of quadratic standard monomials;
! 105: -- cntnonstd updated number of quadratic nonstandard monomials;
! 106: -- nonstd updated polynomial of quadratic nonstandard monomials.
! 107:
! 108: s : natural;
! 109:
! 110: begin
! 111: if k > d
! 112: then if b < accu
! 113: then -- put(b); put(accu);
! 114: s := Brackets.Is_Standard(b,accu);
! 115: if s = 0
! 116: then -- put_line(" is a Standard monomial.");
! 117: cntstd := cntstd + 1;
! 118: else -- put(" is not a Standard monomial with s = ");
! 119: -- put(s,1); new_line;
! 120: cntnonstd := cntnonstd + 1;
! 121: Add(nonstd,Create(b,accu));
! 122: end if;
! 123: end if;
! 124: else for l in start..n loop
! 125: accu(k) := l;
! 126: Enumerate2(n,d,k+1,l+1,b,accu,cntstd,cntnonstd,nonstd);
! 127: end loop;
! 128: end if;
! 129: end Enumerate2;
! 130:
! 131: procedure Enumerate1 ( n,d,k,start : natural; acc1,acc2 : in out Bracket;
! 132: cntstd,cntnonstd : in out natural;
! 133: nonstd : in out Bracket_Polynomial ) is
! 134:
! 135: -- DESCRIPTION :
! 136: -- Lexicographic enumeration of all brackets, with acc1 < acc2 and with
! 137: -- a test whether the pair acc1*acc2 forms a Standard monomial.
! 138: -- Counts the standard and nonstandard monomials and adds every
! 139: -- nonStandard monomial to the polynomial nonstd.
! 140:
! 141: begin
! 142: if k > d
! 143: then Enumerate2(n,d,1,acc1(1),acc1,acc2,cntstd,cntnonstd,nonstd);
! 144: else for l in start..n loop
! 145: acc1(k) := l;
! 146: Enumerate1(n,d,k+1,l+1,acc1,acc2,cntstd,cntnonstd,nonstd);
! 147: end loop;
! 148: end if;
! 149: end Enumerate1;
! 150:
! 151: procedure Enumerate_Pairs ( n,d : in natural;
! 152: nonstd : in out Bracket_Polynomial ) is
! 153:
! 154: -- DESCRIPTION :
! 155: -- Enumerates all pairs (b1,b2), with b1 < b2 and checks whether
! 156: -- the corresponding bracket monomial b1*b2 is standard or not.
! 157:
! 158: b1,b2 : Bracket(1..d);
! 159: cntstd,cntnonstd : natural := 0;
! 160:
! 161: begin
! 162: Enumerate1(n,d,1,1,b1,b2,cntstd,cntnonstd,nonstd);
! 163: put("#Standard bi-monomials : "); put(cntstd,1); new_line;
! 164: put("#nonStandard bi-monomials : "); put(cntnonstd,1); new_line;
! 165: end Enumerate_Pairs;
! 166:
! 167: procedure Enumerate_Brackets ( n,d : in natural ) is
! 168:
! 169: acc : Bracket(1..d);
! 170: cnt : natural := 0;
! 171: nonstd : Bracket_Polynomial;
! 172:
! 173: begin
! 174: Enumerate(n,d,1,1,acc,cnt);
! 175: put("#brackets : "); put(cnt,1); new_line;
! 176: Enumerate_Pairs(n,d,nonstd);
! 177: put_line("The polynomial of all quadratic nonstandard monomials :");
! 178: put(nonstd);
! 179: end Enumerate_Brackets;
! 180:
! 181: procedure Read_Bracket ( b : in out Bracket ) is
! 182:
! 183: d : constant natural := b'last;
! 184: sig : integer;
! 185:
! 186: begin
! 187: put("Give "); put(d,1); put(" row indices : "); get(b,sig);
! 188: put("The sorted bracket : ");
! 189: if sig > 0
! 190: then put("+");
! 191: else put("-");
! 192: end if;
! 193: put(b);
! 194: if Is_Zero(b)
! 195: then put_line(" is zero.");
! 196: else put_line(" is different from zero.");
! 197: end if;
! 198: end Read_Bracket;
! 199:
! 200: procedure Test_Straighten ( n,d : in natural ) is
! 201:
! 202: b,bb : Bracket(1..d);
! 203: s : natural;
! 204: ans : character;
! 205: bp : Bracket_Polynomial;
! 206:
! 207: begin
! 208: loop
! 209: Read_Bracket(b);
! 210: Read_Bracket(bb);
! 211: if Is_Equal(b,bb)
! 212: then put_line("Both brackets are equal.");
! 213: end if;
! 214: if b < bb
! 215: then put(b); put(" < "); put(bb);
! 216: s := Brackets.Is_Standard(b,bb);
! 217: put(" and "); put(b); put(bb);
! 218: bp := Straightening_Syzygy(b,bb);
! 219: else put(b); put(" >= "); put(bb);
! 220: s := Brackets.Is_Standard(bb,b);
! 221: put(" and "); put(bb); put(b);
! 222: bp := Straightening_Syzygy(bb,b);
! 223: end if;
! 224: if s = 0
! 225: then put_line(" is a Standard monomial.");
! 226: else put(" is not a Standard monomial, with s = "); put(s,1); new_line;
! 227: end if;
! 228: put_line("The straightening syzygy : "); put(bp);
! 229: put("Do you want to test more pairs of brackets ? (y/n) "); get(ans);
! 230: exit when ans /= 'y';
! 231: end loop;
! 232: end Test_Straighten;
! 233:
! 234: procedure Test_Laplace ( n,d : natural ) is
! 235:
! 236: p : Bracket_Polynomial := Laplace_Expansion(n,d);
! 237:
! 238: begin
! 239: put("The Laplace expansion of "); put(n,1); put("*"); put(n,1);
! 240: put("-determinant as product of "); put(d,1); put("- and ");
! 241: put(n-d,1); put_line("-blocks : ");
! 242: put(p);
! 243: end Test_Laplace;
! 244:
! 245: function Decompose ( n,d : natural; p : Bracket_Polynomial )
! 246: return Bracket_Polynomial is
! 247:
! 248: -- DESCRIPTION :
! 249: -- Decomposes the ideal of square-free nonStandard monomials.
! 250:
! 251: res : Bracket_Polynomial;
! 252: np : constant natural := Number_of_Monomials(p);
! 253: brackmat : array(1..np,1..2) of Link_to_Bracket;
! 254:
! 255: procedure Initialize is
! 256:
! 257: -- DESCRIPTION :
! 258: -- Initializes the matrix of the brackets that occur in p.
! 259:
! 260: cnt_term : natural := 0;
! 261:
! 262: procedure List_Term ( t : in Bracket_Term; continue : out boolean ) is
! 263:
! 264: cnt_mon : natural := 0;
! 265:
! 266: procedure Store_Bracket ( b : in Bracket; continue : out boolean ) is
! 267: begin
! 268: cnt_mon := cnt_mon + 1;
! 269: brackmat(cnt_term,cnt_mon) := new Bracket'(b);
! 270: continue := true;
! 271: end Store_Bracket;
! 272: procedure Store_Brackets is
! 273: new Bracket_Monomials.Enumerate_Brackets(Store_Bracket);
! 274:
! 275: begin
! 276: cnt_term := cnt_term+1;
! 277: Store_Brackets(t.monom);
! 278: continue := true;
! 279: end List_Term;
! 280: procedure List_Terms is new Enumerate_Terms(List_Term);
! 281:
! 282: begin
! 283: List_Terms(p);
! 284: end Initialize;
! 285:
! 286: procedure Write is
! 287:
! 288: -- DESCRIPTION :
! 289: -- Writes the matrix of brackets.
! 290:
! 291: begin
! 292: for i in 1..np loop
! 293: for j in 1..2 loop
! 294: put(" ");
! 295: put("b("); put(i,1); put(","); put(j,1); put(") :");
! 296: put(brackmat(i,j).all);
! 297: end loop;
! 298: new_line;
! 299: end loop;
! 300: end Write;
! 301:
! 302: function Occured_Yet ( k : natural; b : Bracket ) return boolean is
! 303:
! 304: -- DESCRIPTION :
! 305: -- Returns true if the bracket b occurs in one of the rows <= k.
! 306:
! 307: begin
! 308: for i in 1..k loop
! 309: for j in 1..2 loop
! 310: if Is_Equal(brackmat(i,j).all,b)
! 311: then return true;
! 312: end if;
! 313: end loop;
! 314: end loop;
! 315: return false;
! 316: end Occured_Yet;
! 317:
! 318: procedure Solve is
! 319:
! 320: -- DESCRIPTION :
! 321: -- Solves the initial ideal of quadratic nonStandard monomials.
! 322:
! 323: accu : Bracket_Monomial;
! 324: nzrb : constant natural := Number_of_Zero_Brackets(n,d);
! 325:
! 326: procedure Solve ( k : in natural; acc : in Bracket_Monomial ) is
! 327: begin
! 328: if k > np
! 329: then Add(res,Create(acc));
! 330: else if Divisible(acc,brackmat(k,1).all)
! 331: or else Divisible(acc,brackmat(k,2).all)
! 332: then Solve(k+1,acc);
! 333: else if Number_of_Brackets(acc) < nzrb
! 334: then for i in 1..2 loop
! 335: declare
! 336: newacc : Bracket_Monomial;
! 337: begin
! 338: Copy(acc,newacc);
! 339: Multiply(newacc,brackmat(k,i).all);
! 340: Solve(k+1,newacc);
! 341: end;
! 342: end loop;
! 343: end if;
! 344: end if;
! 345: end if;
! 346: end Solve;
! 347:
! 348: begin
! 349: Solve(1,accu);
! 350: end Solve;
! 351:
! 352: begin
! 353: Initialize;
! 354: -- Write;
! 355: Solve;
! 356: return res;
! 357: end Decompose;
! 358:
! 359: procedure Enumerate_Straightening_Syzygies ( n,d : in natural ) is
! 360:
! 361: -- DESCRIPTION :
! 362: -- Sets up of the straightening syzygies for all nonStandard
! 363: -- quadratic monomials.
! 364:
! 365: nonstd,decom : Bracket_Polynomial;
! 366:
! 367: procedure Write_Syzygy ( s : in Bracket_Polynomial; cont : out boolean ) is
! 368: begin
! 369: put(s); cont := true;
! 370: end Write_Syzygy;
! 371: procedure Write_Syzygies is new Enumerate_Syzygies(Write_Syzygy);
! 372:
! 373: begin
! 374: put("(n,d) = "); put("("); put(n,1); put(","); put(d,1); put(")");
! 375: put(" #Brackets : "); put(Number_of_Brackets(n,d),1);
! 376: put(" #Linear Equations : ");
! 377: put(Number_of_LInear_Equations(n,d),1); new_line;
! 378: put_line("Type of linear equation :");
! 379: put(Laplace_Expansion(n,n-d));
! 380: nonstd := nonStandard_Monomials(n,d);
! 381: put_line("The polynomial with all quadratic nonStandard monomials :");
! 382: put(nonstd);
! 383: put_line("All quadratic straightening syzygies : ");
! 384: Write_Syzygies(nonstd);
! 385: put("#Zero-Brackets : "); put(Number_of_Zero_Brackets(n,d),1); new_line;
! 386: decom := Decompose(n,d,nonstd);
! 387: put_line("Decomposition of the ideal of nonStandard monomials :");
! 388: put(decom);
! 389: put("Number of solutions : "); put(Number_of_Monomials(decom),1);
! 390: end Enumerate_Straightening_Syzygies;
! 391:
! 392: procedure Main is
! 393:
! 394: n,d : natural;
! 395: ans : character;
! 396:
! 397: begin
! 398: new_line;
! 399: put_line("Interactive testing of straightening algorithm.");
! 400: loop
! 401: new_line;
! 402: put_line("Choose one of the following : ");
! 403: put_line(" 0. Exit this program.");
! 404: put_line(" 1. Enumerate all brackets.");
! 405: put_line(" 2. Test the straightening algorithm.");
! 406: put_line(" 3. Test Laplace expansion.");
! 407: put_line(" 4. Enumerate straightening syzygies.");
! 408: put("Make your choice (0,1,2,3,or 4) : "); get(ans);
! 409: exit when ans = '0';
! 410: new_line;
! 411: put("Give the number of entries in bracket : "); get(d);
! 412: put("Give the number of elements to choose from : "); get(n);
! 413: case ans is
! 414: when '1' => Enumerate_Brackets(n,d);
! 415: when '2' => Test_Straighten(n,d);
! 416: when '3' => Test_Laplace(n,d);
! 417: when '4' => Enumerate_Straightening_Syzygies(n,d);
! 418: when others => put_line("Bad answer.");
! 419: end case;
! 420: end loop;
! 421: end Main;
! 422: begin
! 423: Main;
! 424: end ts_straighten;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>