1 | /////////////////////////////////////////////////////////////////////////////// |
---|
2 | |
---|
3 | version="$Id: solve.lib,v 1.6 1999-07-02 14:33:45 wenk Exp $"; |
---|
4 | info=" |
---|
5 | LIBRARY: solve.lib PROCEDURES TO SOLVE POLYNOMIAL SYSTEMS |
---|
6 | AUTHOR: Moritz Wenk, email: wenk@mathematik.uni-kl.de |
---|
7 | |
---|
8 | ures_solve(i,..); find all roots of 0-dimensional ideal i with resultants |
---|
9 | mp_res_mat(i,..); multipolynomial resultant matrix of ideal i |
---|
10 | laguerre_solve(p,..); find all roots of univariate polynom p |
---|
11 | interpolate(i,j,d); interpolate poly from evaluation points i and results j |
---|
12 | "; |
---|
13 | |
---|
14 | /////////////////////////////////////////////////////////////////////////////// |
---|
15 | |
---|
16 | proc ures_solve( ideal gls, list # ) |
---|
17 | "USAGE: ures_solve(i[,k,l,m]); i ideal, k,l,m integers |
---|
18 | k=0: use sparse resultant matrix of Gelfand, Kapranov and Zelevinsky |
---|
19 | k=1: use resultant matrix of Macaulay (k=0 is default) |
---|
20 | l>0: defines precision of fractional part if ground field is Q |
---|
21 | m=0,1,2: number of iterations for approximation of roots (default=2) |
---|
22 | ASSUME: i is a zerodimensional ideal with |
---|
23 | nvars(basering) = ncols(i) = number of vars actually occuring in i |
---|
24 | RETURN: list of all (complex) roots of the polynomial system i = 0, |
---|
25 | of type number if the ground field is the complex numbers, |
---|
26 | of type string if the ground field is the rational or real numbers |
---|
27 | EXAMPLE: example ures_solve; shows an example |
---|
28 | " |
---|
29 | { |
---|
30 | int typ=0; |
---|
31 | int polish=2; |
---|
32 | int prec=30; |
---|
33 | |
---|
34 | if ( size(#) >= 1 ) |
---|
35 | { |
---|
36 | typ= #[1]; |
---|
37 | if ( typ < 0 || typ > 1 ) |
---|
38 | { |
---|
39 | ERROR("Valid values for second parameter k are: |
---|
40 | 0: use sparse Resultant (default) |
---|
41 | 1: use Macaulay Resultant"); |
---|
42 | } |
---|
43 | } |
---|
44 | if ( size(#) >= 2 ) |
---|
45 | { |
---|
46 | prec= #[2]; |
---|
47 | if ( prec == 0 ) { prec = 30; } |
---|
48 | if ( prec < 0 ) |
---|
49 | { |
---|
50 | ERROR("Third parameter l must be positive!"); |
---|
51 | } |
---|
52 | } |
---|
53 | if ( size(#) >= 3 ) |
---|
54 | { |
---|
55 | polish= #[3]; |
---|
56 | if ( polish < 0 || polish > 2 ) |
---|
57 | { |
---|
58 | ERROR("Valid values for fourth parameter m are: |
---|
59 | 0,1,2: number of iterations for approximation of roots"); |
---|
60 | } |
---|
61 | } |
---|
62 | |
---|
63 | if ( size(#) > 3 ) |
---|
64 | { |
---|
65 | ERROR("only three parameters allowed!"); |
---|
66 | } |
---|
67 | |
---|
68 | int digits= system("setFloatDigits",prec); |
---|
69 | |
---|
70 | return(uressolve(gls,typ,polish)); |
---|
71 | |
---|
72 | } |
---|
73 | example |
---|
74 | { |
---|
75 | "EXAMPLE:";echo=2; |
---|
76 | // compute the intersection points of two curves |
---|
77 | ring rsq = 0,(x,y),lp; |
---|
78 | ideal gls= x2 + y2 - 10, x2 + xy + 2y2 - 16; |
---|
79 | ures_solve(gls); |
---|
80 | // result is a list (x,y)-coordinates as strings |
---|
81 | |
---|
82 | pause; |
---|
83 | // now with complex coefficient field, precision is 10 digits |
---|
84 | ring rsc= (real,10,I),(x,y),lp; |
---|
85 | ideal i = x2 + y2 - 8, x2 + xy + 2y2; |
---|
86 | ures_solve(i); |
---|
87 | // result is a list of (x,y)-coordinates of complex numbers |
---|
88 | } |
---|
89 | /////////////////////////////////////////////////////////////////////////////// |
---|
90 | |
---|
91 | proc laguerre_solve( poly f, list # ) |
---|
92 | "USAGE: laguerre_solve( p[,l,m]); f poly, l,m integers |
---|
93 | l>0: defines precision of fractional part if ground field is Q |
---|
94 | m=0,1,2: number of iterations for approximation of roots (default=2) |
---|
95 | ASSUME: p is an univariate polynom |
---|
96 | RETURN: list of all (complex) roots of the polynomial p; |
---|
97 | of type number if the ground field is the complex numbers, |
---|
98 | of type string if the ground field is the rational or real numbers |
---|
99 | EXAMPLE: example laguerre_solve; shows an example |
---|
100 | " |
---|
101 | { |
---|
102 | int polish=2; |
---|
103 | int prec=30; |
---|
104 | |
---|
105 | if ( size(#) >= 1 ) |
---|
106 | { |
---|
107 | prec= #[1]; |
---|
108 | if ( prec == 0 ) { prec = 30; } |
---|
109 | if ( prec < 0 ) |
---|
110 | { |
---|
111 | ERROR("Fisrt parameter must be positive!"); |
---|
112 | } |
---|
113 | } |
---|
114 | if ( size(#) >= 2 ) |
---|
115 | { |
---|
116 | polish= #[2]; |
---|
117 | if ( polish < 0 || polish > 2 ) |
---|
118 | { |
---|
119 | ERROR("Valid values for third parameter are: |
---|
120 | 0,1,2: number of iterations for approximation of roots"); |
---|
121 | } |
---|
122 | } |
---|
123 | if ( size(#) > 2 ) |
---|
124 | { |
---|
125 | ERROR("only two parameters allowed!"); |
---|
126 | } |
---|
127 | |
---|
128 | int digits= system("setFloatDigits",prec); |
---|
129 | |
---|
130 | return(laguerre(f,polish)); |
---|
131 | |
---|
132 | } |
---|
133 | example |
---|
134 | { |
---|
135 | "EXAMPLE:";echo=2; |
---|
136 | // Find all roots of an univariate polynomial using Laguerre's method: |
---|
137 | ring rs1= 0,(x,y),lp; |
---|
138 | poly f = 15x5 + x3 + x2 - 10; |
---|
139 | laguerre_solve(f); |
---|
140 | |
---|
141 | pause; |
---|
142 | // Now with 10 digits precision: |
---|
143 | laguerre_solve(f,10); |
---|
144 | |
---|
145 | pause; |
---|
146 | // Now with complex coefficients, precision is 20 digits: |
---|
147 | ring rsc= (real,20,I),x,lp; |
---|
148 | poly f = (15+I*5)*x^5 + (0.25+I*2)*x^3 + x2 - 10*I; |
---|
149 | list l = laguerre_solve(f); |
---|
150 | l; |
---|
151 | } |
---|
152 | /////////////////////////////////////////////////////////////////////////////// |
---|
153 | |
---|
154 | proc mp_res_mat( ideal i, list # ) |
---|
155 | "USAGE: mp_res_mat(i[,k]); i ideal, k integer |
---|
156 | k=0: sparse resultant matrix of Gelfand, Kapranov and Zelevinsky |
---|
157 | k=1: resultant matrix of Macaulay (k=0 is default) |
---|
158 | ASSUME: nvars(basering) = ncols(i)-1 = number of vars actually occuring in i, |
---|
159 | for k=1 i must be homogeneous |
---|
160 | RETURN: module representing the multipolynomial resultant matrix |
---|
161 | EXAMPLE: example mp_res_mat; shows an example |
---|
162 | " |
---|
163 | { |
---|
164 | int typ=2; |
---|
165 | |
---|
166 | if ( size(#) == 1 ) |
---|
167 | { |
---|
168 | typ= #[1]; |
---|
169 | if ( typ < 0 || typ > 1 ) |
---|
170 | { |
---|
171 | ERROR("Valid values for third parameter are: |
---|
172 | 0: sparse resultant (default) |
---|
173 | 1: Macaulay resultant"); |
---|
174 | } |
---|
175 | } |
---|
176 | if ( size(#) > 1 ) |
---|
177 | { |
---|
178 | ERROR("only two parameters allowed!"); |
---|
179 | } |
---|
180 | |
---|
181 | return(mpresmat(i,typ)); |
---|
182 | |
---|
183 | } |
---|
184 | example |
---|
185 | { |
---|
186 | "EXAMPLE:";echo=2; |
---|
187 | // compute resultant matrix in ring with parameters (sparse resultant matrix) |
---|
188 | ring rsq= (0,u0,u1,u2),(x1,x2),lp; |
---|
189 | ideal i= u0+u1*x1+u2*x2,x1^2 + x2^2 - 10,x1^2 + x1*x2 + 2*x2^2 - 16; |
---|
190 | module m = mp_res_mat(i); |
---|
191 | print(m); |
---|
192 | // computing sparse resultant |
---|
193 | det(m); |
---|
194 | |
---|
195 | pause; |
---|
196 | // compute resultant matrix (Macaulay resultant matrix) |
---|
197 | ring rdq= (0,u0,u1,u2),(x0,x1,x2),lp; |
---|
198 | ideal h= homog(imap(rsq,i),x0); |
---|
199 | hgls; |
---|
200 | module m = mp_res_mat(h,1); |
---|
201 | print(m); |
---|
202 | // computing Macaulay resultant (should be the same as above!) |
---|
203 | det(m); |
---|
204 | |
---|
205 | pause; |
---|
206 | // compute numerical sparse resultant matrix |
---|
207 | setring rsq; |
---|
208 | ideal ir= 15+2*x1+5*x2,x1^2 + x2^2 - 10,x1^2 + x1*x2 + 2*x2^2 - 16; |
---|
209 | module mn = mp_res_mat(ir); |
---|
210 | print(mn); |
---|
211 | // computing sparse resultant |
---|
212 | det(mn); |
---|
213 | } |
---|
214 | /////////////////////////////////////////////////////////////////////////////// |
---|
215 | |
---|
216 | proc interpolate( ideal p, ideal w, int d ) |
---|
217 | "USAGE: interpolate(p,v,d); p,v ideals, d int |
---|
218 | ASSUME: ground field K is the rational numbers, |
---|
219 | p and v consist of numbers of the ground filed K, p must have |
---|
220 | n elements, v must have N=(d+1)^n elements where n=nvars(basering) |
---|
221 | and d=deg(f) (f is the unknown polynomial in K[x1,...,xn]) |
---|
222 | COMPUTE: polynomial f with values given by v at points p1,..,pN derived from p; |
---|
223 | more precisely: consider p as point in K^n and v as N elements in K, |
---|
224 | let p1,..,pN be the points in K^n obtained by evaluating all monomials |
---|
225 | of degree 0,1,...,N at p in lexicographical order, |
---|
226 | then the procedure computes the polynomial f satisfying f(pi) = v[i] |
---|
227 | RETURN: polynomial f of degree d |
---|
228 | NOTE: mainly useful for n=1, with f satisfying f(p^(i-1)) = v[i], i=1..d+1 |
---|
229 | EXAMPLE: example interpolate; shows an example |
---|
230 | " |
---|
231 | { |
---|
232 | return(vandermonde(p,w,d)); |
---|
233 | } |
---|
234 | example |
---|
235 | { |
---|
236 | "EXAMPLE:"; echo=2; |
---|
237 | ring r1 = 0,(x),lp; |
---|
238 | // determine f with deg(f) = 4 and |
---|
239 | // v = values of f at points 3^0, 3^1, 3^2, 3^3, 3^4 |
---|
240 | ideal v=16,0,11376,1046880,85949136; |
---|
241 | interpolate( 3, v, 4 ); |
---|
242 | |
---|
243 | ring r2 = 0,(x,y),dp; |
---|
244 | // determine f with deg(f) = 3 and |
---|
245 | // v = values of f at 16 points (2,3)^0=(1,1),...,(2,3)^15=(2^15,3^15) |
---|
246 | // valuation point (2,3) |
---|
247 | ideal p = 2,3; |
---|
248 | ideal v= 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16; |
---|
249 | interpolate( p,v,3 ); |
---|
250 | } |
---|
251 | /////////////////////////////////////////////////////////////////////////////// |
---|
252 | |
---|
253 | |
---|
254 | |
---|
255 | |
---|
256 | |
---|
257 | |
---|