![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to fewnomial_system_solvers.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 / Implift|
Return to fewnomial_system_solvers.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift |
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_Integer_Vectors;
5: with Standard_Integer_VecVecs; use Standard_Integer_VecVecs;
6: with Standard_Complex_Vectors; use Standard_Complex_Vectors;
7: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
8: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
9: with Standard_Complex_Laur_Polys; use Standard_Complex_Laur_Polys;
10: with Binomial_System_Solvers; use Binomial_System_Solvers;
11:
12: package body Fewnomial_System_Solvers is
13:
14: function Equal ( v : Standard_Integer_Vectors.Link_to_Vector; d : Degrees )
15: return boolean is
16: begin
17: return Standard_Integer_Vectors.Equal
18: (v,Standard_Integer_Vectors.Link_to_Vector(d));
19: end Equal;
20:
21: function Is_Fewnomial_System ( p : Laur_Sys ) return boolean is
22:
23: da : VecVec(p'range);
24: nb : natural;
25: n : natural := p'last - p'first + 1;
26: cnst : Standard_Integer_Vectors.Link_to_Vector;
27: res : boolean;
28:
29: procedure Scan_Term ( t : in Term; cont : out boolean ) is
30:
31: found : boolean;
32:
33: begin
34: found := false;
35: for i in da'range loop
36: exit when i > nb;
37: if Equal(da(i),t.dg)
38: then found := true;
39: exit;
40: end if;
41: end loop;
42: if not found
43: then if nb < n
44: then nb := nb + 1;
45: found := true;
46: da(nb) := new Standard_Integer_Vectors.Vector(t.dg'range);
47: for j in t.dg'range loop
48: da(nb)(j) := t.dg(j);
49: end loop;
50: elsif nb < n + 1
51: then nb := nb + 1;
52: cnst := new Standard_Integer_Vectors.Vector(t.dg'range);
53: found := true;
54: for j in t.dg'range loop
55: cnst(j) := t.dg(j);
56: end loop;
57: elsif Equal(cnst,t.dg)
58: then found := true;
59: else res := false;
60: end if;
61: end if;
62: cont := found;
63: end Scan_Term;
64: procedure Scan_Terms is new Visiting_Iterator (Scan_Term);
65:
66: begin
67: res := true;
68: nb := 0;
69: for i in p'range loop
70: Scan_Terms(p(i));
71: exit when not res;
72: end loop;
73: Clear(da);
74: Standard_Integer_Vectors.Clear(cnst);
75: return res;
76: end Is_Fewnomial_System;
77:
78: procedure Scale ( cnst : in Standard_Integer_Vectors.Link_to_Vector;
79: da : in out VecVec ) is
80:
81: -- DESCRIPTION :
82: -- If cnst is not equal to the zero degrees,
83: -- then all entries in da will be shifted along cnst.
84:
85: zero : boolean;
86:
87: begin
88: zero := true;
89: for i in cnst'range loop
90: if cnst(i) /= 0
91: then zero := false; exit;
92: end if;
93: end loop;
94: if not zero
95: then for i in da'range loop
96: Standard_Integer_Vectors.Sub(da(i),cnst);
97: end loop;
98: end if;
99: end Scale;
100:
101: procedure Initialize ( p : in Laur_Sys; n : in natural;
102: a : in out matrix; b : in out vector;
103: da : in out VecVec; fail : out boolean ) is
104:
105: -- DESCRIPTION :
106: -- initializes the data needed for the solver;
107: -- fail becomes true if the system p has too many different terms.
108:
109: cnst : Standard_Integer_Vectors.Link_to_Vector;
110: da_numb,cnt : natural;
111: fl : boolean;
112:
113: use Standard_Integer_Vectors;
114:
115: procedure Search_Term ( t : in Term; cont : out boolean ) is
116:
117: found : boolean;
118:
119: begin
120: found := false;
121: for i in da'range loop
122: exit when i > da_numb;
123: if Equal(da(i),t.dg)
124: then found := true;
125: a(cnt,i) := t.cf;
126: exit;
127: end if;
128: end loop;
129: if not found
130: then
131: if da_numb < n
132: then da_numb := da_numb + 1;
133: da(da_numb) := new Standard_Integer_Vectors.Vector(t.dg'range);
134: for k in da(da_numb)'range loop
135: da(da_numb)(k) := t.dg(k);
136: end loop;
137: found := true;
138: a(cnt,da_numb) := t.cf;
139: elsif da_numb < n+1
140: then da_numb := da_numb + 1;
141: cnst := new Standard_Integer_Vectors.Vector(t.dg'range);
142: for k in cnst'range loop
143: cnst(k) := t.dg(k);
144: end loop;
145: found := true;
146: b(cnt) := -t.cf;
147: elsif Equal(cnst,t.dg)
148: then found := true;
149: b(cnt) := -t.cf;
150: end if;
151: end if;
152: if not found
153: then cont := false;
154: fl := true;
155: else cont := true;
156: fl := false;
157: end if;
158: end Search_Term;
159: procedure Search is new Visiting_Iterator(Search_Term);
160:
161: begin
162: da_numb := 0;
163: for i in p'range loop
164: for j in a'range(2) loop
165: a(i,j) := Create(0.0);
166: end loop;
167: b(i) := Create(0.0);
168: cnt := i;
169: Search(p(i));
170: exit when fl;
171: end loop;
172: if not fl and then (cnst /= null)
173: then Scale(cnst,da);
174: Standard_Integer_Vectors.Clear(cnst);
175: fail := false;
176: else fail := true;
177: end if;
178: end Initialize;
179:
180: procedure Solve ( p : in Laur_Sys; sols : in out Solution_List;
181: fail : out boolean ) is
182:
183: n : natural := p'length;
184: da : VecVec(p'range);
185: a : Matrix(1..n,1..n);
186: b : Vector(1..n);
187: fl : boolean;
188:
189: begin
190: Initialize(p,n,a,b,da,fl);
191: if fl
192: then fail := true;
193: else fail := false;
194: declare
195: piv : Standard_Natural_Vectors.Vector(1..n);
196: info : integer;
197: begin
198: lufac(a,n,piv,info);
199: if info = 0
200: then lusolve(a,n,piv,b);
201: fl := false;
202: for i in b'range loop
203: if AbsVal(b(i)) + 1.0 = 1.0
204: then fl := true;
205: exit;
206: end if;
207: end loop;
208: if not fl
209: then Solve(da,b,n,sols);
210: end if;
211: end if;
212: end;
213: end if;
214: Clear(da);
215: end Solve;
216:
217: end Fewnomial_System_Solvers;