Annotation of OpenXM_contrib/PHC/Ada/Schubert/pieri_homotopies.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
2: with Standard_Random_Numbers; use Standard_Random_Numbers;
3: with Standard_Natural_Vectors;
4: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
5: with Numeric_Minor_Equations; use Numeric_Minor_Equations;
6: with Determinantal_Systems; use Determinantal_Systems;
7: with Specialization_of_Planes; use Specialization_of_Planes;
8: with Curves_into_Grassmannian; use Curves_into_Grassmannian;
9:
10: package body Pieri_Homotopies is
11:
12: -- AUXILIARIES TO THE QUANTUM CASE :
13:
14: procedure Multiply ( p : in out Poly; var,deg : natural ) is
15:
16: -- DESCRIPTION :
17: -- Multiplies p with x(var)**deg.
18:
19: procedure Multiply_Term ( t : in out Term; continue : out boolean ) is
20: begin
21: t.dg(var) := t.dg(var) + deg;
22: continue := true;
23: end Multiply_Term;
24: procedure Multiply_Terms is new Changing_Iterator(Multiply_Term);
25:
26: begin
27: Multiply_Terms(p);
28: end Multiply;
29:
30: procedure Multiply ( xpm : in out Standard_Complex_Poly_Matrices.Matrix;
31: col,var,deg : in natural ) is
32:
33: -- DESCRIPTION :
34: -- Multiplies all polynomial in the column col of xpm with x(var)**deg.
35:
36: begin
37: for i in xpm'range(1) loop
38: if xpm(i,col) /= Null_Poly
39: then Multiply(xpm(i,col),var,deg);
40: end if;
41: end loop;
42: end Multiply;
43:
44: procedure Add ( p : in out Poly;
45: var : in natural; start : in Complex_Number ) is
46:
47: -- DESCRIPTION :
48: -- Performs p := p + (1-s)*c, where s = x(var) and c = start.
49:
50: n : constant natural := Number_of_Unknowns(p);
51: t : Term;
52:
53: begin
54: t.dg := new Standard_Natural_Vectors.Vector'(1..n => 0);
55: t.cf := start;
56: Add(p,t); -- p = p + c
57: t.dg(var) := 1;
58: Sub(p,t); -- p = p + c - c*s = p + c*(1-s)
59: Clear(t);
60: end Add;
61:
62: procedure Add ( xpm : in out Standard_Complex_Poly_Matrices.Matrix;
63: col,var : in natural; start : in Vector ) is
64:
65: -- DESCRIPTION :
66: -- Adds to every nonzero polynomial in the column col of xpm
67: -- the term c*(1-s), where s = x(var) and c = start(i).
68:
69: begin
70: for i in xpm'range(1) loop
71: if xpm(i,col) /= Null_Poly
72: then Add(xpm(i,col),var,start(i));
73: end if;
74: end loop;
75: end Add;
76:
77: function Degree1 ( nd : Node; n : natural ) return natural is
78:
79: -- DESCRIPTION :
80: -- Returns the maximum of one and the degree of the first column
81: -- of the map that fits the pattern as decribed by pivots in the node.
82: -- The parameter n must equal m+p.
83:
84: d : constant natural := (nd.bottom(1) - nd.top(1))/n;
85:
86: begin
87: if d = 0
88: then return 1;
89: else return d;
90: end if;
91: end Degree1;
92:
93: function Moving_Parameter0 ( n,xk,tk : natural;
94: start,target : Complex_Number ) return Poly is
95:
96: -- DESCRIPTION :
97: -- Returns the equation (x-start)*(1-t) + (x-target)*t = 0 that describes
98: -- the motion of x from start to target as t goes from 0 to 1.
99: -- Note that this equation equals x + (start-target)*t - start = 0.
100: -- This is the older version without using a random constant.
101:
102: -- ON ENTRY :
103: -- n total number of variables, continuation parameter t included;
104: -- xk index of the moving variable x;
105: -- tk index of the continuation parameter
106: -- start starting value for x;
107: -- target target value for x.
108:
109: res : Poly;
110: t : Term;
111:
112: begin
113: t.dg := new Standard_Natural_Vectors.Vector'(1..n => 0);
114: t.cf := -start;
115: res := Create(t); -- res = -start
116: t.cf := Create(1.0);
117: t.dg(xk) := 1;
118: Add(res,t); -- res = x - start
119: t.dg(xk) := 0;
120: t.dg(tk) := 1;
121: t.cf := start - target; -- res = x + (start - target)*t - start
122: Add(res,t);
123: Clear(t);
124: return res;
125: end Moving_Parameter0;
126:
127: function Moving_Parameter ( n,xk,tk : natural;
128: start,target : Complex_Number ) return Poly is
129:
130: -- DESCRIPTION :
131: -- Returns the equation c*(x-start)*(1-t) + (x-target)*t = 0 for
132: -- the motion of x from start to target as t goes from 0 to 1.
133: -- This version uses a constant c, randomly generated from within.
134:
135: -- ON ENTRY :
136: -- n total number of variables, continuation parameter t included;
137: -- xk index of the moving variable x;
138: -- tk index of the continuation parameter t;
139: -- start starting value for x;
140: -- target target value for x.
141:
142: res : Poly;
143: t : Term;
144: c : Complex_Number := Random1;
145:
146: begin
147: t.dg := new Standard_Natural_Vectors.Vector'(1..n => 0);
148: t.cf := -c*start;
149: res := Create(t); -- res = -c*start
150: t.cf := -t.cf;
151: t.dg(tk) := 1;
152: Add(res,t); -- res = -c*start + c*start*t = -c*start*(1-t)
153: t.cf := c;
154: t.dg(tk) := 0;
155: t.dg(xk) := 1;
156: Add(res,t); -- res = -c*start*(1-t) + c*x
157: t.cf := -t.cf;
158: t.dg(tk) := 1;
159: Add(res,t); -- res = -c*start*(1-t) + c*x - c*x*t
160: t.cf := Create(1.0); -- = -c*start*(1-t) + c*x*(1-t)
161: Add(res,t); -- = c*(1-t)*(x-start)
162: t.cf := -target; -- res = c*(1-t)*(x-start) + x*t
163: t.dg(xk) := 0;
164: Add(res,t); -- res = c*(1-t)*(x-start) + x*t - target*t
165: Clear(t); -- = c*(1-t)*(x-start) + (x-target)*t
166: return res;
167: end Moving_Parameter;
168:
169: function Constant_Parameter ( n,i : natural; s : Complex_Number )
170: return Poly is
171:
172: -- DESCRIPTION :
173: -- Returns the polynomial x_i - s; n equals the number of variables.
174:
175: res : Poly;
176: t : Term;
177:
178: begin
179: t.dg := new Standard_Natural_Vectors.Vector'(1..n => 0);
180: t.dg(i) := 1;
181: t.cf := Create(1.0);
182: res := Create(t);
183: t.dg(i) := 0;
184: t.cf := -s;
185: Add(res,t);
186: Clear(t);
187: return res;
188: end Constant_Parameter;
189:
190: -- TARGET ROUTINES :
191:
192: function One_Hypersurface_Pieri_Homotopy
193: ( n : natural; nd : Node; expbp : Bracket_Polynomial;
194: xpm : Standard_Complex_Poly_Matrices.Matrix;
195: planes : VecMat ) return Poly_Sys is
196:
197: res : Poly_Sys(1..nd.level);
198: p : constant natural := nd.p;
199: m : constant natural := n-p;
200: special : Standard_Complex_Matrices.Matrix(1..n,1..m);
201: target,start : Poly;
202:
203: begin
204: case nd.tp is
205: when top => special := Special_Plane(m,nd.top);
206: -- special := Special_Top_Plane(m,nd.top);
207: when bottom => special := Special_Plane(m,nd.bottom);
208: -- special := Special_Bottom_Plane(m,nd.bottom);
209: when others => null; -- mixed case treated separately
210: end case;
211: for i in 1..nd.level-1 loop
212: res(i) := Expanded_Minors(planes(i).all,xpm,expbp);
213: Embed(res(i));
214: end loop;
215: target := Expanded_Minors(planes(nd.level).all,xpm,expbp);
216: start := Expanded_Minors(special,xpm,expbp);
217: res(nd.level) := Linear_Homotopy(target,start);
218: Clear(target); Clear(start);
219: return res;
220: end One_Hypersurface_Pieri_Homotopy;
221:
222: function Two_Hypersurface_Pieri_Homotopy
223: ( n : natural; nd : Node; expbp : Bracket_Polynomial;
224: xpm : Standard_Complex_Poly_Matrices.Matrix;
225: planes : VecMat ) return Poly_Sys is
226:
227: res : Poly_Sys(1..nd.level);
228: p : constant natural := nd.p;
229: m : constant natural := n-p;
230: top_special : constant Standard_Complex_Matrices.Matrix
231: := Special_Plane(m,nd.top);
232: -- := Special_Top_Plane(m,nd.top);
233: bot_special : constant Standard_Complex_Matrices.Matrix
234: := Special_Plane(m,nd.bottom);
235: -- := Special_Bottom_Plane(m,nd.bottom);
236: top_start,bot_start,top_target,bot_target : Poly;
237:
238: begin
239: for i in 1..nd.level-2 loop
240: res(i) := Expanded_Minors(planes(i).all,xpm,expbp);
241: Embed(res(i));
242: end loop;
243: top_target := Expanded_Minors(planes(nd.level).all,xpm,expbp);
244: top_start := Expanded_Minors(top_special,xpm,expbp);
245: res(nd.level) := Linear_Homotopy(top_target,top_start);
246: Clear(top_target); Clear(top_start);
247: bot_target := Expanded_Minors(planes(nd.level-1).all,xpm,expbp);
248: bot_start := Expanded_Minors(bot_special,xpm,expbp);
249: res(nd.level-1) := Linear_Homotopy(bot_target,bot_start);
250: Clear(bot_target); Clear(bot_start);
251: return res;
252: end Two_Hypersurface_Pieri_Homotopy;
253:
254: function One_General_Pieri_Homotopy
255: ( n,ind : natural; nd : Node; bs : Bracket_System;
256: start,target : Standard_Complex_Matrices.Matrix;
257: xpm : Standard_Complex_Poly_Matrices.Matrix;
258: planes : VecMat ) return Link_to_Poly_Sys is
259:
260: res : Link_to_Poly_Sys;
261: nva : constant natural := n*nd.p + 1;
262: moving_plane : Standard_Complex_Poly_Matrices.Matrix(1..n,target'range(2))
263: := Moving_U_Matrix(nva,start,target);
264: moving : Poly_Sys(1..bs'last);
265:
266: begin
267: for i in 1..ind-1 loop
268: Concat(res,Polynomial_Equations(planes(i).all,xpm));
269: end loop;
270: if res /= null
271: then for i in res'range loop
272: Embed(res(i));
273: end loop;
274: end if;
275: moving := Lifted_Expanded_Minors(moving_plane,xpm,bs);
276: Concat(res,moving);
277: Standard_Complex_Poly_Matrices.Clear(moving_plane);
278: return res;
279: end One_General_Pieri_Homotopy;
280:
281: function Two_General_Pieri_Homotopy
282: ( n,ind : natural; nd : Node; top_bs,bot_bs : Bracket_System;
283: top_start,top_target,bot_start,bot_target
284: : Standard_Complex_Matrices.Matrix;
285: xpm : Standard_Complex_Poly_Matrices.Matrix;
286: planes : VecMat ) return Link_to_Poly_Sys is
287:
288: -- DESCRIPTION :
289: -- Returns the homotopy for general linear subspace intersections,
290: -- in case nd.tp = mixed. The parameter ind indicates the plane
291: -- planes(ind) towards the constructed homotopy works.
292:
293: res : Link_to_Poly_Sys;
294: nva : constant natural := n*nd.p + 1;
295: top_moving : Standard_Complex_Poly_Matrices.Matrix(1..n,top_target'range(2))
296: := Moving_U_Matrix(nva,top_start,top_target);
297: bot_moving : Standard_Complex_Poly_Matrices.Matrix(1..n,bot_target'range(2))
298: := Moving_U_Matrix(nva,bot_start,bot_target);
299: top_movsys : Poly_Sys(1..top_bs'last);
300: bot_movsys : Poly_Sys(1..bot_bs'last);
301:
302: begin
303: for i in 1..ind-1 loop
304: Concat(res,Polynomial_Equations(planes(i).all,xpm));
305: end loop;
306: if res /= null
307: then for i in res'range loop
308: Embed(res(i));
309: end loop;
310: end if;
311: top_movsys := Lifted_Expanded_Minors(top_moving,xpm,top_bs);
312: bot_movsys := Lifted_Expanded_Minors(bot_moving,xpm,bot_bs);
313: Concat(res,top_movsys);
314: Concat(res,bot_movsys);
315: Standard_Complex_Poly_Matrices.Clear(top_moving);
316: Standard_Complex_Poly_Matrices.Clear(bot_moving);
317: return res;
318: end Two_General_Pieri_Homotopy;
319:
320: function One_Quantum_Pieri_Homotopy
321: ( n : natural; nd : Node; expbp : Bracket_Polynomial;
322: xpm : Standard_Complex_Poly_Matrices.Matrix;
323: planes : VecMat; s : Vector ) return Poly_Sys is
324:
325: -- DESCRIPTION :
326: -- Returns the Pieri homotopy that corresponds to the node.
327: -- This homotopy is set up to work only in the hypersurface case,
328: -- when the type of the node is either top or bottom.
329:
330: res : Poly_Sys(1..nd.level+1);
331: p : constant natural := nd.p;
332: m : constant natural := n-p;
333: nvars : constant natural := nd.level + p + 2;
334: -- nvars = level #equations in the x_ij's
335: -- + p because not yet fixed the ones
336: -- + 2 for s and t, note that t is continuation parameter
337: special : Standard_Complex_Matrices.Matrix(1..n,1..m);
338: target,start : Poly;
339: map : Standard_Complex_Poly_Matrices.Matrix(1..n,1..p);
340:
341: begin
342: case nd.tp is
343: when top => special := Special_Plane(m,Modulo(nd.top,m+p));
344: Standard_Complex_Poly_Matrices.Copy(xpm,map);
345: Swap(map,nvars-1,nvars);
346: when bottom => special := Special_Plane(m,Modulo(nd.bottom,m+p));
347: map := xpm;
348: when others => null; -- mixed case treated separately
349: end case;
350: for i in 1..nd.level-1 loop
351: declare
352: eva : Standard_Complex_Poly_Matrices.Matrix(xpm'range(1),xpm'range(2))
353: := Eval(xpm,s(i),Create(1.0));
354: begin
355: res(i) := Expanded_Minors(planes(i).all,eva,expbp);
356: Standard_Complex_Poly_Matrices.Clear(eva);
357: end;
358: end loop;
359: target := Expanded_Minors(planes(nd.level).all,map,expbp);
360: start := Expanded_Minors(special,map,expbp);
361: res(nd.level) := Linear_Interpolation(target,start,nvars);
362: if nd.tp = bottom
363: then
364: res(nd.level+1)
365: := Moving_Parameter(nvars,nvars-1,nvars,Create(1.0),s(nd.level));
366: else
367: res(nd.level+1)
368: := Moving_Parameter(nvars,nvars-1,nvars,Create(1.0),
369: (Create(1.0)/s(nd.level)));
370: Divide_Common_Factor(res(nd.level),nvars);
371: end if;
372: if nd.tp = top
373: then Standard_Complex_Poly_Matrices.Clear(map);
374: end if;
375: Clear(target); Clear(start);
376: return res;
377: end One_Quantum_Pieri_Homotopy;
378:
379: function Two_Quantum_Pieri_Homotopy
380: ( n : natural; nd : Node; expbp : Bracket_Polynomial;
381: xpm : Standard_Complex_Poly_Matrices.Matrix;
382: planes : VecMat; s : Vector ) return Poly_Sys is
383:
384: res : Poly_Sys(1..nd.level+2);
385: p : constant natural := nd.p;
386: m : constant natural := n-p;
387: nvars : constant natural := nd.level + p + 2;
388: top_special : constant Standard_Complex_Matrices.Matrix
389: := Special_Plane(m,Modulo(nd.top,m+p));
390: bot_special : constant Standard_Complex_Matrices.Matrix
391: := Special_Plane(m,Modulo(nd.bottom,m+p));
392: top_start,bot_start,top_target,bot_target : Poly;
393: map : Standard_Complex_Poly_Matrices.Matrix(xpm'range(1),xpm'range(2))
394: := Insert(xpm,nvars);
395:
396: begin
397: -- first do bottom pivots with s1 = nvars-1
398: bot_target := Expanded_Minors(planes(nd.level-1).all,map,expbp);
399: bot_start := Expanded_Minors(bot_special,map,expbp);
400: res(nd.level-1) := Linear_Interpolation(bot_target,bot_start,nvars+1);
401: Divide_Common_Factor(res(nd.level-1),nvars+1);
402: res(nd.level+1)
403: := Moving_Parameter(nvars+1,nvars-1,nvars+1,Create(1.0),s(nd.level-1));
404: Clear(bot_target); Clear(bot_start);
405: -- swap s1 with s2 to deal with the fixed equations
406: Swap(map,nvars-1,nvars);
407: for i in 1..nd.level-2 loop
408: declare
409: eva : Standard_Complex_Poly_Matrices.Matrix(map'range(1),map'range(2))
410: := Eval(map,s(i),Create(1.0));
411: begin
412: res(i) := Expanded_Minors(planes(i).all,eva,expbp);
413: Standard_Complex_Poly_Matrices.Clear(eva);
414: end;
415: end loop;
416: -- swap t with s2 for top pivots, s2 = nvars
417: Swap(map,nvars,nvars+1);
418: top_target := Expanded_Minors(planes(nd.level).all,map,expbp);
419: top_start := Expanded_Minors(top_special,map,expbp);
420: res(nd.level) := Linear_Interpolation(top_target,top_start,nvars+1);
421: res(nd.level+2)
422: := Moving_Parameter(nvars+1,nvars,nvars+1,Create(1.0),
423: (Create(1.0)/s(nd.level)));
424: Divide_Common_Factor(res(nd.level),nvars+1);
425: Clear(top_target); Clear(top_start);
426: Standard_Complex_Poly_Matrices.Clear(map);
427: return res;
428: end Two_Quantum_Pieri_Homotopy;
429:
430: function One_General_Quantum_Pieri_Homotopy
431: ( n,ind : natural; nd : Node; s_mode : natural;
432: bs : Bracket_System;
433: start,target : Standard_Complex_Matrices.Matrix;
434: xpm : Standard_Complex_Poly_Matrices.Matrix;
435: planes : VecMat; s : Vector ) return Link_to_Poly_Sys is
436:
437: res : Link_to_Poly_Sys;
438: p : constant natural := nd.p;
439: m : constant natural := n-p;
440: nvars : constant natural := nd.level + p + 2;
441: -- nvars = level #equations in the x_ij's
442: -- + p because not yet fixed the ones
443: -- + 2 for s and t, note that t is continuation parameter
444: eva,map : Standard_Complex_Poly_Matrices.Matrix(1..n,1..p);
445: moving_plane : Standard_Complex_Poly_Matrices.Matrix(1..n,target'range(2))
446: := Moving_U_Matrix(nvars,start,target);
447: moving : Poly_Sys(1..bs'last);
448: movpar : Poly_Sys(1..1);
449: inddiv : natural;
450: starcol : Vector(start'range(1));
451: deg1 : constant natural := Degree1(nd,n);
452:
453: begin
454: if s_mode = 0
455: then for i in starcol'range loop
456: starcol(i) := start(i,1);
457: end loop;
458: end if;
459: -- PART I : fixed equations
460: case nd.tp is
461: when top => Standard_Complex_Poly_Matrices.Copy(xpm,map);
462: Swap(map,nvars-1,nvars);
463: when bottom => map := xpm;
464: when others => null; -- mixed case treated separately
465: end case;
466: for i in 1..ind-1 loop
467: eva := Eval(xpm,s(i),Create(1.0));
468: if s_mode = 0
469: then --Multiply(eva,1,nvars-1,deg1);
470: --Add(eva,1,nvars-1,starcol);
471: Multiply(eva,1,nvars-1,deg1);
472: Add(eva,1,nvars,starcol);
473: end if;
474: Concat(res,Polynomial_Equations(planes(i).all,eva));
475: Standard_Complex_Poly_Matrices.Clear(eva);
476: end loop;
477: if res = null
478: then inddiv := 0;
479: else inddiv := res'last;
480: end if;
481: -- PART II : moving equation for the plane
482: if s_mode = 2
483: then moving := Expanded_Minors(moving_plane,map,bs);
484: else eva := Eval(xpm,Create(1.0),Create(0.0));
485: if s_mode = 0
486: then --Multiply(eva,1,nvars-1,deg1);
487: --Add(eva,1,nvars-1,starcol);
488: Multiply(eva,1,nvars-1,deg1);
489: Add(eva,1,nvars,starcol);
490: end if;
491: moving := Expanded_Minors(moving_plane,eva,bs);
492: Standard_Complex_Poly_Matrices.Clear(eva);
493: end if;
494: Concat(res,moving);
495: if nd.tp = top
496: then Standard_Complex_Poly_Matrices.Clear(map);
497: for i in inddiv+1..res'last loop
498: Divide_Common_Factor(res(i),nvars);
499: end loop;
500: end if;
501: -- PART III : moving equation for the interpolation points
502: case s_mode is
503: when 0 => -- move s from 0 to 1
504: case nd.tp is
505: when bottom =>
506: movpar(1)
507: := Moving_Parameter(nvars,nvars-1,nvars,
508: Create(0.0),Create(1.0));
509: when others => null;
510: end case;
511: when 1 => -- leave s constant at 1
512: movpar(1) := Constant_Parameter(nvars,nvars-1,Create(1.0));
513: when 2 => -- move s from 1 to target
514: case nd.tp is
515: when top =>
516: movpar(1) -- move s from 1 to 1/s_ind (swap of s and t)
517: := Moving_Parameter(nvars,nvars-1,nvars,Create(1.0),
518: (Create(1.0)/s(ind)));
519: when bottom =>
520: movpar(1) -- move s from 1 to s_ind
521: := Moving_Parameter(nvars,nvars-1,nvars,Create(1.0),s(ind));
522: when others => null;
523: end case;
524: when others => null;
525: end case;
526: Concat(res,movpar);
527: return res;
528: end One_General_Quantum_Pieri_Homotopy;
529:
530: function Two_General_Quantum_Pieri_Homotopy
531: ( n,ind : natural; nd : Node; top_bs,bot_bs : Bracket_System;
532: top_start,top_target,bot_start,bot_target
533: : Standard_Complex_Matrices.Matrix;
534: xpm : Standard_Complex_Poly_Matrices.Matrix;
535: planes : VecMat; s : Vector ) return Link_to_Poly_Sys is
536:
537: res : Link_to_Poly_Sys;
538:
539: begin
540: return res;
541: end Two_General_Quantum_Pieri_Homotopy;
542:
543: end Pieri_Homotopies;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to pieri_homotopies.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