Annotation of OpenXM_contrib/PHC/Ada/Schubert/numeric_minor_equations.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_Poly_Functions; use Standard_Complex_Poly_Functions;
4: with Standard_Complex_Poly_SysFun; use Standard_Complex_Poly_SysFun;
5: with Brackets; use Brackets;
6: with Plane_Representations; use Plane_Representations;
7: with Evaluated_Minors; use Evaluated_Minors;
8: with Symbolic_Minor_Equations; use Symbolic_Minor_Equations;
9:
10: package body Numeric_Minor_Equations is
11:
12: tol : constant double_float := 10.0**(-10);
13:
14: -- EXPANDING ACCORDING A BRACKET MONOMIAL :
15:
16: function Expanded_Minors
17: ( cffmat : Standard_Floating_Matrices.Matrix;
18: polmat : Standard_Complex_Poly_Matrices.Matrix;
19: bm : Bracket_Monomial ) return Poly is
20:
21: res : Poly := Null_Poly;
22: first : boolean := true;
23:
24: procedure Visit_Bracket ( b : in Bracket; continue : out boolean ) is
25:
26: factor : double_float;
27: minor : Poly;
28:
29: begin
30: if first
31: then declare
32: bb : Bracket(b'first..b'last-1) := b(b'first+1..b'last);
33: begin
34: factor := Determinant(cffmat,bb);
35: end;
36: first := false;
37: else minor := Expanded_Minor(polmat,b);
38: if (minor /= Null_Poly) and (abs(factor) > tol)
39: then Mul(minor,Create(factor));
40: Add(res,minor);
41: end if;
42: Clear(minor);
43: end if;
44: continue := true;
45: end Visit_Bracket;
46: procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
47:
48: begin
49: Visit_Brackets(bm);
50: return res;
51: end Expanded_Minors;
52:
53: function Expanded_Minors
54: ( cffmat : Standard_Complex_Matrices.Matrix;
55: polmat : Standard_Complex_Poly_Matrices.Matrix;
56: bm : Bracket_Monomial ) return Poly is
57:
58: res : Poly := Null_Poly;
59: first : boolean := true;
60: factor : Complex_Number;
61:
62: procedure Visit_Bracket ( b : in Bracket; continue : out boolean ) is
63:
64: minor : Poly;
65:
66: begin
67: if first
68: then declare
69: bb : Bracket(b'first..b'last-1) := b(b'first+1..b'last);
70: begin
71: factor := Determinant(cffmat,bb);
72: end;
73: first := false;
74: else minor := Expanded_Minor(polmat,b);
75: if (minor /= Null_Poly) and (AbsVal(factor) > tol)
76: then Mul(minor,factor);
77: Add(res,minor);
78: end if;
79: Clear(minor);
80: end if;
81: continue := true;
82: end Visit_Bracket;
83: procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
84:
85: begin
86: Visit_Brackets(bm);
87: return res;
88: end Expanded_Minors;
89:
90: function Expanded_Minors
91: ( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
92: bm : Bracket_Monomial ) return Poly is
93:
94: res : Poly := Null_Poly;
95: first : boolean := true;
96: factor : Poly;
97:
98: procedure Visit_Bracket ( b : in Bracket; continue : out boolean ) is
99:
100: minor : Poly;
101:
102: begin
103: if first
104: then declare
105: bb : Bracket(b'first..b'last-1) := b(b'first+1..b'last);
106: begin
107: factor := Expanded_Minor(cntmat,bb);
108: end;
109: first := false;
110: else minor := Expanded_Minor(polmat,b);
111: if (minor /= Null_Poly) and (factor /= Null_Poly)
112: then Mul(minor,factor);
113: Add(res,minor);
114: end if;
115: Clear(factor); Clear(minor);
116: end if;
117: continue := true;
118: end Visit_Bracket;
119: procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
120:
121: begin
122: Visit_Brackets(bm);
123: return res;
124: end Expanded_Minors;
125:
126: function Lifted_Expanded_Minors
127: ( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
128: bm : Bracket_Monomial ) return Poly is
129:
130: res : Poly := Null_Poly;
131: first : boolean := true;
132: factor : Poly;
133:
134: procedure Visit_Bracket ( b : in Bracket; continue : out boolean ) is
135:
136: minor,extmin : Poly;
137:
138: begin
139: if first
140: then declare
141: bb : Bracket(b'first..b'last-1) := b(b'first+1..b'last);
142: begin
143: factor := Expanded_Minor(cntmat,bb);
144: end;
145: first := false;
146: else minor := Expanded_Minor(polmat,b);
147: if (minor /= Null_Poly) and (factor /= Null_Poly)
148: then extmin := Extend_Zero_Lifting(minor);
149: Mul(extmin,factor);
150: Add(res,extmin);
151: end if;
152: Clear(factor); Clear(minor); Clear(extmin);
153: end if;
154: continue := true;
155: end Visit_Bracket;
156: procedure Visit_Brackets is new Enumerate_Brackets(Visit_Bracket);
157:
158: begin
159: Visit_Brackets(bm);
160: return res;
161: end Lifted_Expanded_Minors;
162:
163: -- EXPANDING ACCORDING A BRACKET POLYNOMIAL :
164:
165: function Expanded_Minors
166: ( cffmat : Standard_Floating_Matrices.Matrix;
167: polmat : Standard_Complex_Poly_Matrices.Matrix;
168: bp : Bracket_Polynomial ) return Poly is
169:
170: res : Poly := Null_Poly;
171:
172: procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
173:
174: minor : Poly := Expanded_Minors(cffmat,polmat,t.monom);
175:
176: begin
177: if REAL_PART(t.coeff) > 0.0
178: then Add(res,minor);
179: else Sub(res,minor);
180: end if;
181: Clear(minor);
182: continue := true;
183: end Visit_Term;
184: procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
185:
186: begin
187: Visit_Terms(bp);
188: return res;
189: end Expanded_Minors;
190:
191: function Expanded_Minors
192: ( cffmat : Standard_Complex_Matrices.Matrix;
193: polmat : Standard_Complex_Poly_Matrices.Matrix;
194: bp : Bracket_Polynomial ) return Poly is
195:
196: res : Poly := Null_Poly;
197:
198: procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
199:
200: minor : Poly := Expanded_Minors(cffmat,polmat,t.monom);
201:
202: begin
203: if REAL_PART(t.coeff) > 0.0
204: then Add(res,minor);
205: else Sub(res,minor);
206: end if;
207: Clear(minor);
208: continue := true;
209: end Visit_Term;
210: procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
211:
212: begin
213: Visit_Terms(bp);
214: return res;
215: end Expanded_Minors;
216:
217: function Expanded_Minors
218: ( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
219: bp : Bracket_Polynomial ) return Poly is
220:
221: res : Poly := Null_Poly;
222:
223: procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
224:
225: minor : Poly := Expanded_Minors(cntmat,polmat,t.monom);
226:
227: begin
228: if REAL_PART(t.coeff) > 0.0
229: then Add(res,minor);
230: else Sub(res,minor);
231: end if;
232: Clear(minor);
233: continue := true;
234: end Visit_Term;
235: procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
236:
237: begin
238: Visit_Terms(bp);
239: return res;
240: end Expanded_Minors;
241:
242: function Lifted_Expanded_Minors
243: ( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
244: bp : Bracket_Polynomial ) return Poly is
245:
246: res : Poly := Null_Poly;
247:
248: procedure Visit_Term ( t : in Bracket_Term; continue : out boolean ) is
249:
250: minor : Poly := Lifted_Expanded_Minors(cntmat,polmat,t.monom);
251:
252: begin
253: if REAL_PART(t.coeff) > 0.0
254: then Add(res,minor);
255: else Sub(res,minor);
256: end if;
257: Clear(minor);
258: continue := true;
259: end Visit_Term;
260: procedure Visit_Terms is new Enumerate_Terms(Visit_Term);
261:
262: begin
263: Visit_Terms(bp);
264: return res;
265: end Lifted_Expanded_Minors;
266:
267: -- EXPANDING TO CONSTRUCT POLYNOMIAL SYSTEMS :
268:
269: function Expanded_Minors
270: ( cffmat : Standard_Floating_Matrices.Matrix;
271: polmat : Standard_Complex_Poly_Matrices.Matrix;
272: bs : Bracket_System ) return Poly_Sys is
273:
274: res : Poly_Sys(bs'first+1..bs'last);
275:
276: begin
277: for i in res'range loop
278: res(i) := Expanded_Minors(cffmat,polmat,bs(i));
279: end loop;
280: return res;
281: end Expanded_Minors;
282:
283: function Expanded_Minors
284: ( cffmat : Standard_Complex_Matrices.Matrix;
285: polmat : Standard_Complex_Poly_Matrices.Matrix;
286: bs : Bracket_System ) return Poly_Sys is
287:
288: res : Poly_Sys(bs'first+1..bs'last);
289:
290: begin
291: for i in res'range loop
292: res(i) := Expanded_Minors(cffmat,polmat,bs(i));
293: end loop;
294: return res;
295: end Expanded_Minors;
296:
297: function Expanded_Minors
298: ( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
299: bs : Bracket_System ) return Poly_Sys is
300:
301: res : Poly_Sys(bs'first+1..bs'last);
302:
303: begin
304: for i in res'range loop
305: res(i) := Expanded_Minors(cntmat,polmat,bs(i));
306: end loop;
307: return res;
308: end Expanded_Minors;
309:
310: function Lifted_Expanded_Minors
311: ( cntmat,polmat : Standard_Complex_Poly_Matrices.Matrix;
312: bs : Bracket_System ) return Poly_Sys is
313:
314: res : Poly_Sys(bs'first+1..bs'last);
315:
316: begin
317: for i in res'range loop
318: res(i) := Lifted_Expanded_Minors(cntmat,polmat,bs(i));
319: end loop;
320: return res;
321: end Lifted_Expanded_Minors;
322:
323: function Evaluate ( p : Poly; x : Standard_Complex_Matrices.Matrix )
324: return Complex_Number is
325:
326: xv : constant Standard_Complex_Vectors.Vector := Vector_Rep(x);
327:
328: begin
329: return Eval(p,xv);
330: end Evaluate;
331:
332: function Evaluate ( p : Poly_Sys; x : Standard_Complex_Matrices.Matrix )
333: return Standard_Complex_Vectors.Vector is
334:
335: xv : constant Standard_Complex_Vectors.Vector := Vector_Rep(x);
336:
337: begin
338: return Eval(p,xv);
339: end Evaluate;
340:
341: procedure Embed ( t : in out Term ) is
342:
343: dg : Degrees
344: := new Standard_Natural_Vectors.Vector(t.dg'first..t.dg'last+1);
345:
346: begin
347: dg(t.dg'range) := t.dg.all;
348: dg(dg'last) := 0;
349: Clear(t.dg);
350: t.dg := dg;
351: end Embed;
352:
353: function Embed ( t : Term ) return Term is
354:
355: res : Term;
356:
357: begin
358: res.dg := new Standard_Natural_Vectors.Vector(t.dg'first..t.dg'last+1);
359: res.dg(t.dg'range) := t.dg.all;
360: res.dg(res.dg'last) := 0;
361: res.cf := t.cf;
362: return res;
363: end Embed;
364:
365: procedure Embed ( p : in out Poly ) is
366:
367: res : Poly := Null_Poly;
368:
369: procedure Embed_Term ( t : in Term; continue : out boolean ) is
370:
371: et : Term := Embed(t);
372:
373: begin
374: Add(res,et);
375: Clear(et);
376: continue := true;
377: end Embed_Term;
378: procedure Embed_Terms is new Visiting_Iterator(Embed_Term);
379:
380: begin
381: Embed_Terms(p);
382: Clear(p);
383: p := res;
384: end Embed;
385:
386: procedure Embed ( p : in out Poly_Sys ) is
387: begin
388: for i in p'range loop
389: Embed(p(i));
390: end loop;
391: end Embed;
392:
393: procedure Embed ( m : in out Standard_Complex_Poly_Matrices.Matrix ) is
394: begin
395: for i in m'range(1) loop
396: for j in m'range(2) loop
397: if m(i,j) /= Null_Poly
398: then Embed(m(i,j));
399: end if;
400: end loop;
401: end loop;
402: end Embed;
403:
404: function Linear_Homotopy ( target,start : Poly ) return Poly is
405:
406: res : Poly := Null_Poly;
407:
408: procedure Embed_Target_Term ( t : in Term; continue : out boolean ) is
409:
410: et : Term := Embed(t);
411:
412: begin
413: et.dg(et.dg'last) := 1;
414: Add(res,et);
415: Clear(et);
416: continue := true;
417: end Embed_Target_Term;
418: procedure Embed_Target_Terms is new Visiting_Iterator(Embed_Target_Term);
419:
420: procedure Embed_Start_Term ( t : in Term; continue : out boolean ) is
421:
422: et : Term := Embed(t);
423:
424: begin
425: Add(res,et);
426: et.dg(et.dg'last) := 1;
427: Sub(res,et);
428: Clear(et);
429: continue := true;
430: end Embed_Start_Term;
431: procedure Embed_Start_Terms is new Visiting_Iterator(Embed_Start_Term);
432:
433: begin
434: Embed_Target_Terms(target);
435: Embed_Start_Terms(start);
436: return res;
437: end Linear_Homotopy;
438:
439: function Linear_Interpolation
440: ( target,start : Poly; k : natural ) return Poly is
441:
442: res : Poly := Null_Poly;
443:
444: procedure Embed_Target_Term ( t : in Term; continue : out boolean ) is
445:
446: et : Term;
447:
448: begin
449: Copy(t,et);
450: et.dg(k) := et.dg(k) + 1; -- multiply with t
451: Add(res,et);
452: Clear(et);
453: continue := true;
454: end Embed_Target_Term;
455: procedure Embed_Target_Terms is new Visiting_Iterator(Embed_Target_Term);
456:
457: procedure Embed_Start_Term ( t : in Term; continue : out boolean ) is
458:
459: et : Term;
460:
461: begin
462: Copy(t,et);
463: Add(res,et); -- res := res + et
464: et.dg(k) := et.dg(k) + 1; -- multiply with t
465: Sub(res,et); -- res := res + et - t*et
466: Clear(et);
467: continue := true;
468: end Embed_Start_Term;
469: procedure Embed_Start_Terms is new Visiting_Iterator(Embed_Start_Term);
470:
471: begin
472: Embed_Target_Terms(target);
473: Embed_Start_Terms(start);
474: return res;
475: end Linear_Interpolation;
476:
477: procedure Divide_Common_Factor ( p : in out Poly; k : in natural ) is
478:
479: first : boolean := true;
480: min : natural;
481:
482: procedure Min_Power ( t : in Term; continue : out boolean ) is
483: begin
484: if first
485: then first := false;
486: min := t.dg(k); -- initialize minimal power
487: else if t.dg(k) < min
488: then min := t.dg(k);
489: end if;
490: end if;
491: continue := true;
492: end Min_Power;
493: procedure Find_Min_Power is new Visiting_Iterator(Min_Power);
494:
495: procedure Divide ( t : in out Term; continue : out boolean ) is
496: begin
497: t.dg(k) := t.dg(k) - min;
498: continue := true;
499: end Divide;
500: procedure Divide_Min_Power is new Changing_Iterator(Divide);
501:
502: begin
503: Find_Min_Power(p);
504: if min > 0
505: then Divide_Min_Power(p);
506: end if;
507: end Divide_Common_Factor;
508:
509: end Numeric_Minor_Equations;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to numeric_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