Annotation of OpenXM_contrib/PHC/Ada/Schubert/symbolic_minor_equations.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
2: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
3: with Standard_Natural_Vectors;
4: with Matrix_Indeterminates;
5: with Bracket_Polynomials; use Bracket_Polynomials;
6: with Straightening_Syzygies; use Straightening_Syzygies;
7:
8: package body Symbolic_Minor_Equations is
9:
10: -- AUXILIARIES TO Minor_Equations :
11:
12: function Substitute ( b,minor : Bracket ) return Bracket is
13:
14: -- DESCRIPTION :
15: -- The ith entry in the bracket b is replaced by minor(b(i)).
16:
17: res : Bracket(b'range);
18: start : integer;
19:
20: begin
21: if b(b'first) = 0
22: then start := b'first+1; -- bracket is a coefficient bracket
23: res(res'first) := 0; -- result is also a coefficient bracket
24: else start := b'first; -- second bracket in expansion
25: end if;
26: for i in start..b'last loop
27: res(i) := minor(b(i));
28: end loop;
29: return res;
30: end Substitute;
31:
32: function Substitute ( bm : Bracket_Monomial; minor : Bracket )
33: return Bracket_Monomial is
34:
35: -- DESCRIPTION :
36: -- Substitutes the entries in the brackets according to the minor.
37:
38: -- REQUIRED : bm is a quadratic monomial.
39:
40: res : Bracket_Monomial;
41: lb : Link_to_Bracket;
42: first : boolean := true;
43:
44: procedure Substitute_Bracket ( b : in Bracket; continue : out boolean ) is
45:
46: sb : Bracket(b'range) := Substitute(b,minor);
47:
48: begin
49: if first
50: then lb := new Bracket'(sb); -- save to preserve order
51: first := false;
52: else res := sb*lb.all;
53: end if;
54: continue := true;
55: end Substitute_Bracket;
56: procedure Substitute_Brackets is
57: new Enumerate_Brackets(Substitute_Bracket);
58:
59: begin
60: Substitute_Brackets(bm);
61: Clear(lb);
62: return res;
63: end Substitute;
64:
65: function Substitute ( bt : Bracket_Term; minor : Bracket )
66: return Bracket_Term is
67:
68: -- DESCRIPTION :
69: -- The ith entry in every bracket of the term according to the minor.
70:
71: res : Bracket_Term;
72:
73: begin
74: res.coeff := bt.coeff;
75: res.monom := Substitute(bt.monom,minor);
76: return res;
77: end Substitute;
78:
79: function Substitute ( bp : Bracket_Polynomial; minor : Bracket )
80: return Bracket_Polynomial is
81:
82: -- DESCRIPTION :
83: -- The labels in the brackets of bp are replaced by those in the minor.
84: -- The bracket polynomial represents the Laplace expansion of that minor.
85:
86: res : Bracket_Polynomial;
87:
88: procedure Substitute_Term ( t : in Bracket_Term; continue : out boolean ) is
89:
90: st : Bracket_Term := Substitute(t,minor);
91:
92: begin
93: Add(res,st);
94: Clear(st);
95: continue := true;
96: end Substitute_Term;
97: procedure Substitute_Terms is new Enumerate_Terms(Substitute_Term);
98:
99: begin
100: Substitute_Terms(bp);
101: return res;
102: end Substitute;
103:
104: -- AUXILIARIES TO Expanded_Minor :
105:
106: function Subtract ( b : Bracket; i : natural ) return Bracket is
107:
108: -- DESCRIPTION :
109: -- Returns a smaller bracket with the ith entry removed.
110:
111: res : Bracket(b'first..b'last-1);
112:
113: begin
114: res(b'first..(i-1)) := b(b'first..(i-1));
115: res(i..res'last) := b((i+1)..b'last);
116: return res;
117: end Subtract;
118:
119: function General_Expanded_Minor ( m : Standard_Complex_Poly_Matrices.Matrix;
120: b : Bracket ) return Poly is
121:
122: -- DESCRIPTION :
123: -- This function treats the case for Expanded_Minor when b'length > 2.
124:
125: res : Poly := Null_Poly;
126: sig : integer;
127: submin,acc : Poly;
128:
129: begin
130: if b'last mod 2 = 0
131: then sig := -1;
132: else sig := +1;
133: end if;
134: for i in b'range loop
135: if m(b(i),b'last) /= Null_Poly
136: then submin := Expanded_Minor(m,Subtract(b,i));
137: if submin /= Null_Poly
138: then acc := m(b(i),b'last)*submin;
139: if sig > 0
140: then Add(res,acc);
141: else Sub(res,acc);
142: end if;
143: Clear(acc);
144: end if;
145: Clear(submin);
146: end if;
147: sig := -sig;
148: end loop;
149: return res;
150: end General_Expanded_Minor;
151:
152: procedure Purify ( p : in out Poly ) is
153:
154: -- DESCRIPTION :
155: -- Eliminates terms that have coefficients less that 10.0**(-10).
156:
157: tol : constant double_float := 10.0**(-10);
158: res : Poly := Null_Poly;
159:
160: procedure Scan_Term ( t : in Term; cont : out boolean ) is
161: begin
162: if AbsVal(t.cf) > tol
163: then Add(res,t);
164: end if;
165: cont := true;
166: end Scan_Term;
167: procedure Scan_Terms is new Visiting_Iterator(Scan_Term);
168:
169: begin
170: Scan_Terms(p);
171: Clear(p);
172: p := res;
173: end Purify;
174:
175: -- TARGET ROUTINES :
176:
177: function Schubert_Pattern ( n : natural; b1,b2 : Bracket )
178: return Standard_Complex_Poly_Matrices.Matrix is
179:
180: res : Standard_Complex_Poly_Matrices.Matrix(1..n,b1'range);
181: d : constant natural := b1'last;
182:
183: begin
184: for i in 1..n loop
185: for j in b1'range loop
186: if ((i < b2(j)) or (i > n+1-b1(b1'last+1-j)))
187: then res(i,j) := Null_Poly;
188: else res(i,j) := Matrix_Indeterminates.Monomial(n,d,i,j);
189: end if;
190: end loop;
191: end loop;
192: return res;
193: end Schubert_Pattern;
194:
195: function Localization_Pattern ( n : natural; top,bottom : Bracket )
196: return Standard_Complex_Poly_Matrices.Matrix is
197:
198: p : constant natural := top'length;
199: res : Standard_Complex_Poly_Matrices.Matrix(1..n,1..p);
200:
201: begin
202: for i in res'range(1) loop
203: for j in res'range(2) loop
204: if i < top(j) or i > bottom(j)
205: then res(i,j) := Null_Poly;
206: else res(i,j) := Matrix_Indeterminates.Monomial(n,p,i,j);
207: end if;
208: end loop;
209: end loop;
210: return res;
211: end Localization_Pattern;
212:
213: -- SYMBOLIC REPRESENTATION OF THE EQUATIONS :
214:
215: function Maximal_Minors ( n,m : natural ) return Bracket_Monomial is
216:
217: res : Bracket_Monomial;
218: first : boolean := true;
219: accu : Bracket(1..m);
220:
221: procedure Enumerate_Brackets ( k,start : natural ) is
222:
223: -- DESCRIPTION :
224: -- Enumerate all m-brackets with entries between 1 and n,
225: -- beginning at start, and currently filling in the kth entry.
226:
227: begin
228: if k > m
229: then if first
230: then res := Create(accu); first := false;
231: else Multiply(res,accu);
232: end if;
233: else for i in start..n loop
234: accu(k) := i;
235: Enumerate_Brackets(k+1,i+1);
236: end loop;
237: end if;
238: end Enumerate_Brackets;
239:
240: begin
241: Enumerate_Brackets(1,1);
242: return res;
243: end Maximal_Minors;
244:
245: function Minor_Equations
246: ( m,d : natural; bm : Bracket_Monomial ) return Bracket_System is
247:
248: res : Bracket_System(0..Number_of_Brackets(bm));
249: gen : Bracket_Polynomial := Laplace_Expansion(m,d);
250: cnt : natural := 0;
251:
252: procedure Substitute_Bracket ( b : in Bracket; continue : out boolean ) is
253: begin
254: cnt := cnt + 1;
255: res(cnt) := Substitute(gen,b);
256: continue := true;
257: end Substitute_Bracket;
258: procedure Substitute_Brackets is
259: new Enumerate_Brackets(Substitute_Bracket);
260:
261: begin
262: res(0) := gen;
263: Substitute_Brackets(bm);
264: return res;
265: end Minor_Equations;
266:
267: function Expanded_Minor ( m : Standard_Complex_Poly_Matrices.Matrix;
268: b : Bracket ) return Poly is
269:
270: res : Poly := Null_Poly;
271: acc : Poly;
272:
273: begin
274: if b'length = 1
275: then Copy(m(b(1),1),res);
276: elsif b'length = 2
277: then if (m(b(1),1) /= Null_Poly) and (m(b(2),2) /= Null_Poly)
278: then res := m(b(1),1)*m(b(2),2);
279: end if;
280: if (m(b(2),1) /= Null_Poly) and (m(b(1),2) /= Null_Poly)
281: then acc := m(b(2),1)*m(b(1),2);
282: Sub(res,acc);
283: Clear(acc);
284: end if;
285: else res := General_Expanded_Minor(m,b);
286: end if;
287: Purify(res);
288: return res;
289: end Expanded_Minor;
290:
291: function Extend_Zero_Lifting ( p : Poly ) return Poly is
292:
293: res : Poly := Null_Poly;
294:
295: procedure Extend_Term ( t : in Term; continue : out boolean ) is
296:
297: et : Term;
298:
299: begin
300: et.dg := new Standard_Natural_Vectors.Vector(t.dg'first..t.dg'last+1);
301: et.dg(t.dg'range) := t.dg.all;
302: et.dg(et.dg'last) := 0;
303: et.cf := t.cf;
304: Add(res,et);
305: Clear(et.dg);
306: continue := true;
307: end Extend_Term;
308: procedure Extend_Terms is new Visiting_Iterator(Extend_Term);
309:
310: begin
311: Extend_Terms(p);
312: return res;
313: end Extend_Zero_Lifting;
314:
315: end Symbolic_Minor_Equations;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to symbolic_minor_equations.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