Annotation of OpenXM_contrib/PHC/Ada/Schubert/ts_straighten.adb, Revision 1.1.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>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ts_straighten.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