![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to random_product_system.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|
Return to random_product_system.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Product |
1.1 maekawa 1: with unchecked_deallocation;
2: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
3: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
4: with Standard_Natural_Vectors;
5: with Standard_Integer_Vectors;
6: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
7: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
8: with Generic_Lists;
9:
10: package body Random_Product_System is
11:
12: -- DATA STRUCTURES :
13:
14: package List_of_Vectors is new Generic_Lists(Link_to_Vector);
15: type Equation_List is new List_of_Vectors.List;
16:
17: type Equation is record
18: first,last : Equation_List;
19: end record;
20:
21: type Equations is array(positive range <>) of Equation;
22: type Link_To_Equations is access Equations;
23: procedure free is new unchecked_deallocation (Equations,Link_To_Equations);
24:
25: -- INTERNAL DATA :
26:
27: rps : Link_To_Equations;
28:
29: Bad_Condition : constant double_float := 10.0**(-12);
30:
31: -- OPERATIONS :
32:
33: procedure Init ( n : in natural ) is
34: begin
35: rps := new Equations(1..n);
36: end Init;
37:
38: procedure Clear ( eql : in out Equation_List ) is
39:
40: temp : Equation_List := eql;
41: lv : Link_To_Vector;
42:
43: begin
44: while not Is_Null(temp) loop
45: lv := Head_Of(temp);
46: Clear(lv);
47: temp := Tail_of(temp);
48: end loop;
49: List_Of_Vectors.Clear(List_Of_Vectors.List(eql));
50: end Clear;
51:
52: procedure Clear ( eq : in out Equation ) is
53: begin
54: Clear(eq.first);
55: -- eq.last is just a pointer to the last element of eq.first;
56: -- if eq.first disappears, then also eq.last does
57: end Clear;
58:
59: procedure Clear ( eqs : in out Equations ) is
60: begin
61: for i in eqs'range loop
62: Clear(eqs(i));
63: end loop;
64: end Clear;
65:
66: procedure Clear is
67: begin
68: if rps /= null
69: then for i in rps'range loop
70: Clear(rps(i));
71: end loop;
72: free(rps);
73: end if;
74: end Clear;
75:
76: procedure Add_Hyperplane ( i : in natural; h : in Vector ) is
77:
78: eqi : Equation renames rps(i);
79: lh : Link_To_Vector := new Vector'(h);
80:
81: begin
82: if Is_Null(eqi.first)
83: then Construct(lh,eqi.first);
84: eqi.last := eqi.first;
85: else declare
86: temp : Equation_List;
87: begin
88: Construct(lh,temp);
89: Swap_Tail(eqi.last,temp);
90: eqi.last := Tail_Of(eqi.last);
91: end;
92: end if;
93: end Add_Hyperplane;
94:
95: function Dimension return natural is
96: begin
97: if rps = null
98: then return 0;
99: else return rps'last;
100: end if;
101: end Dimension;
102:
103: function Number_of_Hyperplanes ( i : natural ) return natural is
104: begin
105: if rps = null
106: then return 0;
107: else return Length_Of(rps(i).first);
108: end if;
109: end Number_Of_Hyperplanes;
110:
111: function Get_Hyperplane ( i,j : in natural ) return Link_to_Vector is
112:
113: nulvec : Link_to_Vector := null;
114:
115: begin
116: if rps = null
117: then return nulvec;
118: else declare
119: eqi : Equation_List := rps(i).first;
120: count : natural := 1;
121: begin
122: while not Is_Null(eqi) loop
123: if count = j
124: then return Head_Of(eqi);
125: else count := count + 1;
126: eqi := Tail_Of(eqi);
127: end if;
128: end loop;
129: end;
130: return nulvec;
131: end if;
132: end Get_Hyperplane;
133:
134: function Get_Hyperplane ( i,j : in natural ) return Vector is
135:
136: lres : Link_to_Vector := Get_Hyperplane(i,j);
137: nulvec : Vector(0..0) := (0..0 => Create(0.0));
138:
139: begin
140: if lres = null
141: then return nulvec;
142: else return lres.all;
143: end if;
144: end Get_Hyperplane;
145:
146: procedure Change_Hyperplane ( i,j : in natural; h : in Vector ) is
147: begin
148: if rps = null
149: then return;
150: else declare
151: eqi : Equation_List := rps(i).first;
152: lv : Link_To_Vector;
153: count : natural := 1;
154: begin
155: while not Is_Null(eqi) loop
156: if count = j
157: then lv := Head_Of(eqi);
158: for k in h'range loop
159: lv(k) := h(k);
160: end loop;
161: return;
162: else count := count + 1;
163: eqi := Tail_Of(eqi);
164: end if;
165: end loop;
166: end;
167: end if;
168: end Change_Hyperplane;
169:
170: procedure Solve ( i,n : in natural; sols,sols_last : in out Solution_List;
171: a : in out Matrix; b : in out Vector;
172: nb : in out natural ) is
173: begin
174: if i > n
175: then declare
176: aa : Matrix(a'range(1),a'range(2));
177: bb : Vector(b'range);
178: rcond : double_float;
179: ipvt : Standard_Natural_Vectors.Vector(b'range);
180: begin
181: for k in a'range(1) loop
182: for l in a'range(2) loop
183: aa(k,l) := a(k,l);
184: end loop;
185: bb(k) := b(k);
186: end loop;
187: lufco(aa,n,ipvt,rcond);
188: nb := nb + 1;
189: if abs(rcond) > Bad_Condition
190: then lusolve(aa,n,ipvt,bb);
191: declare
192: s : Solution(n);
193: begin
194: s.t := Create(0.0);
195: s.m := 1;
196: s.v := bb;
197: s.err := 0.0; s.rco := rcond; s.res := 0.0;
198: Append(sols,sols_last,s);
199: end;
200: end if;
201: exception
202: when numeric_error => return;
203: end;
204: else declare
205: eqi : Equation_List := rps(i).first;
206: h : Vector(0..n);
207: count : natural := 0;
208: begin
209: while not Is_Null(eqi) loop
210: count := count + 1;
211: h := Head_Of(eqi).all;
212: b(i) := -h(0);
213: for j in 1..n loop
214: a(i,j) := h(j);
215: end loop;
216: Solve(i+1,n,sols,sols_last,a,b,nb);
217: eqi := Tail_Of(eqi);
218: end loop;
219: end;
220: end if;
221: end Solve;
222:
223: procedure Solve ( sols : in out Solution_List; nl : out natural ) is
224:
225: n : natural := rps'last;
226: m : Matrix(1..n,1..n);
227: v : Vector(1..n);
228: num : natural := 0;
229: last : Solution_List;
230:
231: begin
232: for i in 1..n loop
233: for j in 1..n loop
234: m(i,j) := Create(0.0);
235: end loop;
236: v(i) := Create(0.0);
237: end loop;
238: Solve(1,n,sols,last,m,v,num);
239: nl := num;
240: end Solve;
241:
242: procedure Solve ( sols : in out Solution_List; nl : out natural;
243: l : in List ) is
244:
245: n : natural := rps'last;
246: m : Matrix(1..n,1..n);
247: v : Vector(1..n);
248: num : natural := 0;
249: temp : List := l;
250: pos : Standard_Integer_Vectors.Vector(1..n);
251: stop : boolean := false;
252: last : Solution_List;
253:
254: procedure PSolve ( i,n : in natural; sols,sols_last : in out Solution_List;
255: a : in out Matrix; b : in out Vector;
256: nb : in out natural ) is
257: begin
258: if i > n
259: then declare
260: aa : Matrix(a'range(1),a'range(2));
261: bb : Vector(b'range);
262: rcond : double_float;
263: ipvt : Standard_Natural_Vectors.Vector(b'range);
264: begin
265: for k in a'range(1) loop
266: for l in a'range(2) loop
267: aa(k,l) := a(k,l);
268: end loop;
269: bb(k) := b(k);
270: end loop;
271: lufco(aa,n,ipvt,rcond);
272: nb := nb + 1;
273: if abs(rcond) > Bad_Condition
274: then lusolve(aa,n,ipvt,bb);
275: declare
276: s : Solution(n);
277: begin
278: s.t := Create(0.0);
279: s.m := 1;
280: s.v := bb;
281: s.err := 0.0; s.rco := rcond; s.res := 0.0;
282: Append(sols,sols_last,s);
283: end;
284: end if;
285: if Is_Null(temp)
286: then stop := true;
287: else pos := Head_Of(temp).all;
288: temp := Tail_Of(temp);
289: end if;
290: exception
291: when numeric_error => return;
292: end;
293: else declare
294: eqi : Equation_List := rps(i).first;
295: h : Vector(0..n);
296: count : natural := 0;
297: begin
298: while not Is_Null(eqi) loop
299: count := count + 1;
300: if count = pos(i)
301: then h := Head_Of(eqi).all;
302: b(i) := -h(0);
303: for j in 1..n loop
304: a(i,j) := h(j);
305: end loop;
306: --put("eq"); put(i,1); put(count,1); put(" ");
307: PSolve(i+1,n,sols,sols_last,a,b,nb);
308: end if;
309: exit when stop;
310: eqi := Tail_Of(eqi);
311: end loop;
312: end;
313: end if;
314: end PSolve;
315:
316: begin
317: if not Is_Null(temp)
318: then pos := Head_Of(temp).all;
319: temp := Tail_Of(temp);
320: for i in 1..n loop
321: for j in 1..n loop
322: m(i,j) := Create(0.0);
323: end loop;
324: v(i) := Create(0.0);
325: end loop;
326: PSolve(1,n,sols,last,m,v,num);
327: nl := num;
328: end if;
329: end Solve;
330:
331: function Polynomial ( h : Vector ) return Poly is
332:
333: res : Poly;
334: t : Term;
335: n : natural := h'last;
336:
337: begin
338: for j in 0..n loop
339: if h(j) /= Create(0.0)
340: then t.dg := new Standard_Natural_Vectors.Vector'(1..n => 0);
341: t.cf := h(j);
342: if j /= 0
343: then t.dg(j) := 1;
344: end if;
345: Add(res,t);
346: Clear(t.dg);
347: end if;
348: end loop;
349: return res;
350: end Polynomial;
351:
352: function Create ( i : in natural ) return Poly is
353:
354: eql : Equation_List := rps(i).first;
355: hyp,res : Poly := Null_Poly;
356:
357: begin
358: while not Is_Null(eql) loop
359: hyp := Polynomial(Head_Of(eql).all);
360: if res = Null_Poly
361: then Copy(hyp,res);
362: else Mul(res,hyp);
363: end if;
364: Clear(hyp);
365: eql := Tail_Of(eql);
366: end loop;
367: return res;
368: end Create;
369:
370: function Polynomial_System return Poly_Sys is
371:
372: res : Poly_Sys(rps'range);
373:
374: begin
375: for i in rps'range loop
376: res(i) := Create(i);
377: end loop;
378: return res;
379: end Polynomial_System;
380:
381: end Random_Product_System;