![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to orbits_of_solutions.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 / Symmetry|
Return to orbits_of_solutions.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Symmetry |
1.1 maekawa 1: with unchecked_deallocation;
2: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
3: with Standard_Complex_Vectors;
4: with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
5: with Permute_Operations; use Permute_Operations;
6:
7: package body Orbits_of_Solutions is
8:
9: -- AUXILIAIRIES :
10:
11: function Equal_upon_Sign
12: ( x,y : Standard_Complex_Vectors.Vector; tol : double_float )
13: return boolean is
14:
15: -- DESCRIPTION :
16: -- Returns true if for all i : | x(i) +- y(i) | <= tol.
17:
18: begin
19: for i in x'range loop
20: if (AbsVal(x(i)-y(i)) > tol) and then (AbsVal(x(i)+y(i)) > tol)
21: then return false;
22: end if;
23: end loop;
24: return true;
25: end Equal_upon_Sign;
26:
27: procedure Is_In ( sols : in out Solution_List; s : in Solution;
28: sign : in boolean; tol : in double_float;
29: isin : out boolean ) is
30:
31: -- DESCRIPTION :
32: -- The output variable isin isin becomes true if one of the solutions
33: -- in the list sols equals s. When this is the case, the multiplicity
34: -- of the corresponding solution is augmented.
35: -- Otherwise, isin equals false on return.
36:
37: s1 : Solution(s.n);
38: temp : Solution_List;
39:
40: begin
41: temp := sols;
42: while not Is_Null(temp) loop
43: s1 := Head_Of(temp).all;
44: if (sign and then Equal_upon_Sign(s.v,s1.v,tol))
45: or else (not sign and then Equal(s.v,s1.v,tol))
46: then Head_Of(temp).m := Head_Of(temp).m + 1;
47: isin := true;
48: return;
49: end if;
50: temp := Tail_Of(temp);
51: end loop;
52: isin := false;
53: end Is_In;
54:
55: -- CONSTRUCTORS :
56:
57: procedure Analyze ( l : in List_of_Permutations; sign : in boolean;
58: tol : in double_float; sols : in out Solution_List ) is
59: begin
60: if not Is_Null(sols)
61: then declare
62: res,res_last,tmpsols : Solution_List;
63: n : natural := Head_Of(sols).n;
64: s1,s2 : Solution(n);
65: lp : Link_to_Permutation;
66: tmpl : List_of_Permutations;
67: found : boolean;
68: begin
69: tmpsols := sols;
70: while not Is_Null(tmpsols) loop
71: s1 := Head_Of(tmpsols).all;
72: tmpl := l;
73: while not Is_Null(tmpl) loop
74: lp := Head_Of(tmpl);
75: s2.t := s1.t;
76: s2.m := 1;
77: s2.v := Permutation(lp.all)*s1.v;
78: Is_In(res,s2,sign,tol,found);
79: exit when found;
80: tmpl := Tail_Of(tmpl);
81: end loop;
82: if not found
83: then Append(res,res_last,s1);
84: end if;
85: tmpsols := Tail_Of(tmpsols);
86: end loop;
87: Clear(sols); sols := res;
88: end;
89: end if;
90: end Analyze;
91:
92: function Permutable ( s : Solution; sign : boolean; tol : double_float;
93: sols : Solution_List ) return boolean is
94:
95: -- DESCRIPTION :
96: -- Returns true if the solutions s can be permuted into one of the
97: -- solutions of the given list.
98:
99: tmp : Solution_List := sols;
100:
101: begin
102: while not Is_Null(tmp) loop
103: if sign and then Sign_Permutable(s.v,Head_Of(tmp).v,tol)
104: then return true;
105: elsif Permutable(s.v,Head_Of(tmp).v,tol)
106: then return true;
107: else tmp := Tail_Of(tmp);
108: end if;
109: end loop;
110: return false;
111: end Permutable;
112:
113: procedure Permutable
114: ( s : in Solution; sign : in boolean; tol : in double_float;
115: sols : in out Solution_List; isin : out boolean ) is
116:
117: -- DESCRIPTION :
118: -- The output variable isin becomes true when the given solution can be
119: -- permuted into one of the solutions of the given list.
120: -- In this case, also the multiplicity of the corresponding solutions is
121: -- augmented. Otherwise, isin is false on return.
122:
123: tmp : Solution_List := sols;
124: ls1 : Link_to_Solution;
125:
126: begin
127: while not Is_Null(tmp) loop
128: ls1 := Head_Of(tmp);
129: if (sign and then Sign_Permutable(s.v,ls1.v,tol))
130: or else (not sign and then Permutable(s.v,ls1.v,tol))
131: then isin := true;
132: ls1.m := ls1.m + 1; Set_Head(tmp,ls1);
133: return;
134: else tmp := Tail_Of(tmp);
135: end if;
136: end loop;
137: isin := false;
138: end Permutable;
139:
140: function Generating ( sols : Solution_List; sign : boolean;
141: tol : double_float ) return Solution_List is
142:
143: tmp,gensols,gensols_last : Solution_List;
144:
145: begin
146: tmp := sols;
147: while not Is_Null(tmp) loop
148: declare
149: ls : Link_to_Solution := Head_Of(tmp);
150: isin : boolean;
151: begin
152: Permutable(ls.all,sign,tol,gensols,isin);
153: if not isin
154: then Append(gensols,gensols_last,ls.all);
155: end if;
156: end;
157: tmp := Tail_Of(tmp);
158: end loop;
159: return gensols;
160: end Generating;
161:
162: procedure Orbit_Structure ( s : in Solution; tol : in double_float;
163: orbit : in out Permutation;
164: nbdiff : out natural ) is
165:
166: nb : natural := 1;
167: found : boolean;
168:
169: begin
170: orbit(1) := nb;
171: for i in 2..s.n loop
172: found := false;
173: for j in 1..(i-1) loop
174: if AbsVal(s.v(i) - s.v(j)) < tol
175: then orbit(i) := orbit(j);
176: found := true;
177: exit;
178: end if;
179: end loop;
180: if not found
181: then nb := nb + 1;
182: orbit(i) := nb;
183: end if;
184: end loop;
185: nbdiff := nb;
186: end Orbit_Structure;
187:
188: function Orbits ( sols : Solution_List; tol : double_float )
189: return permutation is
190:
191: ind : natural;
192: tmp : Solution_List := sols;
193: n : natural := Head_Of(sols).n;
194: o,orb : Permutation(1..n);
195: sol : Solution(n);
196:
197: begin
198: orb := (1..n => 0);
199: o := orb;
200: while not Is_Null(tmp) loop
201: sol := Head_Of(tmp).all;
202: Orbit_Structure(sol,tol,o,ind);
203: orb(ind) := orb(ind) + sol.m;
204: tmp := Tail_Of(tmp);
205: end loop;
206: return orb;
207: end Orbits;
208:
209: procedure Orbits ( grp : in List_of_Permutations; tol : in double_float;
210: sols : in out Solution_List;
211: lorb : in out List_of_Orbits ) is
212:
213: tmp : Solution_List;
214: n : natural := Head_Of(sols).n;
215:
216: begin
217: Analyze(grp,false,tol,sols);
218: tmp := sols;
219: while not Is_Null(tmp) loop
220: declare
221: sol : Solution(n) := Head_Of(tmp).all;
222: orbi : Orbit(n);
223: begin
224: Orbit_Structure(sol,tol,orbi.orb,orbi.nbdiff);
225: orbi.nbgen := 1;
226: orbi.nbsols := sol.m;
227: declare
228: tmp2 : List_of_Orbits := lorb;
229: found : boolean := false;
230: begin
231: while not Is_Null(tmp2) loop
232: declare
233: orbi2 : Orbit(n) := Head_Of(tmp2).all;
234: lorbi2 : Link_to_Orbit := Head_Of(tmp2);
235: begin
236: if orbi2.nbdiff = orbi.nbdiff
237: and then Same_Orbit(orbi.orb,orbi2.orb)
238: then lorbi2.nbgen := orbi2.nbgen + 1;
239: lorbi2.nbsols := orbi2.nbsols + sol.m;
240: found := true;
241: exit;
242: else tmp2 := Tail_Of(tmp2);
243: end if;
244: end;
245: end loop;
246: if not found
247: then declare
248: liorbi : Link_to_Orbit := new Orbit'(orbi);
249: begin
250: Construct(liorbi,lorb);
251: end;
252: end if;
253: end;
254: end;
255: tmp := Tail_Of(tmp);
256: end loop;
257: end Orbits;
258:
259: -- SELECTOR :
260:
261: function Same_Orbit ( orb1,orb2 : Permutation ) return boolean is
262:
263: f1,f2 : Permutation(1..orb1'last);
264: found : boolean;
265:
266: begin
267: -- construct frequencies :
268: f1 := (f1'range => 0);
269: for i in orb1'range loop
270: f1(orb1(i)) := f1(orb1(i)) + 1;
271: end loop;
272: f2 := (f2'range => 0);
273: for i in orb2'range loop
274: f2(orb2(i)) := f2(orb2(i)) + 1;
275: end loop;
276: -- compare the frequencies :
277: for i in f1'range loop
278: found := false;
279: for j in f2'range loop
280: if f1(i) = f2(j)
281: then found := true;
282: exit;
283: end if;
284: end loop;
285: if not found
286: then return false;
287: end if;
288: end loop;
289: for i in f2'range loop
290: found := false;
291: for j in f1'range loop
292: if f2(i) = f1(j)
293: then found := true;
294: exit;
295: end if;
296: end loop;
297: if not found
298: then return false;
299: end if;
300: end loop;
301: return true;
302: end Same_Orbit;
303:
304: -- DESTRUCTOR :
305:
306: procedure Clear ( lorb : in out List_of_Orbits ) is
307:
308: procedure free is new unchecked_deallocation(Orbit,Link_to_Orbit);
309:
310: tmp : List_of_Orbits := lorb;
311: begin
312: while not Is_Null(tmp) loop
313: declare
314: lior : Link_to_Orbit := Head_Of(tmp);
315: begin
316: free(lior);
317: end;
318: tmp := Tail_Of(tmp);
319: end loop;
320: Lists_of_Orbits.Clear(Lists_of_Orbits.List(lorb));
321: end Clear;
322:
323: end Orbits_of_Solutions;