Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Product/random_product_start_systems.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
2: with Standard_Random_Numbers; use Standard_Random_Numbers;
3: with Standard_Natural_Vectors;
4: with Standard_Complex_Vectors; use Standard_Complex_Vectors;
5: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
6: with Random_Product_System;
7: with Set_Structure;
8:
9: package body Random_Product_Start_Systems is
10:
11: type Boolean_Array is array(integer range <>) of boolean;
12:
13: function leq ( d1,d2 : Degrees ) return boolean is
14:
15: -- DESCRIPTION :
16: -- returns true if for all k: d1(k) <= d2(k)
17:
18: begin
19: for k in d1'range loop
20: if d1(k) > d2(k)
21: then return false;
22: end if;
23: end loop;
24: return true;
25: end leq;
26:
27: function Dominated ( d : Degrees; p : Poly ) return boolean is
28:
29: -- DESCRIPTION :
30: -- Returns true if the term indicated by the degrees d
31: -- is already in the set structure;
32: -- the polynomial p contains all terms that are already
33: -- belonging to the set structure.
34:
35: res : boolean := false;
36:
37: procedure Scan ( t : in Term; continue : out boolean ) is
38: begin
39: if leq(d,t.dg)
40: then res := true;
41: continue := false;
42: else continue := true;
43: end if;
44: end Scan;
45: procedure Scan_Terms is new Visiting_Iterator (Scan);
46:
47: begin
48: Scan_Terms(p);
49: return res;
50: end Dominated;
51:
52: procedure Divide ( dg : in out Standard_Natural_Vectors.Vector;
53: i,d : in natural; occupied : in out Boolean_Array ) is
54:
55: -- DESCRIPTION :
56: -- Divides the monomial by every unknown that occurs already
57: -- in one set of the ith component of the set structure.
58: -- Every set that contains already one unknown of the monomial is
59: -- marked by setting the corresponding entry in occupied on true.
60:
61: j : natural;
62:
63: begin
64: for k in dg'range loop
65: j := 1;
66: while (dg(k) /= 0) and (j <= d) loop
67: if not occupied(j) and Set_Structure.Is_In(i,j,k)
68: then occupied(j) := true;
69: dg(k) := dg(k)-1;
70: end if;
71: j := j+1;
72: end loop;
73: end loop;
74: end Divide;
75:
76: procedure Augment ( ba : in Boolean_Array; index : in out integer ) is
77:
78: -- DESCRIPTION :
79: -- Augments the index till ba(index) = false.
80:
81: begin
82: while ba(index) loop
83: index := index + 1;
84: exit when index >= ba'last;
85: end loop;
86: end Augment;
87:
88: procedure Build_Set_Structure ( i,di : in natural; p : in Poly ) is
89:
90: -- DESCRIPTION :
91: -- This procedure builds the set structure for the polynomial p
92:
93: -- ON ENTRY :
94: -- i the number of the equation
95: -- di the degree of the polynomial
96: -- p the polynomial for the i-th equation
97:
98: temp : Poly := Null_Poly;
99:
100: procedure Process ( t : in Term; continue : out boolean ) is
101:
102: ind : natural := 1; -- number of the set
103: processed : Boolean_Array(1..di) := (1..di => false);
104: tdg : Standard_Natural_Vectors.Vector(t.dg'range) := t.dg.all;
105:
106: begin
107: if not Dominated(t.dg,temp)
108: then Divide(tdg,i,di,processed);
109: Augment(processed,ind);
110: for k in tdg'range loop
111: for j in 1..tdg(k) loop
112: if not Set_Structure.Is_In(i,ind,k)
113: then Set_Structure.Add(i,ind,k);
114: processed(ind) := true;
115: end if;
116: Augment(processed,ind);
117: end loop;
118: end loop;
119: Add(temp,t);
120: end if;
121: continue := true;
122: end Process;
123: procedure Process_Terms is new Visiting_Iterator (Process);
124:
125: begin
126: Process_Terms(p);
127: Clear(temp);
128: end Build_Set_Structure;
129:
130: procedure Build_Set_Structure ( p : in Poly_Sys ) is
131:
132: n : natural := p'length;
133: ns : Standard_Natural_Vectors.Vector(1..n);
134:
135: begin
136: for i in p'range loop
137: ns(i) := Degree(p(i));
138: end loop;
139: Set_Structure.Init(ns);
140: for i in p'range loop
141: Build_Set_Structure(i,Degree(p(i)),p(i));
142: end loop;
143: end Build_Set_Structure;
144:
145: procedure Build_Random_Product_System ( n : in natural ) is
146:
147: -- DESCRIPTION :
148: -- Based upon the set structure of p, a random product system
149: -- will be constructed.
150:
151: h : Vector(0..n);
152: first : boolean;
153:
154: begin
155: for i in 1..n loop
156: for j in 1..Set_Structure.Number_Of_Sets(i) loop
157: h := (0..n => Create(0.0));
158: first := true;
159: for k in 1..n loop
160: if Set_Structure.Is_In(i,j,k)
161: then if first
162: then h(k) := Create(1.0); first := false;
163: else h(k) := Random1;
164: end if;
165: end if;
166: end loop;
167: h(0) := Random1;
168: Random_Product_System.Add_Hyperplane(i,h);
169: end loop;
170: end loop;
171: end Build_Random_Product_System;
172:
173: procedure Construct ( p : in Poly_Sys; q : in out Poly_Sys;
174: sols : in out Solution_List ) is
175:
176: nl : natural := 0;
177: n : natural := p'length;
178:
179: begin
180: Build_Set_Structure(p);
181: Random_Product_System.Init(n);
182: Build_Random_Product_System(n);
183: q := Random_Product_System.Polynomial_System;
184: Random_Product_System.Solve(sols,nl);
185: Random_Product_System.Clear;
186: end Construct;
187:
188: end Random_Product_Start_Systems;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to random_product_start_systems.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