[741c669] | 1 | // $Id: rootsmr.lib,v 1.4 2008-03-18 15:51:41 Singular Exp $ |
---|
[d3b3ab] | 2 | // E. Tobis 12.Nov.2004, April 2004 |
---|
[f5ba45] | 3 | // last change 7. May 2005 (G.-M. Greuel) |
---|
[e182c8] | 4 | /////////////////////////////////////////////////////////////////////////////// |
---|
[60d32e] | 5 | category="Teaching" |
---|
[e182c8] | 6 | info=" |
---|
[1c981c] | 7 | LIBRARY: rootsmr.lib Counting the number of real roots of polynomial systems |
---|
[d3b3ab] | 8 | AUTHOR: Enrique A. Tobis, etobis@dc.uba.ar |
---|
[e182c8] | 9 | |
---|
| 10 | OVERVIEW: Routines for counting the number of real roots of a multivariate |
---|
[d3b3ab] | 11 | polynomial system. Two methods are implemented: deterministic |
---|
| 12 | computation of the number of roots, via the signature of a certain |
---|
[f5ba45] | 13 | bilinear form (nrRootsDeterm); and a rational univariate projection, |
---|
| 14 | using a pseudorandom polynomial (nrRootsProbab). It also includes a |
---|
| 15 | command to verify the correctness of the pseudorandom answer. |
---|
[d3b3ab] | 16 | References: Basu, Pollack, Roy, \"Algorithms in Real Algebraic |
---|
| 17 | Geometry\", Springer, 2003. |
---|
[e182c8] | 18 | |
---|
| 19 | PROCEDURES: |
---|
[f5ba45] | 20 | nrRootsProbab(I) Number of real roots of 0-dim ideal (probabilistic) |
---|
| 21 | nrRootsDeterm(I) Number of real roots of 0-dim ideal (deterministic) |
---|
[d3b3ab] | 22 | symsignature(m) Signature of the symmetric matrix m |
---|
| 23 | sturmquery(h,B,I) Sturm query of h on V(I) |
---|
| 24 | matbil(h,B,I) Matrix of the bilinear form on R/I associated to h |
---|
| 25 | matmult(f,B,I) Matrix of multiplication by f (m_f) on R/I in the basis B |
---|
| 26 | tracemult(f,B,I) Trace of m_f (B is an ordered basis of R/I) |
---|
| 27 | coords(f,B,I) Coordinates of f in the ordered basis B |
---|
| 28 | randcharpoly(B,I,n) Pseudorandom charpoly of univ. projection, n optional |
---|
| 29 | verify(p,B,i) Verifies the result of randcharpoly |
---|
| 30 | randlinpoly(n) Pseudorandom linear polynomial, n optional |
---|
| 31 | powersums(f,B,I) Powersums of the roots of a char polynomial |
---|
| 32 | symmfunc(S) Symmetric functions from the powersums S |
---|
| 33 | univarpoly(l) Polynomial with coefficients from l |
---|
| 34 | qbase(i) Like kbase, but the monomials are ordered |
---|
[e182c8] | 35 | |
---|
| 36 | KEYWORDS: real roots, univariate projection |
---|
| 37 | "; |
---|
| 38 | /////////////////////////////////////////////////////////////////// |
---|
| 39 | LIB "linalg.lib"; // We use charpoly |
---|
[67f7108] | 40 | LIB "rootsur.lib"; // We use varsigns |
---|
[e182c8] | 41 | |
---|
[f5ba45] | 42 | proc nrRootsProbab(ideal I, list #) |
---|
| 43 | "USAGE: nrRootsProbab(I,[n]); ideal I, int n |
---|
| 44 | RETURN: int: the number of real roots of the ideal I by a probabilistic |
---|
| 45 | algorithm |
---|
[731e67e] | 46 | ASSUME: If I is not a Groebner basis, then a Groebner basis will be computed |
---|
[f5ba45] | 47 | by using std. If I is already a Groebner basis (i.e. if |
---|
| 48 | attrib(I,"isSB"); returns 1) then this Groebner basis will be |
---|
| 49 | used, hence it must be one w.r.t. (any) global ordering. This may |
---|
| 50 | be useful if the ideal is known to be a Groebner basis or if it |
---|
| 51 | can be computed faster by a different method. |
---|
| 52 | NOTE: If n<10 is given, n is the number of digits being used for |
---|
| 53 | constructing a random characteristic polynomial, a bigger n is |
---|
| 54 | more safe but slower (default: n=5). |
---|
| 55 | If printlevel>0 the number of complex solutions is displayed |
---|
| 56 | (default: printlevel=0). |
---|
| 57 | SEE ALSO: nrroots, nrRootsDeterm, randcharpoly, solve |
---|
| 58 | EXAMPLE: example nrRootsProbab; shows an example" |
---|
| 59 | { |
---|
| 60 | //Note on complexity: Let n = no of complex roots of I (= vdim(std(I)). |
---|
| 61 | //Then the algorithm needs: |
---|
| 62 | //1 std(I) and ~n NF computations (of randcharpoly wrt I) |
---|
| 63 | |
---|
| 64 | if (isparam(I)) { |
---|
| 65 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
| 66 | } |
---|
| 67 | int pr = printlevel-voice+2; |
---|
| 68 | int v; |
---|
| 69 | int n=5; |
---|
| 70 | if (size(#) == 1) { |
---|
| 71 | n=#[1]; |
---|
| 72 | } |
---|
| 73 | if (attrib(I,"isSB")!=1) { |
---|
| 74 | I = std(I); |
---|
| 75 | } |
---|
| 76 | |
---|
| 77 | ideal b = qbase(I); |
---|
| 78 | v = size(b); |
---|
| 79 | if (v == 0) { |
---|
| 80 | ERROR("ideal is not 0-dimensional"); |
---|
| 81 | } |
---|
| 82 | dbprint(pr,"//ideal has " +string(v)+ " complex solutions, counted with multiplicity"); |
---|
| 83 | |
---|
| 84 | poly p = randcharpoly(b,I,n); |
---|
| 85 | |
---|
| 86 | return (nrroots(p)); |
---|
| 87 | } |
---|
| 88 | |
---|
| 89 | example |
---|
| 90 | { |
---|
| 91 | echo = 2; |
---|
| 92 | ring r = 0,(x,y,z),lp; |
---|
| 93 | ideal i = (x-1)*(x-2),(y-1)^3*(x-y),(z-1)*(z-2)*(z-3)^2; |
---|
| 94 | nrRootsProbab(i); //no of real roots (using internally std) |
---|
| 95 | |
---|
| 96 | i = groebner(i); //using the hilbert driven GB computation |
---|
| 97 | int pr = printlevel; |
---|
| 98 | printlevel = 2; |
---|
| 99 | nrRootsProbab(i); |
---|
| 100 | printlevel = pr; |
---|
| 101 | } |
---|
| 102 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 103 | |
---|
| 104 | proc nrRootsDeterm(ideal I) |
---|
| 105 | "USAGE: nrRootsDeterm(I); ideal I |
---|
| 106 | RETURN: int: the number of real roots of the ideal I by a deterministic |
---|
| 107 | algorithm |
---|
[731e67e] | 108 | ASSUME: If I is not a Groebner basis, then a Groebner basis will be computed |
---|
[f5ba45] | 109 | by using std. If I is already a Groebner basis (i.e. if |
---|
| 110 | attrib(I,"isSB"); returns 1) then this Groebner basis will be |
---|
| 111 | used, hence it must be one w.r.t. (any) global ordering. This may |
---|
| 112 | be useful if the ideal is known to be a Groebner basis or if it |
---|
| 113 | can be computed faster by a different method. |
---|
| 114 | NOTE: If printlevel>0 the number of complex solutions is displayed |
---|
[731e67e] | 115 | (default: printlevel=0). The procedure nrRootsProbab is usually faster. |
---|
[f5ba45] | 116 | SEE ALSO: nrroots, nrRootsProbab, sturmquery, solve |
---|
| 117 | EXAMPLE: example nrRootsDeterm; shows an example" |
---|
| 118 | { |
---|
| 119 | //Note on complexity: Let n = no of complex roots of I (= vdim(std(I)). |
---|
| 120 | //Then the algotithm needs: |
---|
| 121 | //1 std(I) and (1/2)n*(n+1)^2 ~ 1/2n^3 NF computations (of monomials wrt I) |
---|
| 122 | |
---|
| 123 | if (isparam(I)) { |
---|
| 124 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
| 125 | } |
---|
| 126 | int pr = printlevel-voice+2; |
---|
| 127 | int v; |
---|
| 128 | |
---|
| 129 | if (attrib(I,"isSB")!=1) { |
---|
| 130 | I = std(I); |
---|
| 131 | } |
---|
| 132 | |
---|
| 133 | ideal b = qbase(I); |
---|
| 134 | v = size(b); |
---|
| 135 | if (v == 0) { |
---|
| 136 | ERROR("ideal is not 0-dimensional"); |
---|
| 137 | } |
---|
| 138 | dbprint(pr,"//ideal has " +string(v)+ " complex solutions, counted with multiplicity"); |
---|
| 139 | |
---|
| 140 | return (sturmquery(1,b,I)); |
---|
| 141 | } |
---|
| 142 | |
---|
| 143 | example |
---|
| 144 | { |
---|
| 145 | echo = 2; |
---|
| 146 | ring r = 0,(x,y,z),lp; |
---|
| 147 | ideal I = (x-1)*(x-2),(y-1),(z-1)*(z-2)*(z-3)^2; |
---|
| 148 | nrRootsDeterm(I); //no of real roots (using internally std) |
---|
| 149 | |
---|
| 150 | I = groebner(I); //using the hilbert driven GB computation |
---|
| 151 | int pr = printlevel; |
---|
| 152 | printlevel = 2; |
---|
| 153 | nrRootsDeterm(I); |
---|
| 154 | printlevel = pr; |
---|
| 155 | } |
---|
| 156 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 157 | |
---|
[e182c8] | 158 | proc symsignature(matrix m) |
---|
| 159 | "USAGE: symsignature(m); m matrix. m must be symmetric. |
---|
[f5ba45] | 160 | RETURN: int: the signature of m |
---|
[e182c8] | 161 | SEE ALSO: matbil,sturmquery |
---|
| 162 | EXAMPLE: example symsignature; shows an example" |
---|
| 163 | { |
---|
| 164 | int positive, negative, i, j; |
---|
| 165 | list l; |
---|
| 166 | poly variable; |
---|
| 167 | |
---|
| 168 | if (isparam(m)) { |
---|
| 169 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
| 170 | } |
---|
| 171 | |
---|
| 172 | if (!isSquare(m)) { |
---|
| 173 | ERROR ("m must be a square matrix"); |
---|
| 174 | } |
---|
| 175 | |
---|
| 176 | // We check whether m is symmetric |
---|
| 177 | for (i = 1;i <= nrows(m);i++) { |
---|
| 178 | for (j = i;j <= nrows(m);j++) { |
---|
| 179 | if (m[i,j] != m[j,i]) { |
---|
[d3b3ab] | 180 | ERROR ("m must be a symmetric matrix"); |
---|
[e182c8] | 181 | } |
---|
| 182 | } |
---|
| 183 | } |
---|
| 184 | |
---|
| 185 | poly f = charpoly(m); // Uses the last variable of the ring |
---|
| 186 | |
---|
| 187 | for (i = size(f);i >= 1;i--) { |
---|
| 188 | l[i] = leadcoef(f[i]); |
---|
| 189 | } |
---|
| 190 | positive = varsigns(l); |
---|
| 191 | |
---|
| 192 | variable = var(nvars(basering)); // charpoly uses the last variable |
---|
| 193 | f = subst(f,variable,-variable); |
---|
| 194 | |
---|
| 195 | for (i = size(f);i >= 1;i--) { |
---|
| 196 | l[i] = leadcoef(f[i]); |
---|
| 197 | } |
---|
| 198 | |
---|
| 199 | negative = varsigns(l); |
---|
| 200 | return (positive - negative); |
---|
| 201 | } |
---|
| 202 | example |
---|
| 203 | { |
---|
| 204 | echo = 2; |
---|
| 205 | ring r = 0,(x,y),dp; |
---|
| 206 | ideal i = x4-y2x,y2-13; |
---|
[d3b3ab] | 207 | i = std(i); |
---|
[e182c8] | 208 | ideal b = qbase(i); |
---|
| 209 | |
---|
| 210 | matrix m = matbil(1,b,i); |
---|
| 211 | symsignature(m); |
---|
| 212 | } |
---|
| 213 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 214 | |
---|
| 215 | proc sturmquery(poly h,ideal B,ideal I) |
---|
| 216 | "USAGE: sturmquery(h,b,i); h poly, b,i ideal |
---|
[f5ba45] | 217 | RETURN: int: the Sturm query of h in V(i) |
---|
[d3b3ab] | 218 | ASSUME: i is a Groebner basis, b is an ordered monomial basis |
---|
| 219 | of r/i, r = basering. |
---|
[e182c8] | 220 | SEE ALSO: symsignature,matbil |
---|
| 221 | EXAMPLE: example sturmquery; shows an example" |
---|
| 222 | { |
---|
| 223 | if (isparam(h) || isparam(B) || isparam(I)) { |
---|
| 224 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
| 225 | } |
---|
| 226 | |
---|
| 227 | return (mysymmsig(matbil(h,B,I))); |
---|
| 228 | } |
---|
| 229 | example |
---|
| 230 | { |
---|
| 231 | echo = 2; |
---|
| 232 | ring r = 0,(x,y),dp; |
---|
| 233 | ideal i = x4-y2x,y2-13; |
---|
[d3b3ab] | 234 | i = std(i); |
---|
[e182c8] | 235 | ideal b = qbase(i); |
---|
| 236 | |
---|
| 237 | sturmquery(1,b,i); |
---|
| 238 | } |
---|
[f5ba45] | 239 | /////////////////////////////////////////////////////////////////////////////// |
---|
[e182c8] | 240 | |
---|
[d3b3ab] | 241 | static proc mysymmsig(matrix m) |
---|
[e182c8] | 242 | // returns the signature of a square symmetric matrix m |
---|
| 243 | { |
---|
| 244 | int positive, negative, i; |
---|
| 245 | list l; |
---|
| 246 | poly variable; |
---|
| 247 | |
---|
| 248 | poly f = charpoly(m); // Uses the last variable of the ring |
---|
| 249 | |
---|
| 250 | for (i = size(f);i >= 1;i--) { |
---|
| 251 | l[i] = leadcoef(f[i]); |
---|
| 252 | } |
---|
| 253 | positive = varsigns(l); |
---|
| 254 | |
---|
| 255 | variable = var(nvars(basering)); // charpoly uses the last variable |
---|
| 256 | f = subst(f,variable,-variable); |
---|
| 257 | |
---|
| 258 | for (i = size(f);i >= 1;i--) { |
---|
| 259 | l[i] = leadcoef(f[i]); |
---|
| 260 | } |
---|
| 261 | |
---|
| 262 | negative = varsigns(l); |
---|
| 263 | return (positive - negative); |
---|
| 264 | } |
---|
| 265 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 266 | |
---|
| 267 | proc matbil(poly h,ideal B,ideal I) |
---|
| 268 | "USAGE: matbil(h,b,i); h poly, b,i ideal |
---|
| 269 | RETURN: matrix: the matrix of the bilinear form (f,g) |-> trace(m_fhg), |
---|
| 270 | m_fhg = multiplication with fhg on r/i |
---|
| 271 | ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, |
---|
| 272 | r = basering |
---|
| 273 | SEE ALSO: matmult,tracemult |
---|
| 274 | EXAMPLE: example matbil; shows an example" |
---|
| 275 | { |
---|
| 276 | matrix m[size(B)][size(B)]; |
---|
| 277 | poly f; |
---|
| 278 | int k,l; |
---|
[f5ba45] | 279 | //h = reduce(h,I); |
---|
[e182c8] | 280 | |
---|
[f5ba45] | 281 | for (k = 1; k <= size(B); k++) { |
---|
| 282 | for (l = 1; l <= k; l++) { |
---|
| 283 | m[k,l] = tracemult(h*B[k]*B[l],B,I)[1]; |
---|
[e182c8] | 284 | m[l,k] = m[k,l]; // The matrix we are trying to compute is symmetric |
---|
| 285 | } |
---|
[f5ba45] | 286 | } |
---|
[e182c8] | 287 | return(m); |
---|
| 288 | } |
---|
| 289 | example |
---|
| 290 | { |
---|
| 291 | echo = 2; |
---|
| 292 | ring r = 0,(x,y),dp; |
---|
| 293 | ideal i = x4-y2x,y2-13; |
---|
[d3b3ab] | 294 | i = std(i); |
---|
[e182c8] | 295 | ideal b = qbase(i); |
---|
| 296 | poly f = x3-xy+y-13+x4-y2x; |
---|
| 297 | |
---|
| 298 | matrix m = matbil(f,b,i); |
---|
| 299 | print(m); |
---|
| 300 | |
---|
| 301 | } |
---|
| 302 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 303 | |
---|
| 304 | proc tracemult(poly f,ideal B,ideal I) |
---|
| 305 | "USAGE: tracemult(f,b,i);f poly, b,i ideal |
---|
[d3b3ab] | 306 | RETURN: number: the trace of the multiplication by f (m_f) on r/i, |
---|
| 307 | written in the monomial basis b of r/i, r = basering |
---|
[e182c8] | 308 | (faster than matmult + trace) |
---|
| 309 | ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i |
---|
| 310 | SEE ALSO: matmult,trace |
---|
| 311 | EXAMPLE: example tracemult; shows an example" |
---|
| 312 | { |
---|
| 313 | int k; // Iterates over the basis monomials |
---|
| 314 | int l; // Iterates over the rows of the matrix |
---|
| 315 | list coordinates; |
---|
| 316 | number m; |
---|
[f5ba45] | 317 | poly g; |
---|
[e182c8] | 318 | |
---|
[f5ba45] | 319 | //f = reduce(f,I); |
---|
| 320 | for (k = 1; k <= size(B); k++) { |
---|
| 321 | l=1; |
---|
| 322 | g = reduce(f*B[k],I); |
---|
| 323 | while (l <= k) { |
---|
| 324 | if (leadmonom(g[l]) == B[k]) { |
---|
| 325 | m = m + leadcoef(g[l]); |
---|
| 326 | break; |
---|
| 327 | } |
---|
| 328 | l++; |
---|
| 329 | } |
---|
[e182c8] | 330 | } |
---|
| 331 | return (m); |
---|
| 332 | } |
---|
| 333 | example |
---|
| 334 | { |
---|
| 335 | echo = 2; |
---|
| 336 | ring r = 0,(x,y),dp; |
---|
| 337 | ideal i = x4-y2x,y2-13; |
---|
[d3b3ab] | 338 | i = std(i); |
---|
[e182c8] | 339 | ideal b = qbase(i); |
---|
| 340 | |
---|
| 341 | poly f = x3-xy+y-13+x4-y2x; |
---|
| 342 | matrix m = matmult(f,b,i); |
---|
| 343 | print(m); |
---|
| 344 | |
---|
[f5ba45] | 345 | tracemult(f,b,i); //the trace of m |
---|
[e182c8] | 346 | } |
---|
| 347 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 348 | |
---|
| 349 | proc matmult(poly f, ideal B, ideal I) |
---|
| 350 | "USAGE: matmult(f,b,i); f poly, b,i ideal |
---|
| 351 | RETURN: matrix: the matrix of the multiplication map by f (m_f) on r/i |
---|
[f5ba45] | 352 | w.r.t. to the monomial basis b of r/i (r = basering) |
---|
| 353 | ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, |
---|
| 354 | as given by qbase(i) |
---|
[e182c8] | 355 | SEE ALSO: coords,matbil |
---|
| 356 | EXAMPLE: example matmult; shows an example" |
---|
| 357 | { |
---|
| 358 | int k; // Iterates over the basis monomials |
---|
| 359 | int l; // Iterates over the rows of the matrix |
---|
| 360 | list coordinates; |
---|
| 361 | matrix m[size(B)][size(B)]; |
---|
| 362 | |
---|
[f5ba45] | 363 | //f = reduce(f,I); |
---|
[e182c8] | 364 | for (k = 1;k <= size(B);k++) { |
---|
[f5ba45] | 365 | coordinates = coords(f*(B[k]),B,I); // f*x_k written on the basis B |
---|
[e182c8] | 366 | for (l = 1;l <= size(B);l++) { |
---|
| 367 | m[l,k] = coordinates[l]; |
---|
| 368 | } |
---|
| 369 | } |
---|
| 370 | return (m); |
---|
| 371 | } |
---|
| 372 | example |
---|
| 373 | { |
---|
| 374 | echo = 2; |
---|
| 375 | ring r = 0,(x,y),dp; |
---|
| 376 | ideal i = x4-y2x,y2-13; |
---|
[d3b3ab] | 377 | i = std(i); |
---|
[e182c8] | 378 | ideal b = qbase(i); |
---|
| 379 | |
---|
| 380 | poly f = x3-xy+y-13+x4-y2x; |
---|
| 381 | matrix m = matmult(f,b,i); |
---|
| 382 | print(m); |
---|
| 383 | } |
---|
| 384 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 385 | |
---|
| 386 | proc coords(poly f,ideal B,ideal I) |
---|
| 387 | "USAGE: coords(f,b,i), f poly, b,i ideal |
---|
[f5ba45] | 388 | RETURN: list of numbers: the coordinates of the class of f (mod i) |
---|
| 389 | in the monomial basis b |
---|
[e182c8] | 390 | ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, |
---|
| 391 | r = basering |
---|
| 392 | SEE ALSO: matmult,matbil |
---|
| 393 | KEYWORDS: coordinates |
---|
| 394 | EXAMPLE: example coords; shows an example" |
---|
| 395 | { |
---|
| 396 | // We assume the basis is sorted according to the ring order |
---|
| 397 | poly g; |
---|
| 398 | int k,l=1,1; |
---|
| 399 | list coordinates; |
---|
| 400 | int N = size(B); |
---|
| 401 | |
---|
| 402 | // We first compute the normal form of f wrt I |
---|
| 403 | g = reduce(f,I); |
---|
[f5ba45] | 404 | int n = size(g); //allways n <= N |
---|
[e182c8] | 405 | |
---|
| 406 | while (k <= N) { |
---|
[f5ba45] | 407 | if (leadmonom(g[l]) == B[k]) { |
---|
[e182c8] | 408 | coordinates[k] = leadcoef(g[l]); |
---|
| 409 | l++; |
---|
| 410 | } else { |
---|
[f5ba45] | 411 | coordinates[k] = number(0); |
---|
[e182c8] | 412 | } |
---|
| 413 | k++; |
---|
| 414 | } |
---|
| 415 | return (coordinates); |
---|
| 416 | } |
---|
| 417 | example |
---|
| 418 | { |
---|
| 419 | echo = 2; |
---|
| 420 | ring r = 0,(x,y),dp; |
---|
| 421 | ideal i = x4-y2x,y2-13; |
---|
[f5ba45] | 422 | poly f = x3-xy+y-13+x4-y2x; |
---|
[d3b3ab] | 423 | i = std(i); |
---|
[e182c8] | 424 | ideal b = qbase(i); |
---|
| 425 | b; |
---|
[f5ba45] | 426 | coords(f,b,i); |
---|
[e182c8] | 427 | } |
---|
| 428 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 429 | |
---|
[d3b3ab] | 430 | static proc isSquare(matrix m) |
---|
[e182c8] | 431 | // returns 1 iff m is a square matrix |
---|
| 432 | { |
---|
| 433 | return (nrows(m)==ncols(m)); |
---|
| 434 | } |
---|
| 435 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 436 | |
---|
[d3b3ab] | 437 | proc randcharpoly(ideal B,ideal I,list #) |
---|
| 438 | "USAGE: randcharpoly(b,i); randcharpoly(b,i,n); b,i ideal; n int |
---|
[e182c8] | 439 | RETURN: poly: the characteristic polynomial of a pseudorandom |
---|
[d3b3ab] | 440 | rational univariate projection having one zero per zero of i. |
---|
[f5ba45] | 441 | If n<10 is given, it is the number of digits being used for the |
---|
| 442 | pseudorandom coefficients (default: n=5) |
---|
[e182c8] | 443 | ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, |
---|
| 444 | r = basering |
---|
[f5ba45] | 445 | NOTE: shows a warning if printlevel>0 (default: printlevel=0) |
---|
[e182c8] | 446 | KEYWORDS: rational univariate projection |
---|
| 447 | EXAMPLE: example randcharpoly; shows an example" |
---|
| 448 | { |
---|
[f5ba45] | 449 | int pr = printlevel - voice + 2; |
---|
[e182c8] | 450 | poly p; |
---|
| 451 | poly generic; |
---|
| 452 | list l; |
---|
| 453 | matrix m; |
---|
| 454 | poly q; |
---|
| 455 | |
---|
[d3b3ab] | 456 | if (size(#) == 1) { |
---|
| 457 | generic = randlinpoly(#[1]); |
---|
| 458 | } else { |
---|
| 459 | generic = randlinpoly(); |
---|
| 460 | } |
---|
[e182c8] | 461 | |
---|
| 462 | p = reduce(generic,I); |
---|
| 463 | m = matmult(p,B,I); |
---|
| 464 | q = charpoly(m); |
---|
| 465 | |
---|
[f5ba45] | 466 | dbprint(pr,"*********************************************************************"); |
---|
| 467 | dbprint(pr,"* WARNING: This polynomial was obtained using pseudorandom numbers.*"); |
---|
| 468 | dbprint(pr,"* If you want to verify the result, please use the command *"); |
---|
| 469 | dbprint(pr,"* *"); |
---|
| 470 | dbprint(pr,"* verify(p,b,i) *"); |
---|
| 471 | dbprint(pr,"* *"); |
---|
| 472 | dbprint(pr,"* where p is the polynomial I returned, b is the monomial basis *"); |
---|
| 473 | dbprint(pr,"* used, and i the Groebner basis of the ideal *"); |
---|
| 474 | dbprint(pr,"*********************************************************************"); |
---|
[e182c8] | 475 | |
---|
| 476 | return(q); |
---|
| 477 | } |
---|
| 478 | example |
---|
| 479 | { |
---|
| 480 | echo = 2; |
---|
| 481 | ring r = 0,(x,y,z),dp; |
---|
| 482 | ideal i = (x-1)*(x-2),(y-1),(z-1)*(z-2)*(z-3)^2; |
---|
[d3b3ab] | 483 | i = std(i); |
---|
[e182c8] | 484 | ideal b = qbase(i); |
---|
| 485 | poly p = randcharpoly(b,i); |
---|
| 486 | p; |
---|
| 487 | nrroots(p); // See nrroots in urrcount.lib |
---|
[f5ba45] | 488 | |
---|
| 489 | int pr = printlevel; |
---|
| 490 | printlevel = pr+2; |
---|
[d3b3ab] | 491 | p = randcharpoly(b,i,5); |
---|
| 492 | nrroots(p); |
---|
[f5ba45] | 493 | printlevel = pr; |
---|
[e182c8] | 494 | } |
---|
[d3b3ab] | 495 | |
---|
[e182c8] | 496 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 497 | |
---|
| 498 | proc verify(poly p,ideal b,ideal i) |
---|
| 499 | "USAGE: verify(p,b,i);p poly, b,i,ideal |
---|
[d3b3ab] | 500 | RETURN: integer: 1 iff the polynomial p splits the points of V(i). |
---|
[e182c8] | 501 | It's used to check the result of randcharpoly |
---|
| 502 | ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, |
---|
| 503 | r = basering |
---|
[f5ba45] | 504 | NOTE: comments the result if printlevel>0 (default: printlevel=0) |
---|
[e182c8] | 505 | SEE ALSO: randcharpoly |
---|
| 506 | EXAMPLE: example verify; shows an example" |
---|
| 507 | { |
---|
[f5ba45] | 508 | int pr = printlevel - voice + 2; |
---|
[741c669] | 509 | poly sqr_free; |
---|
[e182c8] | 510 | int correct; |
---|
| 511 | poly variable; |
---|
| 512 | |
---|
| 513 | if (isparam(p) || isparam(b) || isparam(i)) { |
---|
| 514 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
| 515 | } |
---|
| 516 | |
---|
| 517 | variable = isuni(p); |
---|
[741c669] | 518 | sqr_free = p/gcd(p,diff(p,variable)); |
---|
| 519 | correct = (mat_rk(matbil(1,b,i)) == deg(sqr_free)); |
---|
[e182c8] | 520 | |
---|
| 521 | if (correct) { |
---|
[f5ba45] | 522 | dbprint(pr,"//Verification successful"); |
---|
[e182c8] | 523 | } else { |
---|
[f5ba45] | 524 | dbprint(pr,"//The choice of random numbers was not useful"); |
---|
| 525 | dbprint(pr,"//You might want to try randcharpoly with a larger number of digits"); |
---|
[e182c8] | 526 | } |
---|
| 527 | return (correct); |
---|
| 528 | } |
---|
| 529 | example |
---|
| 530 | { |
---|
| 531 | echo = 2; |
---|
| 532 | ring r = 0,(x,y),dp; |
---|
| 533 | poly f = x3-xy+y-13+x4-y2x; |
---|
| 534 | ideal i = x4-y2x,y2-13; |
---|
[d3b3ab] | 535 | i = std(i); |
---|
[e182c8] | 536 | ideal b = qbase(i); |
---|
| 537 | poly p = randcharpoly(b,i); |
---|
| 538 | verify(p,b,i); |
---|
| 539 | } |
---|
| 540 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 541 | |
---|
[d3b3ab] | 542 | proc randlinpoly(list #) |
---|
| 543 | "USAGE: randlinpoly(); randlinpoly(n); n int |
---|
[f5ba45] | 544 | RETURN: poly: linear combination of the variables of the ring, with |
---|
| 545 | pseudorandom coefficients. If n<10 is given, it is the number of |
---|
| 546 | digits being used for the range of the coefficients (default: n=5) |
---|
[e182c8] | 547 | SEE ALSO: randcharpoly; |
---|
| 548 | EXAMPLE: example randlinpoly; shows an example" |
---|
| 549 | { |
---|
| 550 | int n,i; |
---|
| 551 | poly p = 0; |
---|
[f5ba45] | 552 | int ndigits = 5; |
---|
[d3b3ab] | 553 | |
---|
| 554 | if (size(#) == 1) { |
---|
| 555 | ndigits = #[1]; |
---|
| 556 | } |
---|
| 557 | |
---|
[e182c8] | 558 | n = nvars(basering); |
---|
| 559 | for (i = 1;i <= n;i++) { |
---|
[d3b3ab] | 560 | p = p + var(i)*random(1,10^ndigits); |
---|
[e182c8] | 561 | } |
---|
| 562 | return (p); |
---|
| 563 | } |
---|
| 564 | example |
---|
| 565 | { |
---|
| 566 | echo = 2; |
---|
| 567 | ring r = 0,(x,y,z,w),dp; |
---|
| 568 | poly p = randlinpoly(); |
---|
| 569 | p; |
---|
[f5ba45] | 570 | randlinpoly(5); |
---|
[e182c8] | 571 | } |
---|
| 572 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 573 | |
---|
| 574 | proc powersums(poly f,ideal B,ideal I) |
---|
[f5ba45] | 575 | "USAGE: powersums(f,b,i); f poly; b,i ideal |
---|
[e182c8] | 576 | RETURN: list: the powersums of the results of evaluating f at the zeros of I |
---|
[f5ba45] | 577 | ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, |
---|
| 578 | r = basering |
---|
[e182c8] | 579 | SEE ALSO: symmfunc |
---|
| 580 | EXAMPLE: example symmfunc; shows an example" |
---|
| 581 | { |
---|
[d3b3ab] | 582 | int N,k; |
---|
[e182c8] | 583 | list sums; |
---|
| 584 | |
---|
| 585 | N = size(B); |
---|
| 586 | for (k = 1;k <= N;k++) { |
---|
| 587 | sums = sums + list(leadcoef(trace(matmult(f^k,B,I)))); |
---|
| 588 | } |
---|
| 589 | return (sums); |
---|
| 590 | } |
---|
| 591 | example |
---|
| 592 | { |
---|
| 593 | echo = 2; |
---|
| 594 | ring r = 0,(x,y,z),dp; |
---|
| 595 | |
---|
| 596 | ideal i = (x-1)*(x-2),(y-1),(z+5); // V(I) = {(1,1,-5),(2,1,-5) |
---|
[d3b3ab] | 597 | i = std(i); |
---|
[e182c8] | 598 | |
---|
| 599 | ideal b = qbase(i); |
---|
| 600 | poly f = x+y+z; |
---|
| 601 | list psums = list(-2-3,4+9); // f evaluated at V(I) gives {-3,-2} |
---|
| 602 | list l = powersums(f,b,i); |
---|
| 603 | psums; |
---|
| 604 | l; |
---|
| 605 | } |
---|
| 606 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 607 | |
---|
[d3b3ab] | 608 | proc symmfunc(list S) |
---|
[e182c8] | 609 | "USAGE: symmfunc(s); s list |
---|
| 610 | RETURN: list: the symmetric functions of the roots of a polynomial, given |
---|
| 611 | the power sums of those roots. |
---|
| 612 | SEE ALSO: powersums |
---|
| 613 | EXAMPLE: example symmfunc; shows an example" |
---|
| 614 | { |
---|
[f5ba45] | 615 | // Takes the list of power sums and returns the symmetric functions |
---|
[e182c8] | 616 | list a; |
---|
| 617 | int j,l,N; |
---|
| 618 | number sum; |
---|
| 619 | |
---|
| 620 | N = size(S); |
---|
| 621 | a[N+1] = 1; // We set the length of the list and initialize its last element. |
---|
| 622 | |
---|
| 623 | for (l = N - 1;l >= 0;l--) { |
---|
| 624 | sum = 0; |
---|
| 625 | for (j = l + 1;j <= N;j++) { |
---|
| 626 | sum = sum + ((a[j+1])*(S[j-l])); |
---|
| 627 | } |
---|
| 628 | sum = -sum; |
---|
| 629 | a[l+1] = sum/(N-l); |
---|
| 630 | } |
---|
| 631 | |
---|
| 632 | a = reverse(a); |
---|
| 633 | return (a); |
---|
| 634 | } |
---|
| 635 | example |
---|
| 636 | { |
---|
| 637 | echo = 2; |
---|
| 638 | ring r = 0,x,dp; |
---|
| 639 | poly p = (x-1)*(x-2)*(x-3); |
---|
| 640 | list psums = list(1+2+3,1+4+9,1+8+27); |
---|
| 641 | list l = symmfunc(psums); |
---|
| 642 | l; |
---|
| 643 | p; // Compare p with the elements of l |
---|
| 644 | } |
---|
| 645 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 646 | |
---|
| 647 | proc univarpoly(list l) |
---|
| 648 | "USAGE: univarpoly(l); l list |
---|
[d3b3ab] | 649 | RETURN: poly: a polynomial p on the first variable of basering, say x, |
---|
[e182c8] | 650 | with p = l[1] + l[2]*x + l[3]*x^2 + ... |
---|
| 651 | EXAMPLE: example univarpoly; shows an example" |
---|
| 652 | { |
---|
| 653 | poly p; |
---|
| 654 | int i,n; |
---|
[d3b3ab] | 655 | |
---|
[e182c8] | 656 | n = size(l); |
---|
| 657 | for (i = 1;i <= n;i++) { |
---|
| 658 | p = p + l[i]*var(1)^(n-i); |
---|
| 659 | } |
---|
| 660 | return (p); |
---|
| 661 | } |
---|
| 662 | example |
---|
| 663 | { |
---|
| 664 | echo = 2; |
---|
| 665 | ring r = 0,x,dp; |
---|
| 666 | list l = list(1,2,3,4,5); |
---|
| 667 | poly p = univarpoly(l); |
---|
| 668 | p; |
---|
| 669 | } |
---|
| 670 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 671 | |
---|
| 672 | proc qbase(ideal i) |
---|
| 673 | "USAGE: qbase(I); I zero-dimensional ideal |
---|
| 674 | RETURN: ideal: A monomial basis of the quotient between the basering and the |
---|
| 675 | ideal I, sorted according to the basering order. |
---|
| 676 | SEE ALSO: kbase |
---|
| 677 | KEYWORDS: zero-dimensional |
---|
| 678 | EXAMPLE: example qbase; shows an example" |
---|
| 679 | { |
---|
| 680 | ideal b; |
---|
| 681 | |
---|
| 682 | b = kbase(i); |
---|
| 683 | b = reverseideal(sort(b)[1]); // sort sorts in ascending order |
---|
| 684 | return (b); |
---|
| 685 | } |
---|
| 686 | example |
---|
| 687 | { |
---|
| 688 | echo = 2; |
---|
| 689 | ring r = 0,(x,y,z),dp; |
---|
| 690 | |
---|
| 691 | ideal i = 2x2,-y2,z3; |
---|
[d3b3ab] | 692 | i = std(i); |
---|
[e182c8] | 693 | ideal b = qbase(i); |
---|
| 694 | b; |
---|
| 695 | b = kbase(i); |
---|
| 696 | b; // Compare this with the result of qbase |
---|
| 697 | } |
---|
| 698 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 699 | |
---|
| 700 | static proc reverseideal(ideal b) // Returns b reversed |
---|
| 701 | { |
---|
| 702 | int i; |
---|
| 703 | ideal result; |
---|
| 704 | |
---|
| 705 | result = b[1]; |
---|
| 706 | for (i = 2;i <= size(b);i++) { |
---|
| 707 | result = b[i], result; |
---|
| 708 | } |
---|
| 709 | return (result); |
---|
| 710 | } |
---|
| 711 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 712 | |
---|