[0f5091] | 1 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 2 | |
---|
[1ce62fa] | 3 | version="$Id: solve.lib,v 1.6 1999-07-02 14:33:45 wenk Exp $"; |
---|
[0f5091] | 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 |
---|
[1ce62fa] | 11 | interpolate(i,j,d); interpolate poly from evaluation points i and results j |
---|
[0f5091] | 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) |
---|
[0d69875] | 20 | l>0: defines precision of fractional part if ground field is Q |
---|
[0f5091] | 21 | m=0,1,2: number of iterations for approximation of roots (default=2) |
---|
[0d69875] | 22 | ASSUME: i is a zerodimensional ideal with |
---|
| 23 | nvars(basering) = ncols(i) = number of vars actually occuring in i |
---|
[0f5091] | 24 | RETURN: list of all (complex) roots of the polynomial system i = 0, |
---|
[0d69875] | 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 |
---|
[0f5091] | 27 | EXAMPLE: example ures_solve; shows an example |
---|
| 28 | " |
---|
| 29 | { |
---|
[0d69875] | 30 | int typ=0; |
---|
[0f5091] | 31 | int polish=2; |
---|
| 32 | int prec=30; |
---|
| 33 | |
---|
[0d69875] | 34 | if ( size(#) >= 1 ) |
---|
[0f5091] | 35 | { |
---|
[0d69875] | 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 | } |
---|
[0f5091] | 43 | } |
---|
[0d69875] | 44 | if ( size(#) >= 2 ) |
---|
[0f5091] | 45 | { |
---|
[0d69875] | 46 | prec= #[2]; |
---|
| 47 | if ( prec == 0 ) { prec = 30; } |
---|
| 48 | if ( prec < 0 ) |
---|
| 49 | { |
---|
| 50 | ERROR("Third parameter l must be positive!"); |
---|
| 51 | } |
---|
[0f5091] | 52 | } |
---|
[0d69875] | 53 | if ( size(#) >= 3 ) |
---|
[0f5091] | 54 | { |
---|
[0d69875] | 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 | } |
---|
[0f5091] | 61 | } |
---|
[6035d5] | 62 | |
---|
[0d69875] | 63 | if ( size(#) > 3 ) |
---|
| 64 | { |
---|
| 65 | ERROR("only three parameters allowed!"); |
---|
| 66 | } |
---|
[0f5091] | 67 | |
---|
| 68 | int digits= system("setFloatDigits",prec); |
---|
| 69 | |
---|
| 70 | return(uressolve(gls,typ,polish)); |
---|
[0d69875] | 71 | |
---|
[0f5091] | 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; |
---|
[0d69875] | 83 | // now with complex coefficient field, precision is 10 digits |
---|
[0f5091] | 84 | ring rsc= (real,10,I),(x,y),lp; |
---|
[0d69875] | 85 | ideal i = x2 + y2 - 8, x2 + xy + 2y2; |
---|
[0f5091] | 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 |
---|
[0d69875] | 93 | l>0: defines precision of fractional part if ground field is Q |
---|
[0f5091] | 94 | m=0,1,2: number of iterations for approximation of roots (default=2) |
---|
| 95 | ASSUME: p is an univariate polynom |
---|
[7a4087] | 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 |
---|
[0f5091] | 99 | EXAMPLE: example laguerre_solve; shows an example |
---|
| 100 | " |
---|
| 101 | { |
---|
[0d69875] | 102 | int polish=2; |
---|
[0f5091] | 103 | int prec=30; |
---|
| 104 | |
---|
[0d69875] | 105 | if ( size(#) >= 1 ) |
---|
[0f5091] | 106 | { |
---|
[0d69875] | 107 | prec= #[1]; |
---|
| 108 | if ( prec == 0 ) { prec = 30; } |
---|
| 109 | if ( prec < 0 ) |
---|
| 110 | { |
---|
| 111 | ERROR("Fisrt parameter must be positive!"); |
---|
| 112 | } |
---|
[0f5091] | 113 | } |
---|
[0d69875] | 114 | if ( size(#) >= 2 ) |
---|
[0f5091] | 115 | { |
---|
[0d69875] | 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!"); |
---|
[0f5091] | 126 | } |
---|
| 127 | |
---|
| 128 | int digits= system("setFloatDigits",prec); |
---|
| 129 | |
---|
| 130 | return(laguerre(f,polish)); |
---|
[0d69875] | 131 | |
---|
[0f5091] | 132 | } |
---|
| 133 | example |
---|
| 134 | { |
---|
| 135 | "EXAMPLE:";echo=2; |
---|
| 136 | // Find all roots of an univariate polynomial using Laguerre's method: |
---|
[d866a4] | 137 | ring rs1= 0,(x,y),lp; |
---|
| 138 | poly f = 15x5 + x3 + x2 - 10; |
---|
[0f5091] | 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: |
---|
[d866a4] | 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; |
---|
[0f5091] | 151 | } |
---|
| 152 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 153 | |
---|
[0d69875] | 154 | proc mp_res_mat( ideal i, list # ) |
---|
[0f5091] | 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, |
---|
[0d69875] | 159 | for k=1 i must be homogeneous |
---|
| 160 | RETURN: module representing the multipolynomial resultant matrix |
---|
[0f5091] | 161 | EXAMPLE: example mp_res_mat; shows an example |
---|
| 162 | " |
---|
| 163 | { |
---|
| 164 | int typ=2; |
---|
| 165 | |
---|
[0d69875] | 166 | if ( size(#) == 1 ) |
---|
[0f5091] | 167 | { |
---|
[0d69875] | 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 | } |
---|
[0f5091] | 175 | } |
---|
[0d69875] | 176 | if ( size(#) > 1 ) |
---|
| 177 | { |
---|
| 178 | ERROR("only two parameters allowed!"); |
---|
| 179 | } |
---|
| 180 | |
---|
| 181 | return(mpresmat(i,typ)); |
---|
| 182 | |
---|
[0f5091] | 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; |
---|
[0d69875] | 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); |
---|
[0f5091] | 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; |
---|
[0d69875] | 198 | ideal h= homog(imap(rsq,i),x0); |
---|
[0f5091] | 199 | hgls; |
---|
[0d69875] | 200 | module m = mp_res_mat(h,1); |
---|
[0f5091] | 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; |
---|
[0d69875] | 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); |
---|
[0f5091] | 210 | print(mn); |
---|
| 211 | // computing sparse resultant |
---|
| 212 | det(mn); |
---|
| 213 | } |
---|
| 214 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 215 | |
---|
[0d69875] | 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, |
---|
[0f5091] | 219 | p and v consist of numbers of the ground filed K, p must have |
---|
[0d69875] | 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]) |
---|
[0f5091] | 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 |
---|
[0d69875] | 225 | of degree 0,1,...,N at p in lexicographical order, |
---|
[0f5091] | 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 | { |
---|
[0d69875] | 232 | return(vandermonde(p,w,d)); |
---|
[0f5091] | 233 | } |
---|
| 234 | example |
---|
| 235 | { |
---|
| 236 | "EXAMPLE:"; echo=2; |
---|
| 237 | ring r1 = 0,(x),lp; |
---|
[d866a4] | 238 | // determine f with deg(f) = 4 and |
---|
[0f5091] | 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; |
---|
[0d69875] | 241 | interpolate( 3, v, 4 ); |
---|
[0f5091] | 242 | |
---|
| 243 | ring r2 = 0,(x,y),dp; |
---|
[d866a4] | 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) |
---|
[0f5091] | 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; |
---|
[0d69875] | 249 | interpolate( p,v,3 ); |
---|
[0f5091] | 250 | } |
---|
| 251 | /////////////////////////////////////////////////////////////////////////////// |
---|
[0d69875] | 252 | |
---|
| 253 | |
---|
| 254 | |
---|
| 255 | |
---|
| 256 | |
---|
| 257 | |
---|