Annotation of OpenXM_contrib/PHC/Ada/Schubert/sagbi_homotopies.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; use Standard_Natural_Vectors;
4: with Standard_Complex_Poly_Functions; use Standard_Complex_Poly_Functions;
5: with Brackets; use Brackets;
6: with Bracket_Monomials; use Bracket_Monomials;
7: with Bracket_Polynomials; use Bracket_Polynomials;
8: with Straightening_Syzygies; use Straightening_Syzygies;
9: with Bracket_Expansions; use Bracket_Expansions;
10: with Evaluated_Minors; use Evaluated_Minors;
11:
12: package body SAGBI_Homotopies is
13:
14: function Coordinatize_Hexadecimal ( b : Bracket ) return natural is
15:
16: -- DESCRIPTION :
17: -- Returns the hexadecimal expansion of the entries in the bracket.
18:
19: res : natural := 0;
20:
21: begin
22: for i in b'range loop
23: res := res*16 + b(i);
24: end loop;
25: return res;
26: end Coordinatize_Hexadecimal;
27:
28: function Unsigned ( i : integer ) return natural is
29:
30: -- DESCRIPTION :
31: -- Returns the unsigned integer.
32:
33: n : natural;
34:
35: begin
36: if i < 0
37: then n := -i;
38: else n := i;
39: end if;
40: return n;
41: end Unsigned;
42:
43: function Bracketize_Hexadecimal ( n,d : natural ) return Bracket is
44:
45: -- DESCRIPTION :
46: -- Returns the d-bracket from the hexadecimal expansion n.
47:
48: res : Bracket(1..d);
49: nn : natural := n;
50:
51: begin
52: for i in reverse 1..d loop
53: res(i) := nn mod 16;
54: nn := nn/16;
55: end loop;
56: return res;
57: end Bracketize_Hexadecimal;
58:
59: function Coordinatize ( p : Bracket_Polynomial ) return Bracket_Polynomial is
60:
61: -- DESCRIPTION :
62: -- Replaces the first bracket in every monomial by the decimal expansion.
63:
64: res : Bracket_Polynomial;
65:
66: procedure Coordinatize_Term ( t : in Bracket_Term; cont1 : out boolean ) is
67:
68: first,second : boolean;
69: bm : Bracket_Monomial;
70: bt : Bracket_Term;
71:
72: procedure Coordinatize_Bracket ( b : in Bracket; cont2 : out boolean ) is
73: begin
74: if first
75: then bt.coeff := Create(double_float(Coordinatize_Hexadecimal(b)));
76: first := false;
77: second := true;
78: elsif second
79: then bm := Create(b);
80: else Multiply(bm,b);
81: end if;
82: cont2 := true;
83: end Coordinatize_Bracket;
84: procedure Coordinatize_Brackets is
85: new Enumerate_Brackets(Coordinatize_Bracket);
86:
87: begin
88: first := true; second := false;
89: Coordinatize_Brackets(t.monom);
90: bt.monom := bm;
91: if REAL_PART(t.coeff) < 0.0
92: then Min(res,bt);
93: else Add(res,bt);
94: end if;
95: cont1 := true;
96: end Coordinatize_Term;
97: procedure Coordinatize_Terms is new Enumerate_Terms(Coordinatize_Term);
98:
99: begin
100: Coordinatize_Terms(p);
101: return res;
102: end Coordinatize;
103:
104: procedure Divide ( p : in out Poly; w : in natural ) is
105:
106: -- DESCRIPTION :
107: -- Divides the polynomial by t^w.
108:
109: procedure Divide_Term ( t : in out Term; continue : out boolean ) is
110: begin
111: t.dg(t.dg'last) := t.dg(t.dg'last) - w;
112: continue := true;
113: end Divide_Term;
114: procedure Divide_Terms is new Changing_Iterator(Divide_Term);
115:
116: begin
117: Divide_Terms(p);
118: end Divide;
119:
120: function Weight ( e : Standard_Natural_Vectors.Vector; n,d : natural )
121: return natural is
122:
123: -- DESCRIPTION :
124: -- Returns the weight of the exponent vector for the localization that
125: -- takes the d-by-d identitity matrix in the lower-right of the d-plane.
126: -- The lifting recipe is xij*t^((i-1)*(d-j)).
127:
128: res : natural := 0;
129: jmp,ind : natural;
130:
131: begin
132: for j in 1..d loop
133: jmp := d-j;
134: for i in 1..n-d loop
135: ind := (i-1)*d + j;
136: if e(ind) > 0
137: then res := res + (i-1)*jmp;
138: end if;
139: end loop;
140: end loop;
141: return res;
142: end Weight;
143:
144: function Weight ( locmap : Standard_Natural_Matrices.Matrix;
145: e : Standard_Natural_Vectors.Vector ) return natural is
146:
147: -- DESCRIPTION :
148: -- Returns the weight of the exponent vector as xij*t^((i-1)*(d-j))
149: -- for the localization pattern in locmap.
150:
151: res : natural := 0;
152: d : constant natural := locmap'length(2);
153: jmp,ind : natural;
154:
155: begin
156: ind := 0;
157: for i in locmap'range(1) loop
158: for j in locmap'range(2) loop
159: jmp := d-j;
160: if locmap(i,j) = 2
161: then ind := ind+1;
162: if e(ind) > 0
163: then res := res + (i-1)*jmp;
164: end if;
165: end if;
166: end loop;
167: end loop;
168: return res;
169: end Weight;
170:
171: function Lift ( p : Poly; n,d : natural ) return Poly is
172:
173: -- DESCRIPTION :
174: -- Returns the lifted polynomial, where the xij is lifted according
175: -- to xij*t^((i-1)*(d-j)). The lowest powers of t are divided out.
176: -- The d-by-d identity matrix is the lower-right of the d-plane.
177:
178: res : Poly := Null_Poly;
179: first : boolean := true;
180: minwei : natural;
181:
182: procedure Lift_Term ( t : in Term; continue : out boolean ) is
183:
184: tt : Term;
185: wei : natural;
186:
187: begin
188: tt.cf := t.cf;
189: tt.dg := new Standard_Natural_Vectors.Vector(1..t.dg'last+1);
190: tt.dg(t.dg'range) := t.dg.all;
191: wei := Weight(t.dg.all,n,d);
192: tt.dg(tt.dg'last) := wei;
193: Add(res,tt);
194: Clear(tt.dg);
195: if first
196: then minwei := wei;
197: first := false;
198: elsif wei < minwei
199: then minwei := wei;
200: end if;
201: continue := true;
202: end Lift_Term;
203: procedure Lift_Terms is new Visiting_Iterator(Lift_Term);
204:
205: begin
206: Lift_Terms(p);
207: if minwei /= 0
208: then Divide(res,minwei);
209: end if;
210: return res;
211: end Lift;
212:
213: function Lift ( locmap : Standard_Natural_Matrices.Matrix; p : Poly )
214: return Poly is
215:
216: -- DESCRIPTION :
217: -- Lifts p as to xij*t^((i-1)*(d-j)) and divides by the lowest powers
218: -- of t, respecting the localization pattern in locmap.
219:
220: res : Poly := Null_Poly;
221: first : boolean := true;
222: minwei : natural;
223:
224: procedure Lift_Term ( t : in Term; continue : out boolean ) is
225:
226: tt : Term;
227: wei : natural;
228:
229: begin
230: tt.cf := t.cf;
231: tt.dg := new Standard_Natural_Vectors.Vector(1..t.dg'last+1);
232: tt.dg(t.dg'range) := t.dg.all;
233: wei := Weight(locmap,t.dg.all);
234: tt.dg(tt.dg'last) := wei;
235: Add(res,tt);
236: Clear(tt.dg);
237: if first
238: then minwei := wei;
239: first := false;
240: elsif wei < minwei
241: then minwei := wei;
242: end if;
243: continue := true;
244: end Lift_Term;
245: procedure Lift_Terms is new Visiting_Iterator(Lift_Term);
246:
247: begin
248: Lift_Terms(p);
249: if minwei /= 0
250: then Divide(res,minwei);
251: end if;
252: return res;
253: end Lift;
254:
255: -- TARGET ROUTINES :
256:
257: function Lifted_Localized_Laplace_Expansion ( n,d : natural ) return Poly is
258:
259: res : Poly := Null_Poly;
260: p : Bracket_Polynomial := Laplace_Expansion(n,n-d);
261:
262: procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
263:
264: first : boolean := true;
265: cf : integer;
266:
267: procedure Visit_Bracket ( b : in Bracket; cont : out boolean ) is
268:
269: pb,lp : Poly;
270:
271: begin
272: if first
273: then cf := Coordinatize_Hexadecimal(b);
274: first := false;
275: else pb := Localized_Expand(n,d,b);
276: lp := Lift(pb,n,d); Clear(pb);
277: Mul(lp,Create(double_float(cf)));
278: Add(res,lp);
279: Clear(lp);
280: end if;
281: cont := true;
282: end Visit_Bracket;
283: procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
284:
285: begin
286: Visit_Brackets(t.monom);
287: continue := true;
288: end Visit_Term;
289: procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
290:
291: begin
292: Visit_Terms(p);
293: return res;
294: end Lifted_Localized_Laplace_Expansion;
295:
296: function Lifted_Localized_Laplace_Expansion
297: ( locmap : Standard_Natural_Matrices.Matrix ) return Poly is
298:
299: res : Poly := Null_Poly;
300: n : constant natural := locmap'length(1);
301: d : constant natural := locmap'length(2);
302: p : Bracket_Polynomial := Laplace_Expansion(n,n-d);
303:
304: procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
305:
306: first : boolean := true;
307: cf : integer;
308:
309: procedure Visit_Bracket ( b : in Bracket; cont : out boolean ) is
310:
311: pb,lp : Poly;
312:
313: begin
314: if first
315: then cf := Coordinatize_Hexadecimal(b);
316: first := false;
317: else pb := Expand(locmap,b);
318: Reduce_Variables(locmap,pb);
319: lp := Lift(locmap,pb); Clear(pb);
320: Mul(lp,Create(double_float(cf)));
321: Add(res,lp);
322: Clear(lp);
323: end if;
324: cont := true;
325: end Visit_Bracket;
326: procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
327:
328: begin
329: Visit_Brackets(t.monom);
330: continue := true;
331: end Visit_Term;
332: procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
333:
334: begin
335: Visit_Terms(p);
336: return res;
337: end Lifted_Localized_Laplace_Expansion;
338:
339: function Intersection_Coefficients
340: ( m : Standard_Floating_Matrices.Matrix;
341: c : Standard_Complex_Vectors.Vector )
342: return Standard_Complex_Vectors.Vector is
343:
344: res : Standard_Complex_Vectors.Vector(c'range);
345: nmd : constant natural := m'last(2);
346: ind : integer;
347: b : Bracket(1..nmd);
348:
349: begin
350: for i in c'range loop
351: ind := integer(REAL_PART(c(i)));
352: b := Bracketize_Hexadecimal(Unsigned(ind),nmd);
353: if ind > 0
354: then res(i) := Create(Determinant(m,b));
355: else res(i) := Create(-Determinant(m,b));
356: end if;
357: end loop;
358: return res;
359: end Intersection_Coefficients;
360:
361: function Intersection_Coefficients
362: ( m : Standard_Complex_Matrices.Matrix;
363: c : Standard_Complex_Vectors.Vector )
364: return Standard_Complex_Vectors.Vector is
365:
366: res : Standard_Complex_Vectors.Vector(c'range);
367: nmd : constant natural := m'last(2);
368: ind : integer;
369: b : Bracket(1..nmd);
370:
371: begin
372: for i in c'range loop
373: ind := integer(REAL_PART(c(i)));
374: b := Bracketize_Hexadecimal(Unsigned(ind),nmd);
375: if ind > 0
376: then res(i) := Determinant(m,b);
377: else res(i) := -Determinant(m,b);
378: end if;
379: end loop;
380: return res;
381: end Intersection_Coefficients;
382:
383: function Intersection_Condition
384: ( m : Standard_Floating_Matrices.Matrix; p : Poly ) return Poly is
385:
386: res : Poly := Null_Poly;
387: nmd : constant natural := m'last(2);
388:
389: procedure Visit_Term ( t : in Term; continue : out boolean ) is
390:
391: c : integer := integer(REAL_PART(t.cf));
392: b : Bracket(1..nmd) := Bracketize_Hexadecimal(Unsigned(c),nmd);
393: det : double_float := Determinant(m,b);
394: rt : Term;
395:
396: begin
397: if c > 0
398: then rt.cf := Create(det);
399: else rt.cf := Create(-det);
400: end if;
401: rt.dg := t.dg;
402: Add(res,rt);
403: continue := true;
404: end Visit_Term;
405: procedure Visit_Terms is new Visiting_Iterator(Visit_Term);
406:
407: begin
408: Visit_Terms(p);
409: return res;
410: end Intersection_Condition;
411:
412: function Intersection_Condition
413: ( m : Standard_Complex_Matrices.Matrix; p : Poly ) return Poly is
414:
415: res : Poly := Null_Poly;
416: nmd : constant natural := m'last(2);
417:
418: procedure Visit_Term ( t : in Term; continue : out boolean ) is
419:
420: c : integer := integer(REAL_PART(t.cf));
421: b : Bracket(1..nmd) := Bracketize_Hexadecimal(Unsigned(c),nmd);
422: det : Complex_Number := Determinant(m,b);
423: rt : Term;
424:
425: begin
426: if c > 0
427: then rt.cf := det;
428: else rt.cf := -det;
429: end if;
430: rt.dg := t.dg;
431: Add(res,rt);
432: continue := true;
433: end Visit_Term;
434: procedure Visit_Terms is new Visiting_Iterator(Visit_Term);
435:
436: begin
437: Visit_Terms(p);
438: return res;
439: end Intersection_Condition;
440:
441: end SAGBI_Homotopies;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to sagbi_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