![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to bracket_expansions.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|
Return to bracket_expansions.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Schubert |
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 Bracket_Expansions is
5:
6: -- AUXILIARIES :
7:
8: function Create_Term ( n,d,i,j : natural ) return Term is
9:
10: -- DESCRIPTION :
11: -- Returns the term that corresponds with xij.
12:
13: res : Term;
14:
15: begin
16: res.cf := Create(1.0);
17: res.dg := new Standard_Natural_Vectors.Vector'(1..d*n => 0);
18: res.dg((i-1)*d + j) := 1;
19: return res;
20: end Create_Term;
21:
22: function Create_Localized_Term ( n,d,i,j,k : natural ) return Term is
23:
24: -- DESCRIPTION :
25: -- Creates a term in localized variables:
26: -- if i >= k, then xij equals either 1 if i=j, or 0 if i/=j.
27: -- If i >= k, then the `dg' of the term on return is null.
28:
29: res : Term;
30:
31: begin
32: if i < k
33: then res.cf := Create(1.0);
34: res.dg := new Standard_Natural_Vectors.Vector'(1..d*(n-d) => 0);
35: res.dg((i-1)*d + j) := 1;
36: else if i = j+k-1
37: then res.cf := Create(1.0);
38: else res.cf := Create(0.0);
39: end if;
40: return res;
41: end if;
42: return res;
43: end Create_Localized_Term;
44:
45: function Create_Term ( loc,n,d,i,j : natural ) return Term is
46:
47: -- DESCRIPTION :
48: -- Returns the term that corresponds with xij, with respect to
49: -- the localization information in loc, which is 0,1, or 2.
50:
51: res : Term;
52:
53: begin
54: if loc = 0
55: then res.cf := Create(0.0);
56: else res.cf := Create(1.0);
57: end if;
58: res.dg := new Standard_Natural_Vectors.Vector'(1..d*n => 0);
59: if loc = 2
60: then res.dg((i-1)*d + j) := 1;
61: end if;
62: return res;
63: end Create_Term;
64:
65: -- TARGET ROUTINES :
66:
67: function Expand2 ( n,d : natural; b : Bracket ) return Poly is
68:
69: -- DESCRIPTION :
70: -- Does the expansion in the two-dimensional case.
71:
72: res,subtmp : Poly;
73: b11 : Term := Create_Term(n,d,b(1),1);
74: b12 : Term := Create_Term(n,d,b(1),2);
75: b21 : Term := Create_Term(n,d,b(2),1);
76: b22 : Term := Create_Term(n,d,b(2),2);
77: d1 : Poly := Create(b11);
78: d2 : Poly := Create(b21);
79:
80: begin -- res := b11*b22 - b21*b12;
81: res := b22*d1;
82: subtmp := b12*d2;
83: Sub(res,subtmp);
84: Clear(subtmp);
85: Clear(b22); Clear(b12);
86: Clear(b11); Clear(b21);
87: Clear(d1); Clear(d2);
88: return res;
89: end Expand2;
90:
91: function Expand ( n,d : natural; b : Bracket ) return Poly is
92:
93: res : Poly;
94:
95: begin
96: if b'last = 2
97: then res := Expand2(n,d,b);
98: else res := Null_Poly;
99: declare
100: sig : integer := +1;
101: b1 : Bracket(1..b'last-1);
102: begin
103: b1(1..b'last-1) := b(1..b'last-1);
104: if b'last mod 2 = 0
105: then sig := -1;
106: end if;
107: for i in b'range loop
108: declare
109: xt : Term := Create_Term(n,d,b(i),b'last);
110: expb1,xtexpb1 : Poly;
111: begin
112: b1(1..i-1) := b(1..i-1);
113: b1(i..b1'last) := b(i+1..b'last);
114: expb1 := Expand(n,d,b1);
115: xtexpb1 := xt*expb1;
116: if sig > 0
117: then Add(res,xtexpb1);
118: else Sub(res,xtexpb1);
119: end if;
120: sig := -sig;
121: Clear(expb1); Clear(xt); Clear(xtexpb1);
122: end;
123: end loop;
124: end;
125: end if;
126: return res;
127: end Expand;
128:
129: function Expand ( n,d : natural; b : Bracket_Monomial ) return Poly is
130:
131: res : Poly := Null_Poly;
132: fst : boolean := true;
133:
134: procedure Expand_Bracket ( bb : in Bracket; continue : out boolean ) is
135:
136: expbb : Poly;
137:
138: begin
139: if fst
140: then res := Expand(n,d,bb);
141: fst := false;
142: else expbb := Expand(n,d,bb);
143: Mul(res,expbb);
144: Clear(expbb);
145: end if;
146: continue := true;
147: end Expand_Bracket;
148: procedure Expand_Brackets is new Enumerate_Brackets(Expand_Bracket);
149:
150: begin
151: Expand_Brackets(b);
152: return res;
153: end Expand;
154:
155: function Expand ( n,d : natural; b : Bracket_Term ) return Poly is
156:
157: res : Poly := Expand(n,d,b.monom);
158:
159: begin
160: Mul(res,b.coeff);
161: return res;
162: end Expand;
163:
164: function Expand ( n,d : natural; b : Bracket_Polynomial ) return Poly is
165:
166: res : Poly := Null_Poly;
167:
168: procedure Expand_Term ( t : in Bracket_Term; continue : out boolean ) is
169:
170: expt : Poly := Expand(n,d,t);
171:
172: begin
173: Add(res,expt);
174: Clear(expt);
175: continue := true;
176: end Expand_Term;
177: procedure Expand_Terms is new Enumerate_Terms(Expand_Term);
178:
179: begin
180: Expand_Terms(b);
181: return res;
182: end Expand;
183:
184: function Localized_Expand2
185: ( n,d : natural; b : Bracket; k : natural ) return Poly is
186:
187: -- DESCRIPTION :
188: -- Computes the localized expansion for two dimensions.
189: -- Exploits the fact that b(1) < b(2).
190:
191: res : Poly;
192: b11 : Term := Create_Localized_Term(n,d,b(1),1,k);
193: b12 : Term := Create_Localized_Term(n,d,b(1),2,k);
194: b21 : Term := Create_Localized_Term(n,d,b(2),1,k);
195: b22 : Term := Create_Localized_Term(n,d,b(2),2,k);
196:
197: begin
198: if b(2) < k
199: then declare
200: subtmp : Poly;
201: d1 : Poly := Create(b11);
202: d2 : Poly := Create(b21);
203: begin -- res := b11*b22 - b21*b12;
204: res := b22*d1;
205: subtmp := b12*d2;
206: Sub(res,subtmp);
207: Clear(subtmp);
208: Clear(b22); Clear(b12);
209: Clear(b11); Clear(b21);
210: Clear(d1); Clear(d2);
211: end;
212: else if b(1) < k
213: then b11.cf := b11.cf*b22.cf;
214: b12.cf := b12.cf*b21.cf;
215: res := Create(b11);
216: Sub(res,b12);
217: Clear(b11); Clear(b12);
218: else declare
219: rt : Term;
220: dd : constant natural := d*(n-d);
221: begin
222: rt.cf := b11.cf*b22.cf - b21.cf*b12.cf;
223: rt.dg := new Standard_Natural_Vectors.Vector'(1..dd => 0);
224: res := Create(rt);
225: end;
226: end if;
227: end if;
228: return res;
229: end Localized_Expand2;
230:
231: function Localized_Expand ( n,d : natural; b : Bracket ) return Poly is
232:
233: res : Poly;
234: k : constant natural := n-d+1;
235:
236: begin
237: if b'last = 2
238: then res := Localized_Expand2(n,d,b,k);
239: else res := Null_Poly;
240: declare
241: sig : integer := +1;
242: b1 : Bracket(1..b'last-1);
243: begin
244: b1(1..b'last-1) := b(1..b'last-1);
245: if b'last mod 2 = 0
246: then sig := -1;
247: end if;
248: for i in b'range loop
249: declare
250: xt : Term := Create_Localized_Term(n,d,b(i),b'last,k);
251: expb1,xtexpb1 : Poly;
252: begin
253: if xt.cf /= Create(0.0)
254: then b1(1..i-1) := b(1..i-1);
255: b1(i..b1'last) := b(i+1..b'last);
256: expb1 := Localized_Expand(n,d,b1);
257: xtexpb1 := xt*expb1;
258: if sig > 0
259: then Add(res,xtexpb1);
260: else Sub(res,xtexpb1);
261: end if;
262: Clear(expb1); Clear(xt); Clear(xtexpb1);
263: end if;
264: sig := -sig;
265: end;
266: end loop;
267: end;
268: end if;
269: return res;
270: end Localized_Expand;
271:
272: function Localization_Map ( n,d : natural ) return Matrix is
273:
274: res : Matrix(1..n,1..d);
275: row,col : natural;
276:
277: begin
278: for i in 1..n loop -- initialization
279: for j in 1..d loop
280: res(i,j) := -1; -- means empty space
281: end loop;
282: end loop;
283: col := 0;
284: for i in 1..n loop -- one free element in every row
285: col := col+1;
286: if col > d
287: then col := 1;
288: end if;
289: res(i,col) := 2;
290: end loop;
291: row := 0;
292: for j in 1..d loop -- one 1 in every column
293: loop
294: row := row+1;
295: if row > n
296: then row := 1;
297: end if;
298: exit when (res(row,j) = -1); -- empty space found
299: end loop;
300: res(row,j) := 1;
301: end loop;
302: for j in 1..d loop -- fill in d-1 zeros in every column
303: for i in 1..(d-1) loop
304: row := 0;
305: loop
306: row := row+1;
307: if row > n
308: then row := 1;
309: end if;
310: exit when (res(row,j) = -1); -- empty space found
311: end loop;
312: res(row,j) := 0;
313: end loop;
314: end loop;
315: for i in 1..n loop -- fill rest with free elements
316: for j in 1..d loop
317: if res(i,j) = -1
318: then res(i,j) := 2;
319: end if;
320: end loop;
321: end loop;
322: return res;
323: end Localization_Map;
324:
325: function Expand2 ( locmap : Matrix; b : Bracket ) return Poly is
326:
327: -- DESCRIPTION :
328: -- Expands a 2-by-2 minor of the matrix, selecting the rows with
329: -- entries in the 2-bracket b, respecting the localization map.
330:
331: res,subtmp : Poly;
332: n : constant natural := locmap'length(1);
333: d : constant natural := locmap'length(2);
334: b11 : Term := Create_Term(locmap(b(1),1),n,d,b(1),1);
335: b12 : Term := Create_Term(locmap(b(1),2),n,d,b(1),2);
336: b21 : Term := Create_Term(locmap(b(2),1),n,d,b(2),1);
337: b22 : Term := Create_Term(locmap(b(2),2),n,d,b(2),2);
338: d1 : Poly := Create(b11);
339: d2 : Poly := Create(b21);
340:
341: begin -- res := b11*b22 - b21*b12;
342: res := b22*d1;
343: subtmp := b12*d2;
344: Sub(res,subtmp);
345: Clear(subtmp);
346: Clear(b22); Clear(b12);
347: Clear(b11); Clear(b21);
348: Clear(d1); Clear(d2);
349: return res;
350: end Expand2;
351:
352: function Expand ( locmap : Matrix; b : Bracket ) return Poly is
353:
354: -- DESCRIPTION :
355: -- Expands a d-by-d minor of the matrix, selecting the rows with
356: -- entries in b, respecting the localization map.
357: -- The expansion starts at the last column and proceeds recursively.
358:
359: res : Poly;
360: n : constant natural := locmap'length(1);
361: d : constant natural := locmap'length(2);
362:
363: begin
364: if b'last = 2
365: then res := Expand2(locmap,b);
366: else res := Null_Poly;
367: declare
368: sig : integer := +1;
369: b1 : Bracket(1..b'last-1);
370: begin
371: b1(1..b'last-1) := b(1..b'last-1);
372: if b'last mod 2 = 0
373: then sig := -1;
374: end if;
375: for i in b'range loop
376: declare
377: xt : Term := Create_Term(locmap(b(i),b'last),n,d,b(i),b'last);
378: expb1,xtexpb1 : Poly;
379: begin
380: if xt.cf /= Create(0.0)
381: then b1(1..i-1) := b(1..i-1);
382: b1(i..b1'last) := b(i+1..b'last);
383: expb1 := Expand(locmap,b1);
384: xtexpb1 := xt*expb1;
385: if sig > 0
386: then Add(res,xtexpb1);
387: else Sub(res,xtexpb1);
388: end if;
389: Clear(expb1); Clear(xtexpb1);
390: end if;
391: Clear(xt);
392: end;
393: sig := -sig;
394: end loop;
395: end;
396: end if;
397: return res;
398: end Expand;
399:
400: function Reduce_Variables
401: ( locmap : Matrix; dg : Standard_Natural_Vectors.Vector )
402: return Standard_Natural_Vectors.Vector is
403:
404: -- DESCRIPTION :
405: -- Removes the #variables in the degrees vectors, removing all variables
406: -- that correspond to zeros in the localization map.
407:
408: res : Standard_Natural_Vectors.Vector(dg'range) := dg;
409: cnt : natural := res'last;
410: d : constant natural := locmap'length(2);
411:
412: begin
413: for i in reverse locmap'range(1) loop
414: for j in reverse locmap'range(2) loop
415: if locmap(i,j) /= 2
416: then cnt := cnt - 1;
417: for k in ((i-1)*d+j)..cnt loop
418: res(k) := res(k+1);
419: end loop;
420: end if;
421: end loop;
422: end loop;
423: return res(res'first..cnt);
424: end Reduce_Variables;
425:
426: function Reduce_Variables ( locmap : Matrix; t : Term ) return Term is
427:
428: -- DESCRIPTION :
429: -- Reduces the #variables in the term, removing all variables that
430: -- correspond to zeros in the localization map.
431:
432: res : Term;
433:
434: begin
435: res.cf := t.cf;
436: res.dg := new Standard_Natural_Vectors.Vector'
437: (Reduce_Variables(locmap,t.dg.all));
438: return res;
439: end Reduce_Variables;
440:
441: procedure Reduce_Variables ( locmap : in Matrix; p : in out Poly ) is
442:
443: rp : Poly := Null_Poly;
444:
445: procedure Reduce_Term ( t : in Term; continue : out boolean ) is
446:
447: rt : Term := Reduce_Variables(locmap,t);
448:
449: begin
450: Add(rp,rt);
451: continue := true;
452: end Reduce_Term;
453: procedure Reduce_Terms is new Visiting_Iterator(Reduce_Term);
454:
455: begin
456: Reduce_Terms(p);
457: Clear(p);
458: p := rp;
459: end Reduce_Variables;
460:
461: end Bracket_Expansions;