Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Polynomials/standard_complex_substitutors.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
2: with Standard_Natural_Vectors;
3:
4: package body Standard_Complex_Substitutors is
5:
6: function Substitute ( k : integer; c : Complex_Number; t : Term )
7: return Term is
8:
9: res : Term;
10:
11: begin
12: res.cf := t.cf;
13: for l in 1..t.dg(k) loop
14: res.cf := res.cf*c;
15: end loop;
16: res.dg := new Standard_Natural_Vectors.Vector(t.dg'first..t.dg'last-1);
17: for l in t.dg'range loop
18: if l < k
19: then res.dg(l) := t.dg(l);
20: elsif l > k
21: then res.dg(l-1) := t.dg(l);
22: end if;
23: end loop;
24: return res;
25: end Substitute;
26:
27: function Substitute ( k : integer; c : Complex_Number; p : Poly )
28: return Poly is
29: res : Poly;
30:
31: procedure Substitute_Term ( t : in Term; cont : out boolean ) is
32:
33: st : Term := Substitute(k,c,t);
34:
35: begin
36: Add(res,st);
37: Clear(st);
38: cont := true;
39: end Substitute_Term;
40: procedure Substitute_Terms is new Visiting_Iterator (Substitute_Term);
41:
42: begin
43: Substitute_Terms(p);
44: return res;
45: end Substitute;
46:
47: function Substitute ( k : integer; c : Complex_Number; p : Poly_Sys )
48: return Poly_Sys is
49:
50: res : Poly_Sys(p'range);
51:
52: begin
53: for l in res'range loop
54: res(l) := Substitute(k,c,p(l));
55: end loop;
56: return res;
57: end Substitute;
58:
59: procedure Substitute ( k : in integer; c : in Complex_Number;
60: t : in out Term ) is
61:
62: tmp : Standard_Natural_Vectors.Vector(t.dg'first..t.dg'last-1);
63:
64: begin
65: for l in 1..t.dg(k) loop
66: t.cf := t.cf*c;
67: end loop;
68: for l in t.dg'range loop
69: if l < k
70: then tmp(l) := t.dg(l);
71: elsif l > k
72: then tmp(l-1) := t.dg(l);
73: end if;
74: end loop;
75: Clear(t);
76: t.dg := new Standard_Natural_Vectors.Vector'(tmp);
77: end Substitute;
78:
79: procedure Substitute ( k : in integer; c : in Complex_Number;
80: p : in out Poly ) is
81:
82: -- NOTE :
83: -- An obvious thing to do would be to visit and change all terms,
84: -- and leaving the term order unchanged.
85:
86: procedure Substitute_Term ( t : in out Term; cont : out boolean ) is
87: begin
88: Substitute(k,c,t);
89: cont := true;
90: end Substitute_Term;
91: procedure Substitute_Terms is new Changing_Iterator ( Substitute_Term );
92:
93: begin
94: Substitute_Terms(p);
95: end Substitute;
96:
97: procedure Substitute ( k : in integer; c : in Complex_Number;
98: p : in out Poly_Sys ) is
99: begin
100: for l in p'range loop
101: Substitute(k,c,p(l));
102: end loop;
103: end Substitute;
104:
105: procedure Purify ( p : in out Poly; tol : in double_float ) is
106:
107: -- DESCRIPTION :
108: -- All terms of which the coefficient are in AbsVal smaller
109: -- than tol are deleted.
110:
111: procedure Purify_Term (t : in out Term; continue : out boolean) is
112: begin
113: if AbsVal(t.cf) < tol
114: then t.cf := Create(0.0);
115: end if;
116: continue := true;
117: end Purify_Term;
118: procedure Purify_Terms is new Changing_Iterator(Purify_Term);
119:
120: begin
121: Purify_Terms(p);
122: if Number_Of_Unknowns(p) = 0
123: then Clear(p);
124: p := Null_Poly;
125: end if;
126: end Purify;
127:
128: function Substitute_Factor ( k : integer; h : Vector ) return Poly is
129:
130: -- DESCRIPTION :
131: -- returns the polynomial which will replace the kth unknown.
132:
133: res : Poly;
134: rt : Term;
135:
136: begin
137: rt.dg := new Standard_Natural_Vectors.Vector'((h'first+1)..(h'last-1) => 0);
138: rt.cf := -h(0)/h(k);
139: res := Create(rt);
140: for i in rt.dg'range loop
141: rt.dg(i) := 1;
142: if i < k
143: then rt.cf := -h(i)/h(k);
144: else rt.cf := -h(i+1)/h(k);
145: end if;
146: if AbsVal(rt.cf) > 10.0**(-10)
147: then Add(res,rt);
148: end if;
149: rt.dg(i) := 0;
150: end loop;
151: Clear(rt);
152: return res;
153: end Substitute_Factor;
154:
155: function Substitute ( k : integer; h : Vector; p : Poly ) return Poly is
156:
157: res,sub : Poly;
158:
159: begin
160: sub := Substitute_Factor(k,h);
161: res := Substitute(k,sub,p);
162: Clear(sub);
163: return res;
164: end Substitute;
165:
166: function Substitute ( k : integer; s,p : Poly ) return Poly is
167:
168: res : Poly := Null_Poly;
169:
170: procedure Substitute_Term ( t : in Term; continue : out boolean ) is
171:
172: rt : Term;
173: fac : Poly;
174:
175: begin
176: rt.cf := t.cf;
177: rt.dg := new Standard_Natural_Vectors.Vector(t.dg'first..(t.dg'last-1));
178: for i in rt.dg'range loop
179: if i < k
180: then rt.dg(i) := t.dg(i);
181: else rt.dg(i) := t.dg(i+1);
182: end if;
183: end loop;
184: if t.dg(k) = 0
185: then Add(res,rt);
186: else fac := Create(rt);
187: for i in 1..t.dg(k) loop
188: Mul(fac,s);
189: end loop;
190: Purify(fac,10.0**(-10));
191: Add(res,fac);
192: Clear(fac);
193: end if;
194: Clear(rt);
195: continue := true;
196: end Substitute_Term;
197: procedure Substitute_Terms is new Visiting_Iterator (Substitute_Term);
198:
199: begin
200: Substitute_Terms(p);
201: return res;
202: end Substitute;
203:
204: procedure Substitute ( k : in integer; h : in Vector; p : in out Poly ) is
205:
206: res : Poly;
207:
208: begin
209: res := Substitute(k,h,p);
210: Clear(p); Copy(res,p); Clear(res);
211: end Substitute;
212:
213: procedure Substitute ( k : in integer; s : in Poly; p : in out Poly ) is
214:
215: res : Poly;
216:
217: begin
218: res := Substitute(k,s,p);
219: Clear(p); Copy(res,p); Clear(res);
220: end Substitute;
221:
222: function Substitute ( k : integer; h : Vector; p : Poly_Sys )
223: return Poly_Sys is
224:
225: res : Poly_Sys(p'range);
226: s : Poly := Substitute_Factor(k,h);
227:
228: begin
229: res := Substitute(k,s,p);
230: Clear(s);
231: return res;
232: end Substitute;
233:
234: procedure Substitute ( k : in integer; h : in Vector;
235: p : in out Poly_Sys ) is
236:
237: s : Poly := Substitute_Factor(k,h);
238:
239: begin
240: Substitute(k,s,p);
241: Clear(s);
242: end Substitute;
243:
244: function Substitute ( k : integer; s : Poly; p : Poly_Sys ) return Poly_Sys is
245:
246: res : Poly_Sys(p'range);
247:
248: begin
249: for i in p'range loop
250: res(i) := Substitute(k,s,p(i));
251: end loop;
252: return res;
253: end Substitute;
254:
255: procedure Substitute ( k : in integer; s : in Poly; p : in out Poly_Sys ) is
256: begin
257: for i in p'range loop
258: Substitute(k,s,p(i));
259: end loop;
260: end Substitute;
261:
262: end Standard_Complex_Substitutors;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to standard_complex_substitutors.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Polynomials