![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to scaling.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Homotopy|
Return to scaling.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Homotopy |
1.1 maekawa 1: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
2: with Standard_Mathematical_Functions; use Standard_Mathematical_Functions;
3: with Standard_Natural_Vectors;
4: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
5: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
6:
7: package body Scaling is
8:
9: function log ( b : in natural; x : in double_float ) return double_float is
10: begin
11: return ( LOG10(x)/LOG10(double_float(b)) );
12: end log;
13:
14: procedure Scale ( p : in out Poly ) is
15:
16: sum : double_float := 0.0;
17: number : natural := 0;
18: factor : Complex_Number;
19:
20: procedure Add_To_Sum ( t : in Term; continue : out boolean ) is
21: begin
22: sum := sum + AbsVal(t.cf);
23: number := number + 1;
24: continue := true;
25: end Add_To_Sum;
26: procedure Compute_Sum is new Visiting_Iterator(Add_To_Sum);
27:
28: begin
29: Compute_Sum(p);
30: factor := Create(double_float(number)/sum);
31: Mul(p,factor);
32: end Scale;
33:
34: procedure Scale ( s : in out Poly_Sys ) is
35: begin
36: for i in s'range loop
37: scale(s(i));
38: end loop;
39: end Scale;
40:
41: procedure Scale ( s : in out Poly_Sys; bas : in natural := 2;
42: diff : in boolean; cond : out double_float;
43: sccff : out vector ) is
44:
45: r1,r2,target : Poly;
46: n : natural := s'last - s'first + 1;
47: nm : natural := 2*n;
48: mat : Matrix(1..nm,1..nm);
49: right : Vector(1..nm);
50:
51: function Center_Coefficients ( s : in Poly_Sys; bas : in natural )
52: return Poly is
53: r1,r1i : Poly;
54: n : natural := s'last - s'first + 1;
55:
56: procedure Scan ( p : in Poly; i : in natural ) is
57:
58: init : Poly;
59: t_init : Term;
60:
61: procedure Scan_Term ( t : in Term; continue : out boolean ) is
62:
63: tt : Term;
64: temp : Poly;
65:
66: begin
67: Copy(init,temp);
68: continue := true;
69: for k in t.dg'range loop
70: if t.dg(k) /= 0
71: then tt.cf := Create(double_float(t.dg(k)));
72: tt.dg := new Standard_Natural_Vectors.Vector'(1..2*n => 0);
73: tt.dg(k) := 1;
74: Add(temp,tt);
75: Clear(tt);
76: end if;
77: end loop;
78: tt.cf := Create(log(bas,AbsVal(t.cf)));
79: tt.dg := new Standard_Natural_Vectors.Vector'(1..2*n => 0);
80: Add(temp,tt);
81: Clear(tt);
82: Mul(temp,temp);
83: Add(r1i,temp);
84: Clear(temp);
85: end Scan_Term;
86: procedure Scan_Terms is new Visiting_Iterator(Scan_Term);
87:
88: begin
89: t_init.cf := Create(1.0);
90: t_init.dg := new Standard_Natural_Vectors.Vector'(1..2*n => 0);
91: t_init.dg(n+i) := 1;
92: init := Create(t_init);
93: Clear(t_init);
94: Scan_Terms(p);
95: end Scan;
96:
97: begin
98: r1 := Null_Poly;
99: for i in s'range loop
100: Scan(s(i),i);
101: Add(r1,r1i);
102: Clear(r1i);
103: end loop;
104: return r1;
105: end Center_Coefficients;
106:
107: function Reduce_Diff ( s : in Poly_Sys; bas : in natural ) return Poly is
108:
109: r2,r2i : Poly;
110:
111: procedure Scan2 (p : in Poly; t : in Term; nr : in natural) is
112:
113: count : natural := 0;
114:
115: procedure Scan2_Term (t2 : in Term; continue : out boolean) is
116:
117: tt : Term;
118: temp : Poly := Null_Poly;
119:
120: begin
121: continue := true;
122: count := count + 1;
123: if count > nr
124: then
125: for i in t2.dg'range loop
126: if t.dg(i) /= t2.dg(i)
127: then tt.dg := new Standard_Natural_Vectors.Vector'(1..2*n => 0);
128: tt.dg(i) := 1;
129: tt.cf := Create(double_float(t.dg(i)-t2.dg(i)));
130: Add(temp,tt);
131: Clear(tt);
132: end if;
133: end loop;
134: end if;
135: tt.dg := new Standard_Natural_Vectors.Vector'(1..2*n => 0);
136: tt.cf := Create(log(bas,(AbsVal(t.cf)/AbsVal(t2.cf))));
137: Add(temp,tt);
138: Clear(tt);
139: Mul(temp,temp);
140: Add(r2i,temp);
141: Clear(temp);
142: end Scan2_Term;
143: procedure Scan2_Terms is new Visiting_Iterator(Scan2_Term);
144:
145: begin
146: Scan2_Terms(p);
147: end Scan2;
148:
149: procedure Scan ( p : in Poly ) is
150:
151: nr : natural := 0;
152:
153: procedure Scan_Term ( t : in Term; continue : out boolean ) is
154: begin
155: nr := nr + 1;
156: continue := true;
157: Scan2(p,t,nr);
158: end Scan_Term;
159: procedure Scan_Terms is new Visiting_Iterator(Scan_Term);
160:
161: begin
162: Scan_Terms(p);
163: end Scan;
164:
165: begin
166: r2 := Null_Poly;
167: for i in s'range loop
168: Scan(s(i));
169: Add(r1,r2i);
170: Clear(r2i);
171: end loop;
172: return r2;
173: end Reduce_Diff;
174:
175: procedure Make_Linear_System ( r : in Poly; mat : out Matrix;
176: right : out Vector ) is
177: drj : Poly;
178:
179: procedure Init_Linear_System ( m : out Matrix; r : out Vector ) is
180: begin
181: for i in m'range(1) loop
182: r(i) := Create(0.0);
183: for j in m'range(2) loop
184: m(i,j) := Create(0.0);
185: end loop;
186: end loop;
187: end Init_Linear_System;
188:
189: procedure Scan ( p : in Poly; j : in natural ) is
190:
191: procedure Scan_Term ( t : in Term; continue : out boolean ) is
192: begin
193: continue := true;
194: for i in t.dg'range loop
195: if t.dg(i) = 1
196: then mat(j,i) := t.cf;
197: return;
198: end if;
199: end loop;
200: right(j) := -t.cf;
201: end Scan_Term;
202: procedure Scan_Terms is new Visiting_Iterator (Scan_Term);
203:
204: begin
205: Scan_Terms(p);
206: end Scan;
207:
208: begin
209: Init_Linear_System(mat,right);
210: for j in 1..Number_Of_Unknowns(r) loop
211: drj := Standard_Complex_Polynomials.Diff(r,j);
212: Scan(drj,j);
213: end loop;
214: end Make_Linear_System;
215:
216: procedure Scale ( s : in out Poly_Sys; bas : in natural;
217: mat : in out Matrix; right : in out Vector;
218: cond : out double_float ) is
219:
220: n : natural := s'last - s'first + 1;
221: ipvt : Standard_Natural_Vectors.Vector(1..2*n);
222:
223: procedure Scale ( p : in out Poly; ip : in natural;
224: scalingcoeff : in Vector; bas : in natural ) is
225:
226: procedure Scale_Term ( t : in out Term; continue : out boolean ) is
227:
228: exp : double_float := 0.0;
229:
230: begin
231: exp := REAL_PART(scalingcoeff(n+ip));
232: for k in t.dg'range loop
233: exp := exp + double_float(t.dg(k))*REAL_PART(scalingcoeff(k));
234: end loop;
235: t.cf := t.cf * Create(double_float(bas) ** exp);
236: continue := true;
237: end Scale_Term;
238: procedure Scale_Terms is new Changing_Iterator (Scale_Term);
239:
240: begin
241: Scale_Terms(p);
242: end Scale;
243:
244: begin
245: lufco(mat,2*n,ipvt,cond);
246: lusolve(mat,2*n,ipvt,right);
247: for i in s'range loop
248: Scale(s(i),i,right,bas);
249: end loop;
250: end Scale;
251:
252: begin
253: r1 := center_coefficients(s,bas);
254: if diff
255: then r2 := reduce_diff(s,bas);
256: target := r1 + r2;
257: clear(r1); clear(r2);
258: else copy(r1,target); clear(r1);
259: end if;
260: Make_Linear_System(target,mat,right);
261: clear(target);
262: scale(s,bas,mat,right,cond);
263: sccff := right;
264: end Scale;
265:
266: procedure Scale ( basis : in natural; sccff : in Vector;
267: s : in out Solution ) is
268: begin
269: for i in s.v'range loop
270: s.v(i) := Create(double_float(basis)**REAL_PART(sccff(i))) * s.v(i);
271: end loop;
272: end Scale;
273:
274: procedure Scale ( basis : in natural; sccff : in Vector;
275: sols : in out Solution_List ) is
276: begin
277: if Is_Null(sols)
278: then null;
279: else declare
280: temp : Solution_List := sols;
281: n : natural := Head_Of(sols).n;
282: s : Solution(n);
283: l : Link_To_Solution;
284: begin
285: while not Is_Null(temp) loop
286: l := Head_Of(temp);
287: s := l.all;
288: Scale(basis,sccff,s);
289: Clear(l);
290: l := new Solution'(s);
291: Set_Head(temp,l);
292: temp := Tail_Of(temp);
293: end loop;
294: end;
295: end if;
296: end Scale;
297:
298: end Scaling;