Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Product/multi_homogeneous_start_systems.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
2: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
3: with Standard_Random_Numbers; use Standard_Random_Numbers;
4: with Standard_Natural_Vectors; use Standard_Natural_Vectors;
5: with Standard_Complex_Vectors;
6: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
7: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
8: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
9:
10: with Sets_of_Unknowns; use Sets_of_Unknowns;
11: with Degrees_in_Sets_of_Unknowns; use Degrees_in_Sets_of_Unknowns;
12: with Partitions_of_Sets_of_Unknowns; use Partitions_of_Sets_of_Unknowns;
13: with Degree_Structure; use Degree_Structure;
14: with Random_Product_System;
15:
16: package body Multi_Homogeneous_Start_Systems is
17:
18: procedure Add_Hyperplanes ( i : in natural; p : in Poly;
19: z : in partition; m : in natural;
20: dg : in Standard_Natural_Vectors.Vector) is
21: -- DESCRIPTION :
22: -- This routine adds hyperplanes to the Random Product System
23: -- according to the partition and the degrees of the polynomial
24:
25: -- ON ENTRY :
26: -- i the number of the polynomial in the system;
27: -- p the i-the polynomial of the system;
28: -- z a partition of the set of unknowns of p;
29: -- m the number of sets in z;
30: -- dg the degrees of the polynomial in the sets of z,
31: -- for j in 1..m : dg(j) = Degree(p,z(j)).
32:
33: n : natural := Number_Of_Unknowns(p);
34: h : Standard_Complex_Vectors.Vector(0..n);
35:
36: begin
37: for k in 1..m loop
38: for l in 1..dg(k) loop
39: h(0) := random1;
40: for j in 1..n loop
41: if Is_In(z(k),j)
42: then h(j) := Random1;
43: else h(j) := Create(0.0);
44: end if;
45: end loop;
46: Random_Product_System.Add_Hyperplane(i,h);
47: end loop;
48: end loop;
49: end Add_Hyperplanes;
50:
51: procedure Construct_Random_Product_System ( p : in Poly_Sys ) is
52:
53: -- DESCRIPTION :
54: -- This procedure constructs a random product system by
55: -- finding a good partition for each equation of the system p.
56:
57: n : natural := p'length;
58: m : natural;
59: begin
60: if Degree_Structure.Empty
61: then Degree_Structure.Create(p);
62: end if;
63: for i in p'range loop
64: m := Degree_Structure.Get(i);
65: declare
66: z : Partition(1..m);
67: dg : Standard_Natural_Vectors.Vector(1..m);
68: begin
69: Degree_Structure.Get(i,z,dg);
70: Add_Hyperplanes(i,p(i),z,m,dg);
71: Clear(z);
72: end;
73: end loop;
74: end Construct_Random_Product_System;
75:
76: procedure RPQ ( p : in Poly_Sys; q : out Poly_Sys;
77: sols : in out Solution_List; nl : out natural) is
78:
79: n : natural := p'length;
80:
81: begin
82: Random_Product_System.Init(n);
83: Construct_Random_Product_System(p);
84: q := Random_Product_System.Polynomial_System;
85: Random_Product_System.Solve(sols,nl);
86: Degree_Structure.Clear;
87: Random_Product_System.Clear;
88: end RPQ;
89:
90: procedure Solve ( i,n : in natural; sols : in out Solution_List;
91: a : in out Matrix;
92: b : in out Standard_Complex_Vectors.Vector;
93: acc : in out partition ) is
94: begin
95: if i > n
96: then declare
97: aa : Matrix(a'range(1),a'range(2));
98: bb : Standard_Complex_Vectors.Vector(b'range);
99: rcond : double_float;
100: ipvt : Standard_Natural_Vectors.Vector(b'range);
101: begin
102: for k in a'range(1) loop
103: for l in a'range(2) loop
104: aa(k,l) := a(k,l);
105: end loop;
106: bb(k) := b(k);
107: end loop;
108: lufco(aa,n,ipvt,rcond);
109: -- put("rcond : "); put(rcond); new_line;
110: if rcond + 1.0 /= 1.0
111: then lusolve(aa,n,ipvt,bb);
112: declare
113: s : Solution(n);
114: begin
115: s.t := Create(0.0);
116: s.m := 1;
117: s.v := bb;
118: Add(sols,s);
119: end;
120: end if;
121: exception
122: when numeric_error => return;
123: end;
124: else declare
125: h : Standard_Complex_Vectors.Vector(0..n);
126: count : natural := 0;
127: z : Partition(1..n);
128: m : natural;
129: d : Standard_Natural_Vectors.Vector(1..n);
130: begin
131: Degree_Structure.Get(i,z,d);
132: m := Degree_Structure.Get(i);
133: for j1 in 1..m loop
134: if Degree_Structure.Admissible(acc,i-1,z(j1))
135: then acc(i) := Create(z(j1));
136: for j2 in 1..d(j1) loop
137: count := count + 1;
138: h := Random_Product_System.Get_Hyperplane(i,count);
139: b(i) := -h(0);
140: for k in 1..n loop
141: a(i,k) := h(k);
142: end loop;
143: Solve(i+1,n,sols,a,b,acc);
144: end loop;
145: Clear(acc(i));
146: else count := count + d(j1);
147: end if;
148: end loop;
149: Clear(z);
150: end;
151: end if;
152: end Solve;
153:
154: procedure Solve_Start_System
155: ( n : in natural; sols : in out Solution_List ) is
156:
157: m : Matrix(1..n,1..n);
158: v : Standard_Complex_Vectors.Vector(1..n);
159: acc : Partition(1..n);
160:
161: begin
162: for i in 1..n loop
163: for j in 1..n loop
164: m(i,j) := Create(0.0);
165: end loop;
166: v(i) := Create(0.0);
167: end loop;
168: Solve(1,n,sols,m,v,acc);
169: end Solve_Start_System;
170:
171: procedure GBQ ( p : in Poly_Sys; q : out Poly_Sys;
172: sols : in out Solution_List ) is
173:
174: n : natural := p'length;
175:
176: begin
177: Random_Product_System.Init(n);
178: Construct_Random_Product_System(p);
179: q := Random_Product_System.Polynomial_System;
180: Solve_Start_System(n,sols);
181: Degree_Structure.Clear;
182: Random_Product_System.Clear;
183: end GBQ;
184:
185: end Multi_Homogeneous_Start_Systems;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to multi_homogeneous_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