Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Product/interpolating_homotopies_driver.adb, Revision 1.1.1.1
1.1 maekawa 1: with integer_io; use integer_io;
2: with Communications_with_User; use Communications_with_User;
3: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
4: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
5: with Standard_Complex_Numbers_io; use Standard_Complex_Numbers_io;
6: with Numbers_io; use Numbers_io;
7: with Standard_Random_Numbers; use Standard_Random_Numbers;
8: with Standard_Natural_Vectors;
9: with Standard_Complex_Vectors; use Standard_Complex_Vectors;
10: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
11: with Interpolating_Homotopies; use Interpolating_Homotopies;
12:
13: procedure Interpolating_Homotopies_Driver
14: ( file : in file_type; p : in Poly_Sys; z : in Partition;
15: b : in out natural; q : out Poly_Sys;
16: qsols : in out Solution_List ) is
17:
18: n : constant natural := p'last;
19: interpols : Solution_List;
20: ib,scalind : natural;
21: ans : character;
22:
23: function Random_Interpolating ( n,m : natural ) return Solution_List is
24:
25: -- DESCRIPTION :
26: -- A list of m random n-dimensional vectors will be returned.
27: -- The complex numbers will all have modulus one.
28:
29: res,res_last : Solution_List;
30: s : Solution(n);
31:
32: begin
33: s.t := Create(0.0);
34: s.m := 1;
35: s.err := 0.0; s.rco := 1.0; s.res := 0.0;
36: for i in 1..m loop
37: for j in 1..n loop
38: s.v(j) := random1;
39: end loop;
40: Append(res,res_last,s);
41: end loop;
42: return res;
43: end Random_Interpolating;
44:
45: procedure Random_Linear_Scaler ( n : in natural; p : in Poly;
46: v : out vector; l : out natural ) is
47:
48: -- DESCRIPTION :
49: -- Returns a random vector of dimension n+1, with range 0..n.
50: -- There will be a nonzero entry only for those unknowns that occur in p.
51:
52: res : Vector(0..n);
53: last : natural := 0;
54:
55: begin
56: for i in res'range loop
57: if Degree(p,i) > 0
58: then res(i) := random1;
59: last := last + 1;
60: else res(i) := Create(0.0);
61: end if;
62: end loop;
63: v := res; l := last;
64: end Random_Linear_Scaler;
65:
66: function Scale ( sc,v : Vector; last : integer ) return Complex_Number is
67:
68: -- DESCRIPTION :
69: -- Returns the last component of the vector v, that is v(last),
70: -- such that sc(0) + sum sc(i)*v(i), i in v'range, holds.
71:
72: res : Complex_Number := sc(0);
73:
74: begin
75: for i in v'first..last-1 loop
76: res := res + sc(i)*v(i);
77: end loop;
78: res := -res/sc(last);
79: return res;
80: end Scale;
81:
82: function Random_Interpolating
83: ( n,m : natural; scaler : Vector; scallast : natural )
84: return Solution_List is
85:
86: -- DESCRIPTION :
87: -- A list of m random n-dimensional vectors will be returned.
88: -- The complex numbers will all have modulus one, except the last one,
89: -- indicated by scallast, that has been chosen to satisfy the scaler
90: -- equation, defined by sum of scaler(i)*x_i = 0, with x_0 = 1.
91:
92: res,res_last : Solution_List;
93: s : Solution(n);
94:
95: begin
96: s.t := Create(0.0);
97: s.m := 1;
98: for i in 1..m loop
99: for j in 1..(n-1) loop
100: s.v(j) := random1;
101: end loop;
102: s.v(n) := Scale(scaler,s.v,scallast);
103: Append(res,res_last,s);
104: end loop;
105: return res;
106: end Random_Interpolating;
107:
108: function Create ( v : Vector ) return Poly is
109:
110: res : Poly := Null_Poly;
111: t : Term;
112:
113: begin
114: t.dg := new Standard_Natural_Vectors.Vector'(v'first+1..v'last => 0);
115: for i in v'range loop
116: t.cf := v(i);
117: if i > v'first
118: then t.dg(i) := 1;
119: end if;
120: Add(res,t);
121: if i > v'first
122: then t.dg(i) := 0;
123: end if;
124: end loop;
125: Clear(t);
126: return res;
127: end Create;
128:
129: function Interpolating_by_User ( n,m : natural ) return Solution_List is
130:
131: -- DESCRIPTION :
132: -- A list of m n-dimensional vectors will be read from standard input.
133:
134: res,res_last : Solution_List;
135: s : Solution(n);
136: -- f1,f2 : double_float;
137:
138: begin
139: put("Reading "); put(m,1); put(" "); put(n,1);
140: put_line("-dimensional complex vectors.");
141: for i in 1..m loop
142: s.t := Create(0.0);
143: s.m := 1;
144: s.err := 0.0; s.rco := 1.0; s.res := 0.0;
145: put("Give the components of vector "); put(i,1);
146: put_line(" :");
147: for j in 1..n loop
148: -- Read_Double_Float(f1);
149: -- Read_Double_Float(f2);
150: -- s.v(j) := Create(f1,f2);
151: get(s.v(j));
152: end loop;
153: Append(res,res_last,s);
154: end loop;
155: return res;
156: end Interpolating_by_User;
157:
158: procedure Driver_for_Interpolation is
159:
160: -- DESCRIPTION : interpolation without a scaling equation
161:
162: dp,ip : Poly_Sys(p'range);
163:
164: begin
165: dp := Dense_Representation(p,z);
166: ip := Independent_Representation(dp);
167: ib := Independent_Roots(ip);
168: if ib > b
169: then ib := b;
170: end if;
171: put("The number of independent roots : "); put(ib,1); new_line;
172: put("Do you want to give interpolation vectors by yourself ? (y/n) ");
173: Ask_Yes_or_No(ans);
174: if ans = 'y'
175: then interpols := Interpolating_by_User(n,ib);
176: else put_line("Random interpolating vectors will be generated.");
177: interpols := Random_Interpolating(n,ib);
178: end if;
179: q := Interpolate(ip,ib,interpols);
180: qsols := interpols; b := ib;
181: Clear(dp); Clear(ip);
182: end Driver_for_Interpolation;
183:
184: procedure Driver_for_Scaled_Interpolation is
185:
186: -- DESCRIPTION : interpolation with a scaling equation, p(scalind).
187:
188: dp,ip,pp,qq : Poly_Sys(p'range);
189: scalvec : Vector(0..n);
190: scalveclast : natural;
191:
192: begin
193: Random_Linear_Scaler(n,p(scalind),scalvec,scalveclast);
194: for i in p'range loop
195: if i = scalind
196: then pp(i) := Null_Poly;
197: else pp(i) := p(i);
198: end if;
199: end loop;
200: dp := Dense_Representation(pp,z);
201: ip := Independent_Representation(dp);
202: ib := Independent_Roots(ip,scalveclast);
203: if ib > b
204: then ib := b;
205: end if;
206: put("The number of independent roots : "); put(ib,1); new_line;
207: put("Do you want to give interpolation vectors by yourself ? (y/n) ");
208: Ask_Yes_or_No(ans);
209: if ans = 'y'
210: then interpols := Interpolating_by_User(n,ib);
211: else put_line("Random interpolating vectors will be generated.");
212: interpols := Random_Interpolating(n,ib,scalvec,scalveclast);
213: end if;
214: qq := Interpolate(ip,scalveclast,ib,interpols);
215: qq(scalind) := Create(scalvec);
216: qsols := interpols; b := ib; q := qq;
217: Clear(dp); Clear(ip);
218: end Driver_for_Scaled_Interpolation;
219:
220: begin
221: new_line;
222: put("Give the number of the linear scaling equation (0 if none) : ");
223: Read_Natural(scalind);
224: if scalind = 0
225: then Driver_for_Interpolation;
226: else Driver_for_Scaled_Interpolation;
227: end if;
228: end Interpolating_Homotopies_Driver;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to interpolating_homotopies_driver.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