Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Product/interpolating_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;
4: with Standard_Complex_Vectors;
5: with Standard_Complex_Matrices;
6: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
7: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
8: with Standard_Complex_Poly_Randomizers; use Standard_Complex_Poly_Randomizers;
9: with Sets_of_Unknowns; use Sets_of_Unknowns;
10: with Degrees_in_Sets_of_Unknowns; use Degrees_in_Sets_of_Unknowns;
11:
12: package body Interpolating_Homotopies is
13:
14: -- AUXILIARY OPERATIONS :
15:
16: function Initial_Term ( p : Poly ) return Term is
17:
18: -- DESCRIPTION :
19: -- Returns the initial term of p.
20:
21: res : Term;
22:
23: procedure Init_Term ( t : Term; continue : out boolean ) is
24: begin
25: res := t;
26: continue := false;
27: end Init_Term;
28: procedure Init_Terms is new Visiting_Iterator(Init_Term);
29:
30: begin
31: Init_Terms(p);
32: return res;
33: end Initial_Term;
34:
35: procedure Eliminate_Term ( p : in out Poly; dg : in Degrees ) is
36:
37: -- DESCRIPTION :
38: -- Eliminates the monomial in p that has the given exponent vector.
39:
40: c : Complex_Number := Coeff(p,dg);
41:
42: begin
43: if c /= Create(0.0)
44: then declare
45: t : Term;
46: begin
47: t.cf := c; t.dg := dg;
48: Sub(p,t);
49: end;
50: end if;
51: end Eliminate_Term;
52:
53: function Admitted ( t : Term; i : natural; z : partition;
54: d : Standard_Integer_Matrices.Matrix ) return boolean is
55:
56: -- DESCRIPTION :
57: -- Returns true if the given term does not violate the m-homogeneous
58: -- structure, i.e., if Degree(t,z(j)) <= d(i,j), for j in z'range.
59:
60: begin
61: for j in z'range loop
62: if Degree(t,z(j)) > d(i,j)
63: then return false;
64: end if;
65: end loop;
66: return true;
67: end Admitted;
68:
69: function Dense_Representation ( p : Poly; i : natural; z : Partition;
70: d : Matrix ) return Poly is
71:
72: -- DESCRIPTION :
73: -- Returns the dense representation of the given polynomial p, known as
74: -- p(i) of some polynomial system, of the given m-homogeneous structure.
75: -- The returned polynomial has all its coefficients equal to one.
76:
77: res : Poly := Null_Poly;
78: maxdegs,accu : Standard_Natural_Vectors.Vector(d'range(1));
79:
80: procedure Generate_Monomials
81: ( k : in natural; max : in Standard_Natural_Vectors.Vector;
82: acc : in out Standard_Natural_Vectors.Vector) is
83:
84: -- DESCRIPTION :
85: -- Generates all monomials with exponents in a box, with upper bounds
86: -- given by the vector max. Every monomial that respects the
87: -- m-homogeneous structure will be added to the result res.
88: -- The current unknown is indicated by the parameter k.
89:
90: t : Term;
91:
92: begin
93: for j in 0..max(k) loop
94: acc(k) := j;
95: t.cf := Create(1.0);
96: t.dg := new Standard_Natural_Vectors.Vector'(acc);
97: if Admitted(t,i,z,d)
98: then Add(res,t);
99: end if;
100: Clear(t);
101: if k < max'last
102: then Generate_Monomials(k+1,max,acc);
103: end if;
104: end loop;
105: end Generate_Monomials;
106:
107: begin
108: for j in maxdegs'range loop
109: maxdegs(j) := Degree(p,j); accu(j) := 0;
110: end loop;
111: Generate_Monomials(1,maxdegs,accu);
112: return res;
113: end Dense_Representation;
114:
115: function Evaluate ( t : Term; x : Standard_Complex_Vectors.Vector )
116: return Complex_Number is
117:
118: res : Complex_Number := Create(1.0);
119:
120: begin
121: for i in x'range loop
122: for k in 1..t.dg(i) loop
123: res := res*x(i);
124: end loop;
125: end loop;
126: return res;
127: end Evaluate;
128:
129: procedure Interpolation_System
130: ( p : in Poly; b : in natural; sols : in Solution_List;
131: mat : out Standard_Complex_Matrices.Matrix;
132: rhs : out Standard_Complex_Vectors.Vector ) is
133:
134: -- DESCRIPTION :
135: -- Returns the matrix and right hand side vector of the linear system
136: -- that expresses the interpolationg conditions.
137:
138: m : Standard_Complex_Matrices.Matrix(1..b,1..b);
139: v : Standard_Complex_Vectors.Vector(1..b);
140: cnt : natural := 0;
141:
142: procedure Scan_Term ( t : in Term; continue : out boolean ) is
143:
144: tmp : Solution_List := sols;
145: s : Link_to_Solution;
146:
147: begin
148: for row in m'range(1) loop -- fill in column indicated by cnt
149: s := Head_Of(tmp);
150: if cnt = 0
151: then v(row) := -t.cf*Evaluate(t,s.v);
152: elsif cnt <= b
153: then m(row,cnt) := Evaluate(t,s.v);
154: else v(row) := v(row) - t.cf*Evaluate(t,s.v);
155: end if;
156: tmp := Tail_Of(tmp);
157: end loop;
158: cnt := cnt + 1;
159: continue := true;
160: end Scan_Term;
161: procedure Scan_Poly is new Visiting_Iterator(Scan_Term);
162:
163: begin
164: Scan_Poly(p);
165: mat := m; rhs := v;
166: end Interpolation_System;
167:
168: procedure Interpolation_System
169: ( p : in Poly; i,b : in natural; sols : in Solution_List;
170: mat : out Standard_Complex_Matrices.Matrix;
171: rhs : out Standard_Complex_Vectors.Vector ) is
172:
173: -- DESCRIPTION :
174: -- Returns the matrix and right hand side vector of the linear system
175: -- that expresses the interpolationg conditions. Monomials with degree
176: -- one in x_i will not appear in the interpolation matrix.
177:
178: m : Standard_Complex_Matrices.Matrix(1..b,1..b);
179: v : Standard_Complex_Vectors.Vector(1..b);
180: cnt : natural := 0;
181:
182: procedure Scan_Term ( t : in Term; continue : out boolean ) is
183:
184: tmp : Solution_List := sols;
185: s : Link_to_Solution;
186:
187: begin
188: for row in m'range(1) loop -- fill in column indicated by cnt
189: s := Head_Of(tmp);
190: if cnt = 0
191: then v(row) := -t.cf*Evaluate(t,s.v);
192: else if cnt <= b and t.dg(i) /= 1
193: then m(row,cnt) := Evaluate(t,s.v);
194: else v(row) := v(row) - t.cf*Evaluate(t,s.v);
195: end if;
196: end if;
197: tmp := Tail_Of(tmp);
198: end loop;
199: if cnt = 0 or t.dg(i) /= 1
200: then cnt := cnt + 1;
201: end if;
202: continue := true;
203: end Scan_Term;
204: procedure Scan_Poly is new Visiting_Iterator(Scan_Term);
205:
206: begin
207: Scan_Poly(p);
208: mat := m; rhs := v;
209: end Interpolation_System;
210:
211: procedure Construct_Interpolant ( p : in out Poly;
212: cv : in Standard_Complex_Vectors.Vector ) is
213:
214: -- DESCRIPTION :
215: -- With the coefficients of its monomials, the interpolant will be
216: -- constructed.
217:
218: cnt : natural := 0;
219:
220: procedure Scan_Term ( t : in out Term; continue : out boolean ) is
221: begin
222: if cnt > cv'last
223: then continue := false;
224: else if cnt >= cv'first
225: then t.cf := cv(cnt);
226: end if;
227: cnt := cnt + 1;
228: continue := true;
229: end if;
230: end Scan_Term;
231: procedure Scan_Poly is new Changing_Iterator(Scan_Term);
232:
233: begin
234: Scan_Poly(p);
235: end Construct_Interpolant;
236:
237: procedure Construct_Interpolant ( p : in out Poly; i : in natural;
238: cv : in Standard_Complex_Vectors.Vector ) is
239:
240: -- DESCRIPTION :
241: -- With the coefficients of its monomials, the interpolant will be
242: -- constructed. Monomials that have degree one in x_i will be ignored.
243:
244: cnt : natural := 0;
245:
246: procedure Scan_Term ( t : in out Term; continue : out boolean ) is
247: begin
248: if cnt > cv'last
249: then continue := false;
250: else if cnt = 0 or t.dg(i) /= 1
251: then if cnt >= cv'first
252: then t.cf := cv(cnt);
253: end if;
254: cnt := cnt + 1;
255: end if;
256: continue := true;
257: end if;
258: end Scan_Term;
259: procedure Scan_Poly is new Changing_Iterator(Scan_Term);
260:
261: begin
262: Scan_Poly(p);
263: end Construct_Interpolant;
264:
265: function Interpolate ( p : Poly; b : natural; sols : Solution_List )
266: return Poly is
267:
268: -- DESCRIPTION :
269: -- Constructs the interpolating polynomial with the same monomial
270: -- structure as the given polynomial.
271:
272: res : Poly := Complex_Randomize(p);
273: mat : Standard_Complex_Matrices.Matrix(1..b,1..b);
274: rhs : Standard_Complex_Vectors.Vector(1..b);
275: ipvt : Standard_Natural_Vectors.Vector(1..b);
276: info : integer;
277:
278: begin
279: Interpolation_System(res,b,sols,mat,rhs);
280: for i in ipvt'range loop
281: ipvt(i) := i;
282: end loop;
283: lufac(mat,b,ipvt,info);
284: if info = 0
285: then lusolve(mat,b,ipvt,rhs);
286: end if;
287: Construct_Interpolant(res,rhs);
288: return res;
289: end Interpolate;
290:
291: function Interpolate ( p : Poly; i,b : natural; sols : Solution_List )
292: return Poly is
293:
294: -- DESCRIPTION :
295: -- Constructs the interpolating polynomial with the same monomial
296: -- structure as the given polynomial. The monomials that have degree
297: -- one in x_i will not appear in the interpolation matrix.
298:
299: res : Poly := Complex_Randomize(p);
300: mat : Standard_Complex_Matrices.Matrix(1..b,1..b);
301: rhs : Standard_Complex_Vectors.Vector(1..b);
302: ipvt : Standard_Natural_Vectors.Vector(1..b);
303: info : integer;
304:
305: begin
306: Interpolation_System(res,i,b,sols,mat,rhs);
307: for j in ipvt'range loop
308: ipvt(j) := j;
309: end loop;
310: lufac(mat,b,ipvt,info);
311: if info = 0
312: then lusolve(mat,b,ipvt,rhs);
313: end if;
314: Construct_Interpolant(res,i,rhs);
315: return res;
316: end Interpolate;
317:
318: -- TARGET ROUTINES :
319:
320: function Dense_Representation
321: ( p : Poly_Sys; z : partition ) return Poly_Sys is
322:
323: d : constant Standard_Integer_Matrices.Matrix := Degree_Table(p,z);
324:
325: begin
326: return Dense_Representation(p,z,d);
327: end Dense_Representation;
328:
329: function Dense_Representation
330: ( p : Poly_Sys; z : partition; d : Matrix ) return Poly_Sys is
331:
332: res : Poly_Sys(d'range(1));
333:
334: begin
335: for i in res'range loop
336: if p(i) /= Null_Poly
337: then res(i) := Dense_Representation(p(i),i,z,d);
338: end if;
339: end loop;
340: return res;
341: end Dense_Representation;
342:
343: function Independent_Representation ( p : Poly_Sys ) return Poly_Sys is
344:
345: res : Poly_Sys(p'range);
346: it : Term;
347:
348: begin
349: Copy(p,res);
350: for i in res'range loop
351: if p(i) /= Null_Poly
352: then it := Initial_Term(res(i));
353: for j in res'range loop
354: if j /= i and then (p(j) /= Null_Poly)
355: then Eliminate_Term(res(j),it.dg);
356: end if;
357: end loop;
358: end if;
359: end loop;
360: return res;
361: end Independent_Representation;
362:
363: function Independent_Roots ( p : Poly_Sys ) return natural is
364:
365: ntp : natural := 0;
366: min : natural := ntp;
367:
368: begin
369: for i in p'first..p'last loop
370: if p(i) /= Null_Poly
371: then ntp := Number_of_Terms(p(i));
372: if min = 0
373: then min := ntp;
374: elsif ntp < min
375: then min := ntp;
376: end if;
377: end if;
378: end loop;
379: if min = 0
380: then return 0;
381: else return (min-1);
382: end if;
383: end Independent_Roots;
384:
385: function Number_of_Terms ( p : Poly; i : natural ) return natural is
386:
387: -- DESCRIPTION :
388: -- Returns the number of monomials of p, without those monomials that
389: -- have degree one in x_i.
390:
391: cnt : natural := 0;
392:
393: procedure Count_Term ( t : in Term; continue : out boolean ) is
394: begin
395: if t.dg(i) /= 1
396: then cnt := cnt+1;
397: end if;
398: continue := true;
399: end Count_Term;
400: procedure Count_Terms is new Visiting_Iterator(Count_Term);
401:
402: begin
403: Count_Terms(p);
404: return cnt;
405: end Number_of_Terms;
406:
407: function Independent_Roots ( p : Poly_Sys; i : natural ) return natural is
408:
409: ntp : natural := 0;
410: min : natural := ntp;
411:
412: begin
413: for j in p'first..p'last loop
414: if p(j) /= Null_Poly
415: then ntp := Number_of_Terms(p(j),i);
416: if min = 0
417: then min := ntp;
418: elsif ntp < min
419: then min := ntp;
420: end if;
421: end if;
422: end loop;
423: if min = 0
424: then return 0;
425: else return (min-1);
426: end if;
427: end Independent_Roots;
428:
429: function Interpolate ( p : Poly_Sys; b : natural; sols : Solution_List )
430: return Poly_Sys is
431:
432: res : Poly_Sys(p'range);
433:
434: begin
435: for i in p'range loop
436: if p(i) /= Null_Poly
437: then res(i) := Interpolate(p(i),b,sols);
438: end if;
439: end loop;
440: return res;
441: end Interpolate;
442:
443: function Interpolate ( p : Poly_Sys; i,b : natural; sols : Solution_List )
444: return Poly_Sys is
445:
446: res : Poly_Sys(p'range);
447:
448: begin
449: for j in p'range loop
450: if p(j) /= Null_Poly
451: then res(j) := Interpolate(p(j),i,b,sols);
452: end if;
453: end loop;
454: return res;
455: end Interpolate;
456:
457: end Interpolating_Homotopies;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to interpolating_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 / Root_Counts / Product