// stein01.lib, v 0.01 2005/09/15 //(Arnaud Bodin, last modified 2005/09/30) /////////////////////////////////////////////////////////////////////////////// version="critic01.lib, v 0.01 2005/09/15"; category="Singularities"; // line may be commented for some old version of Singular info=" LIBRARY: stein01.lib Around Stein factorisation of polynomials AUTHORS: Arnaud Bodin : email : Arnaud.Bodin@math.univ-lille1.fr PROCEDURES: proc clean_gaussmatriw(matrix A); clean a Gauss reduced matrix proc findminors(matrix A); find the minors of maximal size proc checkhypo(poly f); check if the algorithm applies proc linearfrompde(poly f); make a linear system (via a pde due to Gao-Ruppert) to find reducible factors proc reducible(poly f); main procedure; give reducible fibers as the zero of a polynomial "; LIB "matrix.lib"; /////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////// //////// Matrix operations ///////// ////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// proc clean_gaussmatrix(matrix A) "USAGE: clean_gaussmatrix(A); A matrix (being gauss-reduced) RETURN: matrix, obtained from A by deleting all non necessary rows and columns EXAMPLE: example clean_gaussmatrix; shows an example " { A = simplify(A,2+4); A = transpose(simplify(transpose(A),2+4)); intvec Li, Lj; int i,j; // A priori we keep all for (i=1; i <= nrows(A); i++) {Li = Li,i;}; for (j=1; j<=ncols(A); j++) {Lj = Lj,j;}; for (j=ncols(A); j>0; j--) { for (i=nrows(A); i>0; i--) { if(A[i,j]==1) {Li = takeoff(Li,i); Lj = takeoff(Lj,j); break;}; if(A[i,j]!=0) {break;}; } } if (size(Li) > 1) { Li = Li[2..size(Li)]; Lj = Lj[2..size(Lj)]; A = submat(A, Li,Lj); A = simplify(A,2+4); A = transpose(simplify(transpose(A),2+4)); } return(A); } example { "EXAMPLE:"; echo = 2; ring r=0,(A,B,C),dp; matrix m[6][8]= 0, 0, 1, 0, B,0, 0,0, 0, 0, -C2, 0, 0,0, 7,0, 0, 0, 0, 0, 0,0, 0,0, 0, 0, B, -A, 0,2A, 0,0, 0, 0, 0, 0, 0,0, 1,0, 0, 0, 0, 0, 0,0, 0,1; print(clean_gaussmatrix(m)); } /////////////////////////////////////////////////////////////////////////////// proc findminors(matrix A) "USAGE: findminors(A); A matrix RETURN: ideal, of minors of maximal size EXAMPLE: example findminors; shows an example " { ideal I; if (ncols(A) > nrows(A)) { I = minor(A,nrows(A)); } else { I = minor(A,ncols(A)); }; I = std(I); return(I); } example { "EXAMPLE:"; echo = 2; ring r=0,(A,B,C),dp; matrix m[6][8]= 0, 0, 1, 0, B,0, 0,0, 0, 1, -C2, 0, 0,0, 7,0, 1, 0, 0, 0, 0,0, 0,0, 0, 0, B, -A, 0,2A, 0,0, 1, 1, 0, 0, 0,0, 1,0, 1, 0, 0, 0, 0,0, 0,1; print(findminors(m)); } /////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////// ////// Linear system from a pde ///// ////////////////////////////////////////// proc checkhypo(poly f) "USAGE: checkhypo(f); f poly RETURN: messages, if hypotheses are not verified EXAMPLE: example checkhypo; shows an example " { def bring = basering; // Number of variables 2 + 1 parameters if (nvars(bring)!= 3) {"Warning : incorrect number of variables !"; "Should be three variables : 2 true variables and 1 for the result ex: (x,y,c).";}; // gcd(f-c,fx)=1 number f0= number(subst(subst(f,x,0),y,0)); ideal J = jacob(f); if ((gcd(f-f0,diff(f,var(1)))!=1) or (dim(std(f))>1)) {"Warning : the polynomial should verified gcd(f-c,fx)=1 for all c !";}; // My warning : should be for all c gcd (f-c, .... // Ok for complex number // // To be changed : in fact it works if some f-c is not reduced // But it doesn't give the right multiplicity // Hence I have to exclude the case of composite polynomials // Characteristic intvec v= 1,0,0; int m = deg(f,v); intvec w = 0,1,1; int n = deg(f, w); if ((char(bring)>0) and (char(bring)<= (2*m-1)*n)) {"Warning : the chararcteristic of the ring is not enough big !"; "char(r) should be zero or greater than (2deg(f,x)-1)deg(f,y)."; }; } example { "EXAMPLE:"; echo = 2; ring r=3, (x,y,c), dp; poly f = x2-y4; print(checkhypo(f)); ring r=0, (x,y,c), dp; poly f = x2*(xy-1)-1; print(checkhypo(f)); } /////////////////////////////////////////////////////////////////////////////// proc linearfrompde(poly f) "USAGE: linearfrompde(f); f poly (ring being r= 0, (x,y,c), dp) RETURN: matrix, whose rank give reducible fibers of f, based on Ruppert-Gao pde EXAMPLE: example linearfrompde; shows an example " { def bring = basering; // Partial degrees intvec v= 1,0,0; int m = deg(f,v); intvec w = 0,1,1; int n = deg(f, w); // New ring with unknowns def bring = basering; // user's ring, of type r=0, (x,y,c), dp; poly theminpoly = minpoly; string thechar = charstr(bring); string thevar = varstr(bring); string rstr = "ring rr = (" + thechar + "), (aa(1..m*(n+1)),bb(1..(m+1)*n),"+ thevar +"), dp;"; execute(rstr); minpoly = number(fetch(bring,theminpoly)); poly ff = imap(bring,f); // f the polynomial in (x,y) poly varx = var(2*m*n+m+n+1); // x variable poly vary = var(2*m*n+m+n+2); // y variable poly varc = var(2*m*n+m+n+3); // c parameter /// g; poly g =0; int i,j; for (i=1; i <= m; i++) { for(j=1; j<= n+1; j++) { g = g + aa(i+(j-1)*m)*varx^(i-1)*vary^(j-1); } }; /// h; poly h =0; for (i=1; i <= m+1; i++) { for(j=1; j<= n; j++) { h = h + bb(i+(j-1)*(m+1))*varx^(i-1)*vary^(j-1); } }; /// PDE Equation; poly E = (ff-varc)*(diff(g,vary)-diff(h,varx))+h*diff(ff,varx)-g*diff(ff,vary); // Make a linear system in aa(1),..., bb(m*n) // Note : not linear in c !! matrix C = coeffs(E,vary); int degy = size(C); ideal EE; for (i=1; i<= degy; i++) { EE = EE + ideal(coeffs(C[i,1],varx)); } // Make a matrix int p = size(EE); int q = 2*m*n+m+n; matrix M[p][q]; for (i=1; i<=p; i++) { for (j=1; j<= m*(n+1); j++) { C = coef(EE[i],aa(j)); if (C[1,1]!=1) {M[i,j] = C[2,1];}; }; for (j=m*(n+1)+1; j<= q; j++) { C = coef(EE[i],bb(j-m*(n+1))); if (C[1,1]!=1) {M[i,j] = C[2,1]; }; }; } // Back to user's ring setring bring; matrix MM; MM = imap(rr,M); return(MM); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,c), dp; poly f = x2-y4; print(linearfrompde(f)); } /////////////////////////////////////////////////////////////////////////////// proc reducible(poly f); main procedure; give reducible fibers as the zero of a polynomial "USAGE: reducible(f); f poly (ring being r= 0, (x,y,c), dp) RETURN: poly, whose zeroes c are corresponds to reducible fibers (f=c) EXAMPLE: example reducible; shows an example " { checkhypo(f); matrix M; M = linearfrompde(f); M = gauss_col(M); M = clean_gaussmatrix(M); // Should be check (by me): M should now really be simplified // otherwise I think it could false ideal I; I = findminors(M); poly g = I[1]; return(g); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,c), dp; poly f = x+y*(x-1)*(x-2)*(x-3)*(x-4); poly g = reducible(f); "Reducible fiber are given by the zeroes of :", g; "After factorization, reducible fibers are :"; factorize(g,1); } ////////////////////////////////////////// ////// Internal procedures ///// ////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// proc takeoff(intvec L, int n); take off n from L, size(L) >1 { int i; intvec LL; for (i=1; i<= size(L); i++) { if (L[i]!= n) {LL = LL,L[i];}; } if (size(LL)>1) {LL = LL[2..size(LL)];}; return(LL); } /////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////// ////// Rational fraction ///// ////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// proc linearfrompdeF(poly p, poly q) "USAGE: linearfrompdeF(polyp, poly q); p/q fraction (ring being r= 0, (x,y,c), dp) RETURN: matrix, whose rank give reducible fibers of p-cq, based on Ruppert-Gao pde EXAMPLE: example linearfrompdeR; shows an example " { def bring = basering; // Partial degrees intvec v= 1,0,0; if (deg(p,v)>deg(q,v)) {int m = deg(p,v);} else {int m = deg(q,v);}; intvec w = 0,1,1; if (deg(p,w)>deg(q,w)) {int n = deg(p,w);} else {int n = deg(q,w);}; // New ring with unknowns def bring = basering; // user's ring, of type r=0, (x,y,c), dp; poly theminpoly = minpoly; string thechar = charstr(bring); string thevar = varstr(bring); string rstr = "ring rr = (" + thechar + "), (aa(1..m*(n+1)),bb(1..(m+1)*n),"+ thevar +"), dp;"; execute(rstr); minpoly = number(fetch(bring,theminpoly)); poly pp = imap(bring,p); // p the polynomial in (x,y) poly qq = imap(bring,q); // q the polynomial in (x,y) poly varx = var(2*m*n+m+n+1); // x variable poly vary = var(2*m*n+m+n+2); // y variable poly varc = var(2*m*n+m+n+3); // c parameter /// g; poly g =0; int i,j; for (i=1; i <= m; i++) { for(j=1; j<= n+1; j++) { g = g + aa(i+(j-1)*m)*varx^(i-1)*vary^(j-1); } }; /// h; poly h =0; for (i=1; i <= m+1; i++) { for(j=1; j<= n; j++) { h = h + bb(i+(j-1)*(m+1))*varx^(i-1)*vary^(j-1); } }; /// PDE Equation; poly E = (pp-varc*qq)*(diff(g,vary)-diff(h,varx))+h*diff(pp-varc*qq,varx)-g*diff(pp-varc*qq,vary); // Make a linear system in aa(1),..., bb(m*n) // Note : not linear in c !! matrix C = coeffs(E,vary); int degy = size(C); ideal EE; for (i=1; i<= degy; i++) { EE = EE + ideal(coeffs(C[i,1],varx)); } // Make a matrix int ip = size(EE); int iq = 2*m*n+m+n; matrix M[ip][iq]; for (i=1; i<=ip; i++) { for (j=1; j<= m*(n+1); j++) { C = coef(EE[i],aa(j)); if (C[1,1]!=1) {M[i,j] = C[2,1];}; }; for (j=m*(n+1)+1; j<= iq; j++) { C = coef(EE[i],bb(j-m*(n+1))); if (C[1,1]!=1) {M[i,j] = C[2,1]; }; }; } // Back to user's ring setring bring; matrix MM; MM = imap(rr,M); return(MM); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,c), dp; poly p = x3-y3; poly q = 1-y^3; print(linearfrompdeF(f)); } /////////////////////////////////////////////////////////////////////////////// proc reducibleF(poly p, poly q); main procedure; give reducible fibers as the zero of a polynomial "USAGE: reducibleF(p,q); p,q poly (ring being r= 0, (x,y,c), dp) RETURN: poly, whose zeroes c are corresponds to reducible fibers (p-cq=0) EXAMPLE: example reducibleF; shows an example " { // checkhypo(f); matrix M; M = linearfrompdeF(p,q); M = gauss_col(M); M = clean_gaussmatrix(M); ideal I; print(M); I = findminors(M); poly g = I[1]; return(g); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,c), dp; poly p = x3-y3; poly q = 1-y3; poly g = reducible(p,q); "Reducible fiber are given by the zeroes of :", g; "After factorization, reducible fibers are :"; factorize(g,1); }