Annotation of OpenXM_contrib/PHC/Ada/Schubert/curves_into_grassmannian.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
2: with Standard_Natural_Vectors;
3: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
4:
5: package body Curves_into_Grassmannian is
6:
7: -- CREATOR :
8:
9: function Symbolic_Create ( m,p,q : natural; top,bottom : Bracket )
10: return Standard_Complex_Poly_Matrices.Matrix is
11:
12: res : Standard_Complex_Poly_Matrices.Matrix(1..(m+p),1..p);
13: rws : constant natural := (m+p)*(q+1);
14: n : constant natural := Number_of_Variables(top,bottom) + 2;
15: row,ind,s_deg,t_deg : natural;
16: t : Term;
17:
18: begin
19: for i in res'range(1) loop -- initialization
20: for j in res'range(2) loop
21: res(i,j) := Null_Poly;
22: end loop;
23: end loop;
24: t.dg := new Standard_Natural_Vectors.Vector'(1..n => 0);
25: t.cf := Create(1.0);
26: ind := 0; -- ind counts #variables
27: for j in 1..p loop -- assign columnwise
28: t_deg := (bottom(j)-1)/(m+p); -- degree in t to homogenize
29: row := 0; s_deg := 0;
30: for i in 1..rws loop
31: row := row + 1;
32: if i >= top(j) and i <= bottom(j)
33: then ind := ind+1;
34: t.dg(n-1) := s_deg; t.dg(n) := t_deg;
35: t.dg(ind) := 1;
36: Add(res(row,j),t);
37: t.dg(n-1) := 0; t.dg(n) := 0;
38: t.dg(ind) := 0;
39: end if;
40: if i mod (m+p) = 0
41: then row := 0; s_deg := s_deg+1;
42: if t_deg > 0
43: then t_deg := t_deg-1;
44: end if;
45: end if;
46: end loop;
47: end loop;
48: Clear(t);
49: return res;
50: end Symbolic_Create;
51:
52: -- SELECTORS :
53:
54: function Number_of_Variables ( top,bottom : Bracket ) return natural is
55:
56: cnt : natural := 0;
57:
58: begin
59: for j in top'range loop
60: cnt := cnt + (bottom(j) - top(j) + 1);
61: end loop;
62: return cnt;
63: end Number_of_Variables;
64:
65: function Standard_Coordinate_Frame
66: ( m,p,q : natural; top,bottom : Bracket;
67: coeff : Standard_Complex_Matrices.Matrix )
68: return Standard_Natural_Matrices.Matrix is
69:
70: rws : constant natural := (m+p)*(q+1);
71: res : Standard_Natural_Matrices.Matrix(1..rws,1..p);
72: tol : constant double_float := 10.0**(-10);
73: first : boolean;
74:
75: begin
76: for j in 1..p loop
77: first := true;
78: for i in 1..rws loop
79: if i < top(j) or i > bottom(j)
80: then res(i,j) := 0;
81: elsif (first and (AbsVal(coeff(i,j)) > tol))
82: then res(i,j) := 1; first := false;
83: else res(i,j) := 2;
84: end if;
85: end loop;
86: end loop;
87: return res;
88: end Standard_Coordinate_Frame;
89:
90: function Eval ( c : Term; s,t : Complex_Number ) return Term is
91:
92: res : Term;
93:
94: begin
95: Copy(c,res);
96: for i in 1..res.dg(res.dg'last-1) loop -- evaluate s
97: res.cf := res.cf*s;
98: end loop;
99: res.dg(res.dg'last-1) := 0;
100: for i in 1..res.dg(res.dg'last) loop -- evaluate t
101: res.cf := res.cf*t;
102: end loop;
103: res.dg(res.dg'last) := 0;
104: return res;
105: end Eval;
106:
107: function Eval ( p : Poly; s,t : Complex_Number ) return Poly is
108:
109: res : Poly := Null_Poly;
110:
111: procedure Eval_Term ( ct : in Term; continue : out boolean ) is
112:
113: et : Term := Eval(ct,s,t);
114:
115: begin
116: Add(res,et);
117: Clear(et);
118: continue := true;
119: end Eval_Term;
120: procedure Eval_Terms is new Visiting_Iterator(Eval_Term);
121:
122: begin
123: Eval_Terms(p);
124: return res;
125: end Eval;
126:
127: function Eval ( c : Standard_Complex_Poly_Matrices.Matrix;
128: s,t : Complex_Number )
129: return Standard_Complex_Poly_Matrices.Matrix is
130:
131: res : Standard_Complex_Poly_Matrices.Matrix(c'range(1),c'range(2));
132:
133: begin
134: for i in c'range(1) loop
135: for j in c'range(2) loop
136: if c(i,j) = Null_Poly
137: then res(i,j) := Null_Poly;
138: else res(i,j) := Eval(c(i,j),s,t);
139: end if;
140: end loop;
141: end loop;
142: return res;
143: end Eval;
144:
145: function Elim ( c : Term; s,t : Complex_Number ) return Term is
146:
147: res : Term;
148:
149: begin
150: res.dg := new Standard_Natural_Vectors.Vector'
151: (c.dg(c.dg'first..c.dg'last-2));
152: res.cf := c.cf;
153: for i in 1..c.dg(c.dg'last-1) loop -- evaluate s
154: res.cf := res.cf*s;
155: end loop;
156: for i in 1..c.dg(c.dg'last) loop -- evaluate t
157: res.cf := res.cf*t;
158: end loop;
159: return res;
160: end Elim;
161:
162: function Elim ( p : Poly; s,t : Complex_Number ) return Poly is
163:
164: res : Poly := Null_Poly;
165:
166: procedure Elim_Term ( ct : in Term; continue : out boolean ) is
167:
168: et : Term := Elim(ct,s,t);
169:
170: begin
171: Add(res,et);
172: Clear(et);
173: continue := true;
174: end Elim_Term;
175: procedure Elim_Terms is new Visiting_Iterator(Elim_Term);
176:
177: begin
178: Elim_Terms(p);
179: return res;
180: end Elim;
181:
182: function Elim ( c : Standard_Complex_Poly_matrices.Matrix;
183: s,t : Complex_Number )
184: return Standard_Complex_Poly_Matrices.Matrix is
185:
186: res : Standard_Complex_Poly_Matrices.Matrix(c'range(1),c'range(2));
187:
188: begin
189: for i in c'range(1) loop
190: for j in c'range(2) loop
191: if c(i,j) = Null_Poly
192: then res(i,j) := Null_Poly;
193: else res(i,j) := Elim(c(i,j),s,t);
194: end if;
195: end loop;
196: end loop;
197: return res;
198: end Elim;
199:
200: function Column_Localize ( top,bottom : Bracket;
201: locmap : Standard_Natural_Matrices.Matrix;
202: t : Term ) return Term is
203:
204: -- DESCRIPTION :
205: -- Applies the localization map to the term, eliminating those xij's
206: -- xij for which the corresponding entry in locmap is either 0 or 1.
207:
208: -- NOTE :
209: -- This localization assumes that t.dg(k) = 0 with k for which the
210: -- corresponding (i,j) with locmap(i,j) = 0.
211: -- The localization pattern is traversed columnwise.
212:
213: res : Term;
214: ndg : Standard_Natural_Vectors.Vector(t.dg'range);
215: cnt : natural := t.dg'first-1;
216: ind : natural := cnt;
217:
218: begin
219: for j in locmap'range(2) loop -- columnwise order of the variables
220: for i in top(j)..bottom(j) loop -- restricted range skips the zeros
221: ind := ind+1;
222: if locmap(i,j) = 2 -- skip the ones
223: then cnt := cnt + 1;
224: ndg(cnt) := t.dg(ind);
225: end if;
226: end loop;
227: end loop;
228: for i in ind+1..t.dg'last loop -- leave the lifting !
229: cnt := cnt+1;
230: ndg(cnt) := t.dg(i);
231: end loop;
232: res.cf := t.cf;
233: res.dg := new Standard_Natural_Vectors.Vector'(ndg(1..cnt));
234: return res;
235: end Column_Localize;
236:
237: function Column_Localize ( top,bottom : Bracket;
238: locmap : Standard_Natural_Matrices.Matrix;
239: p : Poly ) return Poly is
240:
241: -- DESCRIPTION :
242: -- Applies the localization map to the polynomial, eliminating
243: -- those xij's for which locmap(i,j) is either 0 or 1.
244:
245: res : Poly := Null_Poly;
246:
247: procedure Column_Localize_Term ( t : in Term; continue : out boolean ) is
248:
249: lt : Term := Column_Localize(top,bottom,locmap,t);
250:
251: begin
252: Add(res,lt);
253: Clear(lt.dg);
254: continue := true;
255: end Column_Localize_Term;
256: procedure Column_Localize_Terms is
257: new Visiting_Iterator(Column_Localize_Term);
258:
259: begin
260: Column_Localize_Terms(p);
261: return res;
262: end Column_Localize;
263:
264: function Column_Localize ( top,bottom : Bracket;
265: locmap : Standard_Natural_Matrices.Matrix;
266: p : Poly_Sys ) return Poly_Sys is
267:
268: res : Poly_Sys(p'range);
269:
270: begin
271: for i in p'range loop
272: res(i) := Column_Localize(top,bottom,locmap,p(i));
273: end loop;
274: return res;
275: end Column_Localize;
276:
277: function Column_Vector_Rep
278: ( top,bottom : Bracket;
279: cffmat : Standard_Complex_Matrices.Matrix )
280: return Standard_Complex_Vectors.Vector is
281:
282: dim : constant natural := Number_of_Variables(top,bottom);
283: res : Standard_Complex_Vectors.Vector(1..dim);
284: cnt : natural := 0;
285:
286: begin
287: for j in cffmat'range(2) loop
288: for i in top(j)..bottom(j) loop
289: cnt := cnt + 1;
290: res(cnt) := cffmat(i,j);
291: end loop;
292: end loop;
293: return res;
294: end Column_Vector_Rep;
295:
296: function Column_Vector_Rep ( locmap : Standard_Natural_Matrices.Matrix;
297: cffmat : Standard_Complex_Matrices.Matrix )
298: return Standard_Complex_Vectors.Vector is
299:
300: dim : constant natural := cffmat'length(1)*cffmat'length(2);
301: res : Standard_Complex_Vectors.Vector(1..dim);
302: cnt : natural := 0;
303:
304: begin
305: for j in cffmat'range(2) loop
306: for i in cffmat'range(1) loop
307: if locmap(i,j) = 2
308: then cnt := cnt + 1;
309: res(cnt) := cffmat(i,j);
310: end if;
311: end loop;
312: end loop;
313: return res(1..cnt);
314: end Column_Vector_Rep;
315:
316: function Column_Matrix_Rep
317: ( locmap : Standard_Natural_Matrices.Matrix;
318: cffvec : Standard_Complex_Vectors.Vector )
319: return Standard_Complex_Matrices.Matrix is
320:
321: res : Standard_Complex_Matrices.Matrix(locmap'range(1),locmap'range(2));
322: cnt : natural := 0;
323:
324: begin
325: for j in locmap'range(2) loop
326: for i in locmap'range(1) loop
327: if locmap(i,j) = 0
328: then res(i,j) := Create(0.0);
329: elsif locmap(i,j) = 1
330: then res(i,j) := Create(1.0);
331: else cnt := cnt + 1;
332: res(i,j) := cffvec(cnt);
333: end if;
334: end loop;
335: end loop;
336: return res;
337: end Column_Matrix_Rep;
338:
339: procedure Swap ( p : in out Poly; k,l : in natural ) is
340:
341: procedure Swap_Term ( t : in out Term; continue : out boolean ) is
342:
343: lval : natural := t.dg(l);
344:
345: begin
346: t.dg(l) := t.dg(k);
347: t.dg(k) := lval;
348: continue := true;
349: end Swap_Term;
350: procedure Swap_Terms is new Changing_Iterator(Swap_Term);
351:
352: begin
353: Swap_Terms(p);
354: end Swap;
355:
356: procedure Swap ( c : in out Standard_Complex_Poly_Matrices.Matrix;
357: k,l : in natural ) is
358: begin
359: for i in c'range(1) loop
360: for j in c'range(2) loop
361: if c(i,j) /= Null_Poly
362: then Swap(c(i,j),k,l);
363: end if;
364: end loop;
365: end loop;
366: end Swap;
367:
368: function Insert ( p : Poly; k : natural ) return Poly is
369:
370: res : Poly := Null_Poly;
371:
372: procedure Insert_Term ( t : in Term; continue : out boolean ) is
373:
374: rt : Term;
375:
376: begin
377: rt.cf := t.cf;
378: rt.dg := new Standard_Natural_Vectors.Vector(t.dg'first..t.dg'last+1);
379: rt.dg(t.dg'first..k-1) := t.dg(t.dg'first..k-1);
380: rt.dg(k) := 0;
381: rt.dg(k+1..rt.dg'last) := t.dg(k..t.dg'last);
382: Add(res,rt);
383: Clear(rt);
384: continue := true;
385: end Insert_Term;
386: procedure Insert_Terms is new Visiting_Iterator(Insert_Term);
387:
388: begin
389: Insert_Terms(p);
390: return res;
391: end Insert;
392:
393: function Insert ( c : Standard_Complex_Poly_Matrices.Matrix; k : natural )
394: return Standard_Complex_Poly_Matrices.Matrix is
395:
396: res : Standard_Complex_Poly_Matrices.Matrix(c'range(1),c'range(2));
397:
398: begin
399: for i in c'range(1) loop
400: for j in c'range(2) loop
401: if c(i,j) = Null_Poly
402: then res(i,j) := Null_Poly;
403: else res(i,j) := Insert(c(i,j),k);
404: end if;
405: end loop;
406: end loop;
407: return res;
408: end Insert;
409:
410: procedure Duplicate ( p : in out Poly; k,l : in natural ) is
411:
412: procedure Duplicate_Term ( t : in out Term; continue : out boolean ) is
413: begin
414: t.dg(l) := t.dg(k);
415: continue := true;
416: end Duplicate_Term;
417: procedure Duplicate_Terms is new Changing_Iterator(Duplicate_Term);
418:
419: begin
420: Duplicate_Terms(p);
421: end Duplicate;
422:
423: end Curves_into_Grassmannian;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to curves_into_grassmannian.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