Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Symmetry/templates.adb, Revision 1.1.1.1
1.1 maekawa 1: with unchecked_deallocation;
2: with Generic_Lists;
3: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
4: with Standard_Random_Numbers; use Standard_Random_Numbers;
5: with Standard_Complex_Vectors;
6: with Random_Product_System;
7:
8: package body Templates is
9:
10: -- DATA STRUCTURES :
11:
12: package List_of_Vectors is new Generic_Lists(Link_to_Vector);
13: type Equation_List is new List_of_Vectors.List;
14:
15: type Equation is record
16: first,last : Equation_List;
17: end record;
18:
19: type Equations is array(positive range <>) of Equation;
20: type Link_To_Equations is access Equations;
21: procedure free is new unchecked_deallocation (Equations,Link_To_Equations);
22:
23: -- INTERNAL DATA :
24:
25: rps : Link_To_Equations;
26:
27: --------------------
28: -- CONSTRUCTORS --
29: --------------------
30:
31: procedure Create ( n : in natural ) is
32: begin
33: rps := new Equations(1..n);
34: end Create;
35:
36: procedure Add_Hyperplane ( i : in natural; h : in Vector ) is
37:
38: eqi : Equation renames rps(i);
39: lh : Link_To_Vector := new Vector'(h);
40:
41: begin
42: if Is_Null(eqi.first)
43: then Construct(lh,eqi.first);
44: eqi.last := eqi.first;
45: else declare
46: temp : Equation_List;
47: begin
48: Construct(lh,temp);
49: Swap_Tail(eqi.last,temp);
50: eqi.last := Tail_Of(eqi.last);
51: end;
52: end if;
53: end Add_Hyperplane;
54:
55: procedure Change_Hyperplane ( i,j : in natural; h : in Vector ) is
56: begin
57: if rps = null
58: then return;
59: else declare
60: eqi : Equation_List := rps(i).first;
61: lv : Link_To_Vector;
62: count : natural := 1;
63: begin
64: while not Is_Null(eqi) loop
65: if count = j
66: then lv := Head_Of(eqi);
67: for k in h'range loop
68: lv(k) := h(k);
69: end loop;
70: return;
71: else count := count + 1;
72: eqi := Tail_Of(eqi);
73: end if;
74: end loop;
75: end;
76: end if;
77: end Change_Hyperplane;
78:
79: -----------------
80: -- SELECTORS --
81: -----------------
82:
83: function Number_of_Hyperplanes ( i : natural ) return natural is
84: begin
85: if rps = null
86: then return 0;
87: else return Length_Of(rps(i).first);
88: end if;
89: end Number_of_Hyperplanes;
90:
91: procedure Get_Hyperplane ( i,j : in natural; h : out Vector ) is
92: begin
93: h := (h'range => 0);
94: if rps = null
95: then return;
96: else declare
97: eqi : Equation_List := rps(i).first;
98: count : natural := 1;
99: begin
100: while not Is_Null(eqi) loop
101: if count = j
102: then h := Head_Of(eqi).all;
103: return;
104: else count := count + 1;
105: eqi := Tail_Of(eqi);
106: end if;
107: end loop;
108: end;
109: end if;
110: end Get_Hyperplane;
111:
112: procedure Polynomial_System ( n,nbfree : in natural ) is
113:
114: rndms : Standard_Complex_Vectors.Vector(0..nbfree);
115:
116: begin
117: -- GENERATE THE FREE COEFFICIENTS :
118: rndms(0) := Create(0.0);
119: for i in rndms'first+1..rndms'last loop
120: rndms(i) := Random1; -- random complex number with radius one
121: end loop;
122: -- BUILD THE RANDOM PRODUCT SYSTEM :
123: Random_Product_System.Init(n);
124: for i in 1..n loop
125: for j in 1..Number_of_Hyperplanes(i) loop
126: declare
127: ih : Standard_Natural_Vectors.Vector(0..n);
128: h : Standard_Complex_Vectors.Vector(0..n);
129: begin
130: Get_Hyperplane(i,j,ih);
131: for k in h'range loop
132: h(k) := rndms(ih(k));
133: end loop;
134: Random_Product_System.Add_Hyperplane(i,h);
135: end;
136: end loop;
137: end loop;
138: end Polynomial_System;
139:
140: function Verify ( n : natural; lp : List ) return natural is
141:
142: temp : List := lp;
143: stop : boolean := false;
144: matrix : array (1..n,1..n) of natural;
145: nb : natural;
146:
147: function Degenerate return boolean is
148: degen : boolean;
149: first : natural;
150: begin
151: for i in 1..n loop
152: first := matrix(i,1);
153: degen := true;
154: for j in 2..n loop
155: if matrix(i,j) /= first
156: then degen := false;
157: end if;
158: exit when not degen;
159: end loop;
160: if degen
161: then return true;
162: end if;
163: end loop;
164: for j in 1..n loop
165: first := matrix(1,j);
166: degen := true;
167: for i in 2..n loop
168: if matrix(i,j) /= first
169: then degen := false;
170: end if;
171: exit when not degen;
172: end loop;
173: if degen
174: then return true;
175: end if;
176: end loop;
177: return false;
178: end Degenerate;
179:
180: procedure PVerify ( i,n : in natural; sum : in out natural ) is
181: begin
182: if i > n
183: then if Is_Null(temp)
184: then sum := sum + 1;
185: stop := true;
186: elsif Degenerate
187: then stop := true;
188: else temp := Tail_Of(temp);
189: sum := sum + 1;
190: end if;
191: else declare
192: eqi : Equation_List := rps(i).first;
193: h : Vector(0..n);
194: count : natural := 0;
195: begin
196: while not Is_Null(eqi) loop
197: count := count + 1;
198: if count = Head_Of(temp)(i)
199: then h := Head_Of(eqi).all;
200: for j in 1..n loop
201: matrix(i,j) := h(j);
202: end loop;
203: PVerify(i+1,n,sum);
204: end if;
205: exit when stop;
206: eqi := Tail_Of(eqi);
207: end loop;
208: end;
209: end if;
210: end PVerify;
211:
212: begin
213: nb := 0;
214: if not Is_Null(temp)
215: then PVerify(1,n,nb);
216: end if;
217: return nb;
218: end Verify;
219:
220: ------------------
221: -- DESTRUCTOR --
222: ------------------
223:
224: procedure Clear ( eql : in out Equation_List ) is
225:
226: temp : Equation_List := eql;
227: lv : Link_To_Vector;
228:
229: begin
230: while not Is_Null(temp) loop
231: lv := Head_Of(temp);
232: Clear(lv);
233: temp := Tail_of(temp);
234: end loop;
235: List_Of_Vectors.Clear(List_Of_Vectors.List(eql));
236: end Clear;
237:
238: procedure Clear ( eq : in out Equation ) is
239: begin
240: Clear(eq.first);
241: -- eq.last is just a pointer to the last element of eq.first;
242: -- if eq.first disappears, then also eq.last does
243: end Clear;
244:
245: procedure Clear ( eqs : in out Equations ) is
246: begin
247: for i in eqs'range loop
248: Clear(eqs(i));
249: end loop;
250: end Clear;
251:
252: procedure Clear is
253: begin
254: if rps /= null
255: then for i in rps'range loop
256: Clear(rps(i));
257: end loop;
258: free(rps);
259: end if;
260: end Clear;
261:
262: end Templates;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>