/////////////////////////////////////////////////////////////////////////////// version="$Id$"; category="Teaching"; info=" LIBRARY: rootsmr.lib Counting the number of real roots of polynomial systems AUTHOR: Enrique A. Tobis, etobis@dc.uba.ar OVERVIEW: Routines for counting the number of real roots of a multivariate polynomial system. Two methods are implemented: deterministic computation of the number of roots, via the signature of a certain bilinear form (nrRootsDeterm); and a rational univariate projection, using a pseudorandom polynomial (nrRootsProbab). It also includes a command to verify the correctness of the pseudorandom answer. References: Basu, Pollack, Roy, \"Algorithms in Real Algebraic Geometry\", Springer, 2003. PROCEDURES: nrRootsProbab(I) Number of real roots of 0-dim ideal (probabilistic) nrRootsDeterm(I) Number of real roots of 0-dim ideal (deterministic) symsignature(m) Signature of the symmetric matrix m sturmquery(h,B,I) Sturm query of h on V(I) matbil(h,B,I) Matrix of the bilinear form on R/I associated to h matmult(f,B,I) Matrix of multiplication by f (m_f) on R/I in the basis B tracemult(f,B,I) Trace of m_f (B is an ordered basis of R/I) coords(f,B,I) Coordinates of f in the ordered basis B randcharpoly(B,I,n) Pseudorandom charpoly of univ. projection, n optional verify(p,B,i) Verifies the result of randcharpoly randlinpoly(n) Pseudorandom linear polynomial, n optional powersums(f,B,I) Powersums of the roots of a char polynomial symmfunc(S) Symmetric functions from the powersums S univarpoly(l) Polynomial with coefficients from l qbase(i) Like kbase, but the monomials are ordered KEYWORDS: real roots, univariate projection "; /////////////////////////////////////////////////////////////////// LIB "linalg.lib"; // We use charpoly LIB "rootsur.lib"; // We use varsigns proc nrRootsProbab(ideal I, list #) "USAGE: nrRootsProbab(I,[n]); ideal I, int n RETURN: int: the number of real roots of the ideal I by a probabilistic algorithm ASSUME: If I is not a Groebner basis, then a Groebner basis will be computed by using std. If I is already a Groebner basis (i.e. if attrib(I,"isSB"); returns 1) then this Groebner basis will be used, hence it must be one w.r.t. (any) global ordering. This may be useful if the ideal is known to be a Groebner basis or if it can be computed faster by a different method. NOTE: If n<10 is given, n is the number of digits being used for constructing a random characteristic polynomial, a bigger n is more safe but slower (default: n=5). If printlevel>0 the number of complex solutions is displayed (default: printlevel=0). SEE ALSO: nrroots, nrRootsDeterm, randcharpoly, solve EXAMPLE: example nrRootsProbab; shows an example" { //Note on complexity: Let n = no of complex roots of I (= vdim(std(I)). //Then the algorithm needs: //1 std(I) and ~n NF computations (of randcharpoly w.r.t. I) if (isparam(I)) { ERROR("This procedure cannot operate with parametric arguments"); } int pr = printlevel-voice+2; int v; int n=5; if (size(#) == 1) { n=#[1]; } if (attrib(I,"isSB")!=1) { I = std(I); } ideal b = qbase(I); v = size(b); if (v == 0) { ERROR("ideal is not 0-dimensional"); } dbprint(pr,"//ideal has " +string(v)+ " complex solutions, counted with multiplicity"); poly p = randcharpoly(b,I,n); return (nrroots(p)); } example { echo = 2; ring r = 0,(x,y,z),lp; ideal i = (x-1)*(x-2),(y-1)^3*(x-y),(z-1)*(z-2)*(z-3)^2; nrRootsProbab(i); //no of real roots (using internally std) i = groebner(i); //using the hilbert driven GB computation int pr = printlevel; printlevel = 2; nrRootsProbab(i); printlevel = pr; } /////////////////////////////////////////////////////////////////////////////// proc nrRootsDeterm(ideal I) "USAGE: nrRootsDeterm(I); ideal I RETURN: int: the number of real roots of the ideal I by a deterministic algorithm ASSUME: If I is not a Groebner basis, then a Groebner basis will be computed by using std. If I is already a Groebner basis (i.e. if attrib(I,"isSB"); returns 1) then this Groebner basis will be used, hence it must be one w.r.t. (any) global ordering. This may be useful if the ideal is known to be a Groebner basis or if it can be computed faster by a different method. NOTE: If printlevel>0 the number of complex solutions is displayed (default: printlevel=0). The procedure nrRootsProbab is usually faster. SEE ALSO: nrroots, nrRootsProbab, sturmquery, solve EXAMPLE: example nrRootsDeterm; shows an example" { //Note on complexity: Let n = no of complex roots of I (= vdim(std(I)). //Then the algotithm needs: //1 std(I) and (1/2)n*(n+1)^2 ~ 1/2n^3 NF computations (of monomials w.r.t. I) if (isparam(I)) { ERROR("This procedure cannot operate with parametric arguments"); } int pr = printlevel-voice+2; int v; if (attrib(I,"isSB")!=1) { I = std(I); } ideal b = qbase(I); v = size(b); if (v == 0) { ERROR("ideal is not 0-dimensional"); } dbprint(pr,"//ideal has " +string(v)+ " complex solutions, counted with multiplicity"); return (sturmquery(1,b,I)); } example { echo = 2; ring r = 0,(x,y,z),lp; ideal I = (x-1)*(x-2),(y-1),(z-1)*(z-2)*(z-3)^2; nrRootsDeterm(I); //no of real roots (using internally std) I = groebner(I); //using the hilbert driven GB computation int pr = printlevel; printlevel = 2; nrRootsDeterm(I); printlevel = pr; } /////////////////////////////////////////////////////////////////////////////// proc symsignature(matrix m) "USAGE: symsignature(m); m matrix. m must be symmetric. RETURN: int: the signature of m SEE ALSO: matbil,sturmquery EXAMPLE: example symsignature; shows an example" { int positive, negative, i, j; list l; poly variable; if (isparam(m)) { ERROR("This procedure cannot operate with parametric arguments"); } if (!isSquare(m)) { ERROR ("m must be a square matrix"); } // We check whether m is symmetric for (i = 1;i <= nrows(m);i++) { for (j = i;j <= nrows(m);j++) { if (m[i,j] != m[j,i]) { ERROR ("m must be a symmetric matrix"); } } } poly f = charpoly(m); // Uses the last variable of the ring for (i = size(f);i >= 1;i--) { l[i] = leadcoef(f[i]); } positive = varsigns(l); variable = var(nvars(basering)); // charpoly uses the last variable f = subst(f,variable,-variable); for (i = size(f);i >= 1;i--) { l[i] = leadcoef(f[i]); } negative = varsigns(l); return (positive - negative); } example { echo = 2; ring r = 0,(x,y),dp; ideal i = x4-y2x,y2-13; i = std(i); ideal b = qbase(i); matrix m = matbil(1,b,i); symsignature(m); } /////////////////////////////////////////////////////////////////////////////// proc sturmquery(poly h,ideal B,ideal I) "USAGE: sturmquery(h,b,i); h poly, b,i ideal RETURN: int: the Sturm query of h in V(i) ASSUME: i is a Groebner basis, b is an ordered monomial basis of r/i, r = basering. SEE ALSO: symsignature,matbil EXAMPLE: example sturmquery; shows an example" { if (isparam(h) || isparam(B) || isparam(I)) { ERROR("This procedure cannot operate with parametric arguments"); } return (mysymmsig(matbil(h,B,I))); } example { echo = 2; ring r = 0,(x,y),dp; ideal i = x4-y2x,y2-13; i = std(i); ideal b = qbase(i); sturmquery(1,b,i); } /////////////////////////////////////////////////////////////////////////////// static proc mysymmsig(matrix m) // returns the signature of a square symmetric matrix m { int positive, negative, i; list l; poly variable; poly f = charpoly(m); // Uses the last variable of the ring for (i = size(f);i >= 1;i--) { l[i] = leadcoef(f[i]); } positive = varsigns(l); variable = var(nvars(basering)); // charpoly uses the last variable f = subst(f,variable,-variable); for (i = size(f);i >= 1;i--) { l[i] = leadcoef(f[i]); } negative = varsigns(l); return (positive - negative); } /////////////////////////////////////////////////////////////////////////////// proc matbil(poly h,ideal B,ideal I) "USAGE: matbil(h,b,i); h poly, b,i ideal RETURN: matrix: the matrix of the bilinear form (f,g) |-> trace(m_fhg), m_fhg = multiplication with fhg on r/i ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, r = basering SEE ALSO: matmult,tracemult EXAMPLE: example matbil; shows an example" { matrix m[size(B)][size(B)]; poly f; int k,l; //h = reduce(h,I); for (k = 1; k <= size(B); k++) { for (l = 1; l <= k; l++) { m[k,l] = tracemult(h*B[k]*B[l],B,I)[1]; m[l,k] = m[k,l]; // The matrix we are trying to compute is symmetric } } return(m); } example { echo = 2; ring r = 0,(x,y),dp; ideal i = x4-y2x,y2-13; i = std(i); ideal b = qbase(i); poly f = x3-xy+y-13+x4-y2x; matrix m = matbil(f,b,i); print(m); } /////////////////////////////////////////////////////////////////////////////// proc tracemult(poly f,ideal B,ideal I) "USAGE: tracemult(f,B,I);f poly, B,I ideal RETURN: number: the trace of the multiplication by f (m_f) on r/I, written in the monomial basis B of r/I, r = basering (faster than matmult + trace) ASSUME: I is given by a Groebner basis and B is an ordered monomial basis of r/I SEE ALSO: matmult,trace EXAMPLE: example tracemult; shows an example" { int k; // Iterates over the basis monomials int l; // Iterates over the rows of the matrix list coordinates; number m; poly g; //f = reduce(f,I); for (k = 1; k <= size(B); k++) { l=1; g = reduce(f*B[k],I); while (l <= k) { if (leadmonom(g[l]) == B[k]) { m = m + leadcoef(g[l]); break; } l++; } } return (m); } example { echo = 2; ring r = 0,(x,y),dp; ideal i = x4-y2x,y2-13; i = std(i); ideal b = qbase(i); poly f = x3-xy+y-13+x4-y2x; matrix m = matmult(f,b,i); print(m); tracemult(f,b,i); //the trace of m } /////////////////////////////////////////////////////////////////////////////// proc matmult(poly f, ideal B, ideal I) "USAGE: matmult(f,b,i); f poly, b,i ideal RETURN: matrix: the matrix of the multiplication map by f (m_f) on r/i w.r.t. to the monomial basis b of r/i (r = basering) ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, as given by qbase(i) SEE ALSO: coords,matbil EXAMPLE: example matmult; shows an example" { int k; // Iterates over the basis monomials int l; // Iterates over the rows of the matrix list coordinates; matrix m[size(B)][size(B)]; //f = reduce(f,I); for (k = 1;k <= size(B);k++) { coordinates = coords(f*(B[k]),B,I); // f*x_k written on the basis B for (l = 1;l <= size(B);l++) { m[l,k] = coordinates[l]; } } return (m); } example { echo = 2; ring r = 0,(x,y),dp; ideal i = x4-y2x,y2-13; i = std(i); ideal b = qbase(i); poly f = x3-xy+y-13+x4-y2x; matrix m = matmult(f,b,i); print(m); } /////////////////////////////////////////////////////////////////////////////// proc coords(poly f,ideal B,ideal I) "USAGE: coords(f,b,i), f poly, b,i ideal RETURN: list of numbers: the coordinates of the class of f (mod i) in the monomial basis b ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, r = basering SEE ALSO: matmult,matbil KEYWORDS: coordinates EXAMPLE: example coords; shows an example" { // We assume the basis is sorted according to the ring order poly g; int k,l=1,1; list coordinates; int N = size(B); // We first compute the normal form of f w.r.t. I g = reduce(f,I); int n = size(g); //allways n <= N while (k <= N) { if (leadmonom(g[l]) == B[k]) { coordinates[k] = leadcoef(g[l]); l++; } else { coordinates[k] = number(0); } k++; } return (coordinates); } example { echo = 2; ring r = 0,(x,y),dp; ideal i = x4-y2x,y2-13; poly f = x3-xy+y-13+x4-y2x; i = std(i); ideal b = qbase(i); b; coords(f,b,i); } /////////////////////////////////////////////////////////////////////////////// static proc isSquare(matrix m) // returns 1 if and only if m is a square matrix { return (nrows(m)==ncols(m)); } /////////////////////////////////////////////////////////////////////////////// proc randcharpoly(ideal B,ideal I,list #) "USAGE: randcharpoly(b,i); randcharpoly(b,i,n); b,i ideal; n int RETURN: poly: the characteristic polynomial of a pseudorandom rational univariate projection having one zero per zero of i. If n<10 is given, it is the number of digits being used for the pseudorandom coefficients (default: n=5) ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, r = basering NOTE: shows a warning if printlevel>0 (default: printlevel=0) KEYWORDS: rational univariate projection EXAMPLE: example randcharpoly; shows an example" { int pr = printlevel - voice + 2; poly p; poly generic; list l; matrix m; poly q; if (size(#) == 1) { generic = randlinpoly(#[1]); } else { generic = randlinpoly(); } p = reduce(generic,I); m = matmult(p,B,I); q = charpoly(m); dbprint(pr,"*********************************************************************"); dbprint(pr,"* WARNING: This polynomial was obtained using pseudorandom numbers.*"); dbprint(pr,"* If you want to verify the result, please use the command *"); dbprint(pr,"* *"); dbprint(pr,"* verify(p,b,i) *"); dbprint(pr,"* *"); dbprint(pr,"* where p is the polynomial I returned, b is the monomial basis *"); dbprint(pr,"* used, and i the Groebner basis of the ideal *"); dbprint(pr,"*********************************************************************"); return(q); } example { echo = 2; ring r = 0,(x,y,z),dp; ideal i = (x-1)*(x-2),(y-1),(z-1)*(z-2)*(z-3)^2; i = std(i); ideal b = qbase(i); poly p = randcharpoly(b,i); p; nrroots(p); // See nrroots in urrcount.lib int pr = printlevel; printlevel = pr+2; p = randcharpoly(b,i,5); nrroots(p); printlevel = pr; } /////////////////////////////////////////////////////////////////////////////// proc verify(poly p,ideal B,ideal I) "USAGE: verify(p,B,I); p poly, B,I,ideal RETURN: integer: 1 if and only if the polynomial p splits the points of V(I). It's used to check the result of randcharpoly ASSUME: I is given by a Groebner basis and B is an ordered monomial basis of r/I, r = basering NOTE: comments the result if printlevel>0 (default: printlevel=0) SEE ALSO: randcharpoly EXAMPLE: example verify; shows an example" { int pr = printlevel - voice + 2; poly sqr_free; int correct; poly variable; if (isparam(p) || isparam(B) || isparam(I)) { ERROR("This procedure cannot operate with parametric arguments"); } variable = isuni(p); sqr_free = p/gcd(p,diff(p,variable)); correct = (mat_rk(matbil(1,B,I)) == deg(sqr_free)); if (correct) { dbprint(pr,"//Verification successful"); } else { dbprint(pr,"//The choice of random numbers was not useful"); dbprint(pr,"//You might want to try randcharpoly with a larger number of digits"); } return (correct); } example { echo = 2; ring r = 0,(x,y),dp; poly f = x3-xy+y-13+x4-y2x; ideal i = x4-y2x,y2-13; i = std(i); ideal b = qbase(i); poly p = randcharpoly(b,i); verify(p,b,i); } /////////////////////////////////////////////////////////////////////////////// proc randlinpoly(list #) "USAGE: randlinpoly(); randlinpoly(n); n int RETURN: poly: linear combination of the variables of the ring, with pseudorandom coefficients. If n<10 is given, it is the number of digits being used for the range of the coefficients (default: n=5) SEE ALSO: randcharpoly; EXAMPLE: example randlinpoly; shows an example" { int n,i; poly p = 0; int ndigits = 5; if (size(#) == 1) { ndigits = #[1]; } n = nvars(basering); for (i = 1;i <= n;i++) { p = p + var(i)*random(1,10^ndigits); } return (p); } example { echo = 2; ring r = 0,(x,y,z,w),dp; poly p = randlinpoly(); p; randlinpoly(5); } /////////////////////////////////////////////////////////////////////////////// proc powersums(poly f,ideal B,ideal I) "USAGE: powersums(f,b,i); f poly; b,i ideal RETURN: list: the powersums of the results of evaluating f at the zeros of I ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, r = basering SEE ALSO: symmfunc EXAMPLE: example symmfunc; shows an example" { int N,k; list sums; N = size(B); for (k = 1;k <= N;k++) { sums = sums + list(leadcoef(trace(matmult(f^k,B,I)))); } return (sums); } example { echo = 2; ring r = 0,(x,y,z),dp; ideal i = (x-1)*(x-2),(y-1),(z+5); // V(I) = {(1,1,-5),(2,1,-5)} i = std(i); ideal b = qbase(i); poly f = x+y+z; list psums = list(-2-3,4+9); // f evaluated at V(I) gives {-3,-2} list l = powersums(f,b,i); psums; l; } /////////////////////////////////////////////////////////////////////////////// proc symmfunc(list S) "USAGE: symmfunc(s); s list RETURN: list: the symmetric functions of the roots of a polynomial, given the power sums of those roots. SEE ALSO: powersums EXAMPLE: example symmfunc; shows an example" { // Takes the list of power sums and returns the symmetric functions list a; int j,l,N; number sum; N = size(S); a[N+1] = 1; // We set the length of the list and initialize its last element. for (l = N - 1;l >= 0;l--) { sum = 0; for (j = l + 1;j <= N;j++) { sum = sum + ((a[j+1])*(S[j-l])); } sum = -sum; a[l+1] = sum/(N-l); } a = reverse(a); return (a); } example { echo = 2; ring r = 0,x,dp; poly p = (x-1)*(x-2)*(x-3); list psums = list(1+2+3,1+4+9,1+8+27); list l = symmfunc(psums); l; p; // Compare p with the elements of l } /////////////////////////////////////////////////////////////////////////////// proc univarpoly(list l) "USAGE: univarpoly(l); l list RETURN: poly: a polynomial p on the first variable of basering, say x, with p = l[1] + l[2]*x + l[3]*x^2 + ... EXAMPLE: example univarpoly; shows an example" { poly p; int i,n; n = size(l); for (i = 1;i <= n;i++) { p = p + l[i]*var(1)^(n-i); } return (p); } example { echo = 2; ring r = 0,x,dp; list l = list(1,2,3,4,5); poly p = univarpoly(l); p; } /////////////////////////////////////////////////////////////////////////////// proc qbase(ideal i) "USAGE: qbase(I); I zero-dimensional ideal RETURN: ideal: A monomial basis of the quotient between the basering and the ideal I, sorted according to the basering order. SEE ALSO: kbase KEYWORDS: zero-dimensional EXAMPLE: example qbase; shows an example" { ideal b; b = kbase(i); b = reverseideal(sort(b)[1]); // sort sorts in ascending order return (b); } example { echo = 2; ring r = 0,(x,y,z),dp; ideal i = 2x2,-y2,z3; i = std(i); ideal b = qbase(i); b; b = kbase(i); b; // Compare this with the result of qbase } /////////////////////////////////////////////////////////////////////////////// static proc reverseideal(ideal b) // Returns b reversed { int i; ideal result; result = b[1]; for (i = 2;i <= size(b);i++) { result = b[i], result; } return (result); } ///////////////////////////////////////////////////////////////////////////////