Annotation of OpenXM_contrib/PHC/Ada/Schubert/ts_expand.adb, Revision 1.1.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