Annotation of OpenXM_contrib/PHC/Ada/Schubert/determinantal_systems.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 Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
5: with Bracket_Monomials; use Bracket_Monomials;
6: with Bracket_Systems; use Bracket_Systems;
7: with Evaluated_Minors; use Evaluated_Minors;
8: with Symbolic_Minor_Equations; use Symbolic_Minor_Equations;
9: with Numeric_Minor_Equations; use Numeric_Minor_Equations;
10:
11: package body Determinantal_Systems is
12:
13: -- AUXILIARY TO EVALUATION OF MAXIMAL MINORS :
14:
15: function Number_of_Maximal_Minors ( nrows,ncols : natural ) return natural is
16:
17: -- DESCRIPTION :
18: -- Returns the number of maximal minors of a matrix with nrows and ncols.
19:
20: -- REQUIRED : nrows > ncols.
21:
22: begin
23: if ncols = 1
24: then return nrows;
25: else return Number_of_Maximal_Minors(nrows-1,ncols-1)*nrows/ncols;
26: end if;
27: end Number_of_Maximal_Minors;
28:
29: -- LOCALIZATION MAPS :
30:
31: function Localize ( locmap : Standard_Natural_Matrices.Matrix;
32: t : Term ) return Term is
33:
34: -- DESCRIPTION :
35: -- Applies the localization map to the term, eliminating those xij's
36: -- xij for which the corresponding entry in locmap is either 0 or 1.
37:
38: -- NOTE :
39: -- This localization assumes that t.dg(k) = 0 with k for which the
40: -- corresponding (i,j) with locmap(i,j) = 0.
41:
42: res : Term;
43: ndg : Standard_Natural_Vectors.Vector(t.dg'range);
44: cnt : natural := t.dg'first-1;
45: ind : natural := cnt;
46:
47: begin
48: for i in locmap'range(1) loop -- lexicographic order of variables
49: for j in locmap'range(2) loop
50: ind := ind+1;
51: if locmap(i,j) = 2
52: then cnt := cnt + 1;
53: ndg(cnt) := t.dg(ind);
54: end if;
55: end loop;
56: end loop;
57: for i in ind+1..t.dg'last loop -- leave the lifting !
58: cnt := cnt+1;
59: ndg(cnt) := t.dg(i);
60: end loop;
61: res.cf := t.cf;
62: res.dg := new Standard_Natural_Vectors.Vector'(ndg(1..cnt));
63: return res;
64: end Localize;
65:
66: function Localize ( locmap : Standard_Natural_Matrices.Matrix;
67: p : Poly ) return Poly is
68:
69: -- DESCRIPTION :
70: -- Applies the localization map to the polynomial, eliminating
71: -- those xij's for which locmap(i,j) is either 0 or 1.
72:
73: res : Poly := Null_Poly;
74:
75: procedure Localize_Term ( t : in Term; continue : out boolean ) is
76:
77: lt : Term := Localize(locmap,t);
78:
79: begin
80: Add(res,lt);
81: Clear(lt.dg);
82: continue := true;
83: end Localize_Term;
84: procedure Localize_Terms is new Visiting_Iterator(Localize_Term);
85:
86: begin
87: Localize_Terms(p);
88: return res;
89: end Localize;
90:
91: -- TARGET ROUTINES :
92:
93: function Standard_Coordinate_Frame
94: ( x : Standard_Complex_Poly_Matrices.Matrix;
95: plane : Standard_Complex_Matrices.Matrix )
96: return Standard_Natural_Matrices.Matrix is
97:
98: res : Standard_Natural_Matrices.Matrix(x'range(1),x'range(2));
99: tol : constant double_float := 10.0**(-10);
100: first : boolean;
101:
102: begin
103: for j in res'range(2) loop
104: first := true;
105: for i in res'range(1) loop
106: if x(i,j) = Null_Poly
107: then res(i,j) := 0;
108: elsif (first and (AbsVal(plane(i,j)) > tol))
109: then res(i,j) := 1; first := false;
110: else res(i,j) := 2;
111: end if;
112: end loop;
113: end loop;
114: return res;
115: end Standard_Coordinate_Frame;
116:
117: function Maximal_Coordinate_Frame
118: ( x : Standard_Complex_Poly_Matrices.Matrix;
119: plane : Standard_Complex_Matrices.Matrix )
120: return Standard_Natural_Matrices.Matrix is
121:
122: res : Standard_Natural_Matrices.Matrix(x'range(1),x'range(2));
123: max,tmp : double_float;
124: ind : natural;
125:
126: begin
127: for j in res'range(2) loop
128: for i in res'range(1) loop
129: if x(i,j) = Null_Poly
130: then res(i,j) := 0;
131: else res(i,j) := 2;
132: end if;
133: end loop;
134: max := 0.0; ind := 0;
135: for i in res'range(1) loop
136: tmp := AbsVal(plane(i,j));
137: if tmp > max
138: then max := tmp; ind := i;
139: end if;
140: end loop;
141: res(ind,j) := 1;
142: end loop;
143: return res;
144: end Maximal_Coordinate_Frame;
145:
146: function Localize ( locmap : Standard_Natural_Matrices.Matrix;
147: p : Poly_Sys ) return Poly_Sys is
148:
149: res : Poly_Sys(p'range);
150:
151: begin
152: for i in p'range loop
153: res(i) := Localize(locmap,p(i));
154: end loop;
155: return res;
156: end Localize;
157:
158: -- CONSTRUCT THE POLYNOMIAL EQUATIONS :
159:
160: function Polynomial_Equations
161: ( l : Standard_Complex_Matrices.Matrix;
162: x : Standard_Complex_Poly_Matrices.Matrix ) return Poly_Sys is
163:
164: n : constant natural := x'length(1);
165: p : constant natural := x'length(2);
166: kd : constant natural := p + l'length(2);
167: bm : Bracket_Monomial := Maximal_Minors(n,kd);
168: bs : Bracket_System(0..Number_of_Brackets(bm))
169: := Minor_Equations(kd,kd-p,bm);
170: res : Poly_Sys(bs'first+1..bs'last) := Expanded_Minors(l,x,bs);
171:
172: begin
173: Clear(bm); Clear(bs);
174: return res;
175: end Polynomial_Equations;
176:
177: procedure Concat ( l : in out Link_to_Poly_Sys; p : Poly_Sys ) is
178: begin
179: if l = null
180: then declare
181: newsys : Poly_Sys(p'range);
182: cnt : natural := p'first-1;
183: begin
184: for i in p'range loop
185: if p(i) /= Null_Poly
186: then cnt := cnt+1;
187: newsys(cnt) := p(i);
188: end if;
189: end loop;
190: l := new Poly_Sys'(newsys(p'first..cnt));
191: end;
192: else declare
193: newsys : Poly_Sys(l'first..l'last+p'length);
194: ind : natural := l'last;
195: begin
196: newsys(l'range) := l(l'range);
197: for i in p'range loop
198: if p(i) /= Null_Poly
199: then ind := ind+1;
200: newsys(ind) := p(i);
201: end if;
202: end loop;
203: l := new Poly_Sys'(newsys(l'first..ind));
204: end;
205: end if;
206: end Concat;
207:
208: function Polynomial_Equations
209: ( l : Standard_Complex_VecMats.VecMat;
210: x : Standard_Complex_Poly_Matrices.Matrix ) return Poly_Sys is
211:
212: n : constant natural := x'length(1);
213: p : constant natural := x'length(2);
214: kd : constant natural := p + l(l'first)'length(2);
215: bm : Bracket_Monomial := Maximal_Minors(n,kd);
216: bs : Bracket_System(0..Number_of_Brackets(bm))
217: := Minor_Equations(kd,kd-p,bm);
218: nb : constant natural := l'length;
219: res : Poly_Sys(bs'first+1..nb*bs'last);
220: cnt : natural := bs'first;
221:
222: begin
223: for i in 1..nb loop
224: declare
225: sys : constant Poly_Sys := Expanded_Minors(l(i).all,x,bs);
226: begin
227: for j in sys'range loop
228: cnt := cnt + 1;
229: res(cnt) := sys(j);
230: end loop;
231: end;
232: end loop;
233: Clear(bm); Clear(bs);
234: return res;
235: end Polynomial_Equations;
236:
237: -- EVALUATORS AND DIFFERENTIATORS :
238:
239: function Eval ( l,x : Standard_Complex_Matrices.Matrix )
240: return Complex_Number is
241:
242: wrk : Standard_Complex_Matrices.Matrix(x'range(1),x'range(1));
243:
244: begin
245: for j in l'range(2) loop
246: for i in l'range(1) loop
247: wrk(i,j) := l(i,j);
248: end loop;
249: end loop;
250: for j in x'range(2) loop
251: for i in x'range(1) loop
252: wrk(i,l'last(2)+j) := x(i,j);
253: end loop;
254: end loop;
255: return Determinant(wrk);
256: end Eval;
257:
258: function Eval ( l,x : Standard_Complex_Matrices.Matrix ) return Vector is
259:
260: n : constant natural := l'length(2) + x'length(2);
261: dim : constant natural := Number_of_Maximal_Minors(l'length(1),n);
262: res : Standard_Complex_Vectors.Vector(1..dim);
263: -- := (1..dim => Create(0.0));
264: cnt : natural := 0;
265: rws : Standard_Natural_Vectors.Vector(1..n);
266: wrk : Standard_Complex_Matrices.Matrix(1..n,1..n);
267:
268: procedure Choose_next_Row ( k,start : in natural ) is
269:
270: -- DESCRIPTION :
271: -- Chooses next k-th row within the range start..l'last(1), if k <= n.
272: -- If k > n, then the minor defined by the current row selection
273: -- is computed and added to the result.
274:
275: begin
276: if k > n
277: then for j in l'range(2) loop -- select rows
278: for i in 1..n loop
279: wrk(i,j) := l(rws(i),j);
280: end loop;
281: end loop;
282: for j in x'range(2) loop
283: for i in 1..n loop
284: wrk(i,l'last(2)+j) := x(rws(i),j);
285: end loop;
286: end loop;
287: cnt := cnt + 1;
288: res(cnt) := Determinant(wrk);
289: else for i in start..l'last(1) loop
290: rws(k) := i;
291: Choose_next_Row(k+1,i+1);
292: end loop;
293: end if;
294: end Choose_next_Row;
295:
296: begin
297: Choose_next_Row(1,1);
298: return res;
299: end Eval;
300:
301: function Minor ( l,x : Standard_Complex_Matrices.Matrix; row,col : natural )
302: return Complex_Number is
303:
304: -- DESCRIPTION :
305: -- Returns the minor [l|x] with row and column removed and with
306: -- the appropriate sign.
307:
308: wrk : Standard_Complex_Matrices.Matrix
309: (x'first(1)..x'last(1)-1,x'first(1)..x'last(1)-1);
310: ind : natural;
311:
312: begin
313: for j in l'range(2) loop
314: for i in l'range(1) loop
315: if i < row
316: then wrk(i,j) := l(i,j);
317: elsif i > row
318: then wrk(i-1,j) := l(i,j);
319: end if;
320: end loop;
321: end loop;
322: for j in x'range(2) loop
323: if j < col
324: then ind := l'last(2) + j;
325: elsif j > col
326: then ind := l'last(2) + j - 1;
327: end if;
328: if j /= col
329: then for i in x'range(1) loop
330: if i < row
331: then wrk(i,ind) := x(i,j);
332: elsif i > row
333: then wrk(i-1,ind) := x(i,j);
334: end if;
335: end loop;
336: end if;
337: end loop;
338: if (row + l'last(2) + col) mod 2 = 0
339: then return Determinant(wrk);
340: else return -Determinant(wrk);
341: end if;
342: end Minor;
343:
344: function Diff ( l,x : Standard_Complex_Matrices.Matrix; i : natural )
345: return Complex_Number is
346:
347: p : constant natural := x'length(2);
348: row,col : natural;
349:
350: begin
351: row := i/p + 1;
352: col := i mod p;
353: if col = 0
354: then col := p;
355: row := row - 1;
356: end if;
357: return Minor(l,x,row,col);
358: end Diff;
359:
360: function Diff ( l,x : Standard_Complex_Matrices.Matrix;
361: locmap : Standard_Natural_Matrices.Matrix; i : natural )
362: return Complex_Number is
363:
364: row,col : natural;
365: cnt : natural := 0;
366:
367: begin
368: for k1 in locmap'range(1) loop
369: for k2 in locmap'range(2) loop
370: if locmap(k1,k2) = 2
371: then cnt := cnt+1;
372: if cnt = i
373: then row := k1;
374: col := k2;
375: end if;
376: end if;
377: exit when (cnt = i);
378: end loop;
379: exit when (cnt = i);
380: end loop;
381: return Minor(l,x,row,col);
382: end Diff;
383:
384: function Eval ( l : Standard_Complex_VecMats.VecMat;
385: x : Standard_Complex_Matrices.Matrix )
386: return Standard_Complex_Vectors.Vector is
387:
388: res : Standard_Complex_Vectors.Vector(l'range);
389:
390: begin
391: for i in l'range loop
392: res(i) := Eval(l(i).all,x);
393: end loop;
394: return res;
395: end Eval;
396:
397: function Diff ( l : Standard_Complex_VecMats.VecMat;
398: x : Standard_Complex_Matrices.Matrix )
399: return Standard_Complex_Matrices.Matrix is
400:
401: res : Standard_Complex_Matrices.Matrix(l'range,1..x'last(1)*x'last(2));
402:
403: begin
404: for i in res'range(1) loop
405: for j in res'range(2) loop
406: res(i,j) := Diff(l(i).all,x,j);
407: end loop;
408: end loop;
409: return res;
410: end Diff;
411:
412: function Old_Diff ( l : Standard_Complex_VecMats.VecMat;
413: x : Standard_Complex_Matrices.Matrix; nvars : natural;
414: locmap : Standard_Natural_Matrices.Matrix )
415: return Standard_Complex_Matrices.Matrix is
416:
417: res : Standard_Complex_Matrices.Matrix(l'range,1..nvars);
418:
419: begin
420: for i in res'range(1) loop
421: for j in res'range(2) loop
422: res(i,j) := Diff(l(i).all,x,locmap,j);
423: end loop;
424: end loop;
425: return res;
426: end Old_Diff;
427:
428: function Diff ( l : Standard_Complex_VecMats.VecMat;
429: x : Standard_Complex_Matrices.Matrix; nvars : natural;
430: locmap : Standard_Natural_Matrices.Matrix )
431: return Standard_Complex_Matrices.Matrix is
432:
433: -- NOTE :
434: -- This Diff is organized according to the localization map,
435: -- to avoid multiple searches.
436:
437: res : Standard_Complex_Matrices.Matrix(l'range,1..nvars);
438: ind : natural := 0;
439:
440: begin
441: for i in locmap'range loop
442: for j in locmap'range loop
443: if locmap(i,j) = 2
444: then ind := ind+1;
445: for k in l'range loop
446: res(k,ind) := Minor(l(k).all,x,i,j);
447: end loop;
448: end if;
449: end loop;
450: end loop;
451: return res;
452: end Diff;
453:
454: end Determinantal_Systems;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to determinantal_systems.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