Annotation of OpenXM_contrib/PHC/Ada/Schubert/ts_expand.adb, Revision 1.1
1.1 ! maekawa 1: with text_io,integer_io; use text_io,integer_io;
! 2: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
! 3: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
! 4: with Standard_Natural_Vectors; use Standard_Natural_Vectors;
! 5: with Standard_Natural_Vectors_io; use Standard_Natural_Vectors_io;
! 6: with Standard_Complex_Vectors; use Standard_Complex_Vectors;
! 7: with Symbol_Table; use Symbol_Table;
! 8: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
! 9: with Standard_Complex_Polynomials_io; use Standard_Complex_Polynomials_io;
! 10: with Standard_Complex_Poly_Functions; use Standard_Complex_Poly_Functions;
! 11: with Matrix_Indeterminates;
! 12: with Brackets,Brackets_io; use Brackets,Brackets_io;
! 13: with Bracket_Monomials; use Bracket_Monomials;
! 14: with Bracket_Polynomials; use Bracket_Polynomials;
! 15: with Bracket_Polynomials_io; use Bracket_Polynomials_io;
! 16: with Straightening_Syzygies; use Straightening_Syzygies;
! 17: with Bracket_Expansions; use Bracket_Expansions;
! 18:
! 19: procedure ts_expand is
! 20:
! 21: -- DESCRIPTION :
! 22: -- Test on the implementation of bracket expansion.
! 23: --- This procedure pre-dates the package SAGBI_Homotopies.
! 24:
! 25: function Number_of_Brackets ( n,d : natural ) return natural is
! 26:
! 27: -- DESCIPTION :
! 28: -- Returns the number of brackets of d entries chosen from n numbers.
! 29:
! 30: a,b : natural;
! 31:
! 32: begin
! 33: a := 1;
! 34: for i in d+1..n loop
! 35: a := a*i;
! 36: end loop;
! 37: b := 1;
! 38: for i in 1..n-d loop
! 39: b := b*i;
! 40: end loop;
! 41: return a/b;
! 42: end Number_of_Brackets;
! 43:
! 44: procedure Write_Laplace_Expansion
! 45: ( n,d : in natural; p : in Bracket_Polynomial ) is
! 46:
! 47: -- DESCRIPTION :
! 48: -- Writes the Laplace expansion in expanded form, i.e.: with xij's.
! 49:
! 50: procedure Write_Term ( t : in Bracket_Term; continue : out boolean ) is
! 51:
! 52: first : boolean;
! 53:
! 54: procedure Write_Bracket ( b : in Bracket; cont : out boolean ) is
! 55: begin
! 56: if first
! 57: then if REAL_PART(t.coeff) > 0.0
! 58: then put("+");
! 59: else put("-");
! 60: end if;
! 61: put(b); put("*(");
! 62: first := false;
! 63: else put(Expand(n,d,b));
! 64: put(")");
! 65: new_line;
! 66: end if;
! 67: cont := true;
! 68: end Write_Bracket;
! 69: procedure Write_Brackets is new Enumerate_Brackets(Write_Bracket);
! 70:
! 71: begin
! 72: first := true;
! 73: Write_Brackets(t.monom);
! 74: continue := true;
! 75: end Write_Term;
! 76: procedure Write_Terms is new Enumerate_Terms(Write_Term);
! 77:
! 78: begin
! 79: Write_Terms(p);
! 80: end Write_Laplace_Expansion;
! 81:
! 82: function Localize ( p : Poly; n,d : natural ) return Poly is
! 83:
! 84: -- DESCRIPTION :
! 85: -- The last (n-d)*d variables are replaced by ones or zeros.
! 86:
! 87: res,tmp : Poly;
! 88: last : natural := d*n;
! 89:
! 90: begin
! 91: Copy(p,res);
! 92: for i in 1..d loop
! 93: for j in 1..d loop
! 94: if i = j
! 95: then tmp := Eval(res,Create(1.0),last);
! 96: else tmp := Eval(res,Create(0.0),last);
! 97: end if;
! 98: Copy(tmp,res); Clear(tmp);
! 99: last := last-1;
! 100: end loop;
! 101: end loop;
! 102: return res;
! 103: end Localize;
! 104:
! 105: procedure Write_Localized_Laplace_Expansion
! 106: ( p : in Bracket_Polynomial; n,d : in natural ) is
! 107:
! 108: -- DESCRIPTION :
! 109: -- Writes the Laplace expansion in expanded form, i.e.: with xij's
! 110: -- and with the last d*d variables set to zeros or ones.
! 111:
! 112: procedure Write_Term ( t : in Bracket_Term; continue : out boolean ) is
! 113:
! 114: first : boolean;
! 115:
! 116: procedure Write_Bracket ( b : in Bracket; cont : out boolean ) is
! 117: begin
! 118: if first
! 119: then if REAL_PART(t.coeff) > 0.0
! 120: then put("+");
! 121: else put("-");
! 122: end if;
! 123: put(b); put("*(");
! 124: first := false;
! 125: else put(Localize(Expand(n,d,b),n,d));
! 126: put(")");
! 127: new_line;
! 128: end if;
! 129: cont := true;
! 130: end Write_Bracket;
! 131: procedure Write_Brackets is new Enumerate_Brackets(Write_Bracket);
! 132:
! 133: begin
! 134: first := true;
! 135: Write_Brackets(t.monom);
! 136: continue := true;
! 137: end Write_Term;
! 138: procedure Write_Terms is new Enumerate_Terms(Write_Term);
! 139:
! 140: begin
! 141: Write_Terms(p);
! 142: end Write_Localized_Laplace_Expansion;
! 143:
! 144: function Coordinatize ( b : Bracket ) return natural is
! 145:
! 146: -- DESCRIPTION :
! 147: -- Returns the decimal expansion of the bracket.
! 148:
! 149: res : natural := 0;
! 150:
! 151: begin
! 152: for i in b'range loop
! 153: res := res*10 + b(i);
! 154: end loop;
! 155: return res;
! 156: end Coordinatize;
! 157:
! 158: function Coordinatize ( p : Bracket_Polynomial ) return Bracket_Polynomial is
! 159:
! 160: -- DESCRIPTION :
! 161: -- Replaces the first bracket in every monomial by the decimal expansion.
! 162:
! 163: res : Bracket_Polynomial;
! 164:
! 165: procedure Coordinatize_Term ( t : in Bracket_Term; cont1 : out boolean ) is
! 166:
! 167: first,second : boolean;
! 168: bm : Bracket_Monomial;
! 169: bt : Bracket_Term;
! 170:
! 171: procedure Coordinatize_Bracket ( b : in Bracket; cont2 : out boolean ) is
! 172: begin
! 173: if first
! 174: then bt.coeff := Create(double_float(Coordinatize(b)));
! 175: first := false;
! 176: second := true;
! 177: elsif second
! 178: then bm := Create(b);
! 179: else Multiply(bm,b);
! 180: end if;
! 181: cont2 := true;
! 182: end Coordinatize_Bracket;
! 183: procedure Coordinatize_Brackets is
! 184: new Enumerate_Brackets(Coordinatize_Bracket);
! 185:
! 186: begin
! 187: first := true; second := false;
! 188: Coordinatize_Brackets(t.monom);
! 189: bt.monom := bm;
! 190: if REAL_PART(t.coeff) < 0.0
! 191: then Min(res,bt);
! 192: else Add(res,bt);
! 193: end if;
! 194: cont1 := true;
! 195: end Coordinatize_Term;
! 196: procedure Coordinatize_Terms is new Enumerate_Terms(Coordinatize_Term);
! 197:
! 198: begin
! 199: Coordinatize_Terms(p);
! 200: return res;
! 201: end Coordinatize;
! 202:
! 203: function Weight ( e : Standard_Natural_Vectors.Vector; n,d : natural )
! 204: return natural is
! 205:
! 206: -- DESCRIPTION :
! 207: -- Returns the weight of the exponent vector.
! 208:
! 209: res : natural := 0;
! 210: jmp,ind : natural;
! 211:
! 212: begin
! 213: for j in 1..d loop
! 214: jmp := d-j;
! 215: for i in 1..n-d loop
! 216: ind := (i-1)*d + j;
! 217: if e(ind) > 0
! 218: then res := res + e(ind)*(i-1)*jmp;
! 219: end if;
! 220: end loop;
! 221: end loop;
! 222: return res;
! 223: end Weight;
! 224:
! 225: procedure Divide ( p : in out Poly; w : in natural ) is
! 226:
! 227: -- DESCRIPTION :
! 228: -- Divides the polynomial by t^w.
! 229:
! 230: procedure Divide_Term ( t : in out Term; continue : out boolean ) is
! 231: begin
! 232: t.dg(t.dg'last) := t.dg(t.dg'last) - w;
! 233: continue := true;
! 234: end Divide_Term;
! 235: procedure Divide_Terms is new Changing_Iterator(Divide_Term);
! 236:
! 237: begin
! 238: Divide_Terms(p);
! 239: end Divide;
! 240:
! 241: function Lift ( p : Poly; n,d : natural ) return Poly is
! 242:
! 243: -- DESCRIPTION :
! 244: -- Returns the lifted polynomial, which should be an expanded bracket.
! 245:
! 246: res : Poly := Null_Poly;
! 247: minwei : natural := 10000;
! 248:
! 249: procedure Lift_Term ( t : in Term; continue : out boolean ) is
! 250:
! 251: tt : Term;
! 252: wei : natural;
! 253:
! 254: begin
! 255: tt.cf := t.cf;
! 256: tt.dg := new Standard_Natural_Vectors.Vector(1..n*d+1);
! 257: for i in t.dg'range loop
! 258: tt.dg(i) := t.dg(i);
! 259: end loop;
! 260: tt.dg(t.dg'last+1..tt.dg'last) := (t.dg'last+1..tt.dg'last => 0);
! 261: wei := Weight(tt.dg.all,n,d);
! 262: tt.dg(tt.dg'last) := wei;
! 263: Add(res,tt);
! 264: Clear(tt.dg);
! 265: if wei < minwei
! 266: then minwei := wei;
! 267: end if;
! 268: continue := true;
! 269: end Lift_Term;
! 270: procedure Lift_Terms is new Visiting_Iterator(Lift_Term);
! 271:
! 272: begin
! 273: Lift_Terms(p);
! 274: if minwei /= 0
! 275: then Divide(res,minwei);
! 276: end if;
! 277: return res;
! 278: end Lift;
! 279:
! 280: procedure Write_Lifted_Localized_Laplace_Expansion
! 281: ( p : in Bracket_Polynomial; n,d : in natural; l : out Poly ) is
! 282:
! 283: -- DESCRIPTION :
! 284: -- Writes the Laplace expansion in expanded form, i.e.: with xij's
! 285: -- and with the last d*d variables set to zeros or ones.
! 286: -- Returns the lifted coordinatized polynomial.
! 287:
! 288: res : Poly := Null_Poly;
! 289:
! 290: procedure Write_Term ( t : in Bracket_Term; continue : out boolean ) is
! 291:
! 292: first : boolean;
! 293: cf : integer;
! 294:
! 295: procedure Write_Bracket ( b : in Bracket; cont : out boolean ) is
! 296:
! 297: lp : Poly;
! 298:
! 299: begin
! 300: if first
! 301: then cf := Coordinatize(b);
! 302: if REAL_PART(t.coeff) > 0.0
! 303: then put("+");
! 304: else put("-"); cf := -cf;
! 305: end if;
! 306: put(b); put("*(");
! 307: first := false;
! 308: else lp := Lift(Localize(Expand(n,d,b),n,d),n,d);
! 309: put(lp);
! 310: put(")");
! 311: new_line;
! 312: Mul(lp,Create(double_float(cf)));
! 313: Add(res,lp);
! 314: Clear(lp);
! 315: end if;
! 316: cont := true;
! 317: end Write_Bracket;
! 318: procedure Write_Brackets is new Enumerate_Brackets(Write_Bracket);
! 319:
! 320: begin
! 321: first := true;
! 322: Write_Brackets(t.monom);
! 323: continue := true;
! 324: end Write_Term;
! 325: procedure Write_Terms is new Enumerate_Terms(Write_Term);
! 326:
! 327: begin
! 328: Write_Terms(p);
! 329: l := res;
! 330: end Write_Lifted_Localized_Laplace_Expansion;
! 331:
! 332: procedure Expand_Laplace ( n,d : natural ) is
! 333:
! 334: -- DESCRIPTION :
! 335: -- Writes the expanded bracket polynomial.
! 336:
! 337: p : Bracket_Polynomial := Laplace_Expansion(n,n-d);
! 338: q : Bracket_Polynomial; -- := Coordinatize(p);
! 339: lq : Poly; -- := Expand(q);
! 340: ld : Poly; -- := Localize(lq,n,d);
! 341: lt,l0 : Poly;
! 342: tsb : Symbol_Table.Symbol;
! 343:
! 344: begin
! 345: put("The Laplace expansion of "); put(n,1); put("*"); put(n,1);
! 346: put("-determinant as product of "); put(d,1); put("- and ");
! 347: put(n-d,1); put_line("-blocks : ");
! 348: put(p);
! 349: return;
! 350: -- put_line("The coordinatized Laplace expansion : ");
! 351: -- put(q);
! 352: put_line("Expanded in terms of xij's : ");
! 353: Write_Laplace_Expansion(n,d,p);
! 354: -- put_line("The coordinatized Laplace expansion in terms of xij's : ");
! 355: -- put(lq); new_line;
! 356: put_line("The localized version : ");
! 357: Write_Localized_Laplace_Expansion(p,n,d);
! 358: -- put(ld); new_line;
! 359: Symbol_Table.Enlarge(1);
! 360: tsb(1) := 't';
! 361: for i in 2..tsb'last loop
! 362: tsb(i) := ' ';
! 363: end loop;
! 364: Symbol_Table.Add(tsb);
! 365: -- lt := Lift(ld,n,d);
! 366: put_line("The lifted localized version :");
! 367: Write_Lifted_Localized_Laplace_Expansion(p,n,d,lt);
! 368: put_line("The coordinatized lifted localized polynomial : ");
! 369: put(lt); new_line;
! 370: put_line("The polynomial in the start system : ");
! 371: l0 := Eval(lt,Create(0.0),n*d+1);
! 372: put(l0); new_line;
! 373: end Expand_Laplace;
! 374:
! 375: procedure Memory_Consumption ( n,d : natural ) is
! 376:
! 377: nb : natural;
! 378: bp : Bracket_Polynomial;
! 379:
! 380: begin
! 381: put("Give number of expansions : "); get(nb);
! 382: for i in 1..nb loop
! 383: for j in 1..n loop
! 384: bp := Laplace_Expansion(n,j);
! 385: Clear(bp);
! 386: end loop;
! 387: end loop;
! 388: end Memory_Consumption;
! 389:
! 390: procedure Expand_Brackets ( n,d : in natural ) is
! 391:
! 392: -- DESCRIPTION :
! 393: -- Reads a bracket and writes the bracket expansion.
! 394:
! 395: ans : character;
! 396: b : Bracket(1..d);
! 397:
! 398: begin
! 399: loop
! 400: put("Give "); put(d,1); put(" entries for the bracket : ");
! 401: get(b);
! 402: put("The expansion of the bracket "); put(b); put_line(" :");
! 403: put(Expand(n,d,b)); new_line;
! 404: put("Do you want to expand other brackets ? (y/n) "); get(ans);
! 405: exit when ans /= 'y';
! 406: end loop;
! 407: end Expand_Brackets;
! 408:
! 409: procedure Localized_Expand_Brackets ( n,d : in natural ) is
! 410:
! 411: -- DESCRIPTION :
! 412: -- Reads a bracket and writes the bracket expansion.
! 413:
! 414: ans : character;
! 415: b : Bracket(1..d);
! 416:
! 417: begin
! 418: loop
! 419: put("Give "); put(d,1); put(" entries for the bracket : "); get(b);
! 420: put("The expansion of the bracket "); put(b); put_line(" :");
! 421: put(Localized_Expand(n,d,b)); new_line;
! 422: put("Do you want to expand other brackets ? (y/n) "); get(ans);
! 423: exit when ans /= 'y';
! 424: end loop;
! 425: end Localized_Expand_Brackets;
! 426:
! 427: procedure Main is
! 428:
! 429: n,d : natural;
! 430: ans : character;
! 431:
! 432: begin
! 433: new_line;
! 434: put_line("Testing bracket expansions");
! 435: loop
! 436: new_line;
! 437: put("Give number of elements to choose from : "); get(n);
! 438: put("Give the number of entries in bracket : "); get(d);
! 439: Matrix_Indeterminates.Initialize_Symbols(n,d);
! 440: put_line("Choose one of the following :");
! 441: put_line(" 1. Expand single brackets.");
! 442: put_line(" 2. Expand single brackets in local coordinates.");
! 443: put_line(" 3. Apply Laplace expansion.");
! 444: put_line(" 4. Test memory consumption.");
! 445: put("Make your choice (1,2,3, or 4) : "); get(ans);
! 446: put("(n,d) = ("); put(n,1); put(","); put(d,1); put(")");
! 447: put(" #brackets : "); put(Number_of_Brackets(n,d),1);
! 448: put(" #equations : "); put((n-d)*d,1); new_line;
! 449: case ans is
! 450: when '1' => Expand_Brackets(n,d);
! 451: when '2' => Localized_Expand_Brackets(n,d);
! 452: when '3' => Expand_Laplace(n,d);
! 453: when '4' => Memory_Consumption(n,d);
! 454: when others => put_line("option not available");
! 455: end case;
! 456: Matrix_Indeterminates.Clear_Symbols;
! 457: put("Do you want more tests for other n and d ? (y/n) "); get(ans);
! 458: exit when ans /= 'y';
! 459: end loop;
! 460: end Main;
! 461:
! 462: begin
! 463: Main;
! 464: end ts_expand;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ts_expand.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