version="$Id: equising.lib,v 1.7.2.5 2003/02/25 16:22:38 lossen Exp $"; category="Singularities"; info=" LIBRARY: equising.lib Equisingularity Stratum of a Family of Plane Curves AUTHOR: Christoph Lossen, lossen@mathematik.uni-kl.de Andrea Mindnich, mindnich@mathematik.uni-kl.de MAIN PROCEDURES: tau_es(f); codim of mu-const stratum in semi-universal def. base esIdeal(f); (Wahl's) equisingularity ideal of f esStratum(F[,m,L]); equisingularity stratum of a family F isEquising(F[,m,L]); tests if a given deformation is equisingular AUXILIARY PROCEDURE: control_Matrix(M); computes list of blowing-up data "; LIB "hnoether.lib"; LIB "poly.lib"; LIB "elim.lib"; LIB "deform.lib"; LIB "sing.lib"; //////////////////////////////////////////////////////////////////////////////// // // The following (static) procedures are used by esComputation // //////////////////////////////////////////////////////////////////////////////// // COMPUTES a weight vector. x and y get weight 1 and all other // variables get weight 0. static proc xyVector() { intvec iv ; iv[nvars(basering)]=0 ; iv[rvar(x)] =1; iv[rvar(y)] =1; return (iv); } //////////////////////////////////////////////////////////////////////////////// // exchanges the variables x and y in the polynomial f static proc swapXY(poly f) { def r_base = basering; ideal MI = maxideal(1); MI[rvar(x)]=y; MI[rvar(y)]=x; map phi = r_base, MI; f=phi(f); return (f); } //////////////////////////////////////////////////////////////////////////////// // computes m-jet w.r.t. the variables x,y (other variables weighted 0 static proc m_Jet(poly F,int m); { intvec w=xyVector(); poly Fd=jet(F,m,w); return(Fd); } //////////////////////////////////////////////////////////////////////////////// // computes the 4 control matrices (input is multsequence(L)) proc control_Matrix(list M); "USAGE: control_Matrix(L); L list ASSUME: L is the output of multsequence(reddevelop(f)). RETURN: list M of 4 intmat's: @format M[1] contains the multiplicities at the respective infinitely near points p[i,j] (i=step of blowup+1, j=branch) -- if branches j=k,...,k+m pass through the same p[i,j] then the multiplicity is stored in M[1][k,j], while M[1][k+1]=...=M[1][k+m]=0. M[2] contains the number of branches meeting at p[i,j] (again, the information is stored according to the above rule) M[3] contains the information about the splitting of M[1][i,j] with respect to different tangents of branches at p[i,j] (information is stored only for minimal j>=k corresponding to a new tangent direction). The entries are the sum of multiplicities of all branches with the respective tangent. M[4] contains the maximal sum of higher multiplicities for a branch passing through p[i,j] ( = degree Bound for blowing up) @end format NOTE: the branches are ordered in such a way that only consecutive branches can meet at an infinitely near point. @* the final rows of the matrices M[1],...,M[3] is (1,1,1,...,1), and correspond to infinitely near points such that the strict transforms of the branches are smooth and intersect the exceptional divisor transversally. SEE ALSO: multsequence EXAMPLE: example control_Matrix; shows an example " { int i,j,k,dummy; dummy=0; for (j=1;j<=ncols(M[2]);j++) { dummy=dummy+M[1][nrows(M[1])-1,j]-M[1][nrows(M[1]),j]; } intmat S[nrows(M[1])+dummy][ncols(M[1])]; intmat T[nrows(M[1])+dummy][ncols(M[1])]; intmat U[nrows(M[1])+dummy][ncols(M[1])]; intmat maxDeg[nrows(M[1])+dummy][ncols(M[1])]; for (i=1;i<=nrows(M[2]);i++) { dummy=1; for (j=1;j<=ncols(M[2]);j++) { for (k=dummy;k1) { U[i-1,dummy]=U[i-1,dummy]+M[1][i-1,k]; } } dummy=k; } } // adding an extra row (in some cases needed to control ES-Stratum // computation) for (i=nrows(M[1]);i<=nrows(S);i++) { for (j=1;j<=ncols(M[2]);j++) { S[i,j]=1; T[i,j]=1; U[i,j]=1; } } // Computing the degree Bounds to be stored in M[4]: for (i=1;i<=nrows(S);i++) { dummy=1; for (j=1;j<=ncols(S);j++) { for (k=dummy;k=2;i--) { for (j=1;j<=ncols(S);j++) { maxDeg[i-1,j]=maxDeg[i-1,j]+maxDeg[i,j]; } } list L=S,T,U,maxDeg; return(L); } //////////////////////////////////////////////////////////////////////////////// // matrix of higher tangent directions: // returns list: 1) tangent directions // 2) swapping information (x <--> y) static proc inf_Tangents(list L,int s); // L aus reddevelop, { matrix M; matrix B[s][size(L)]; intvec V; intmat Mult=multsequence(L)[1]; int i,j,k,counter,e; for (k=1;k<=size(L);k++) { V[k]=L[k][3]; // switch: 0 --> tangent 2nd parameter // 1 --> tangent 1st parameter e=0; M=L[k][1]; counter=1; B[counter,k]=M[1,1]; for (i=1;i<=nrows(M);i++) { for (j=2;j<=ncols(M);j++) { counter=counter+1; if (M[i,j]==x) { if (i<>nrows(M)) { B[counter,k]=M[i,j]; j=ncols(M)+1; // goto new row of HNmatrix... if (counter<>s) { if (counter+1<=nrows(Mult)) { e=Mult[counter-1,k]-Mult[counter,k]-Mult[counter+1,k]; } else { e=Mult[counter-1,k]-Mult[counter,k]-1; } } } else { B[counter,k]=0; j=ncols(M)+1; // goto new row of HNmatrix... } } else { if (e<=0) { B[counter,k]=M[i,j]; } else // point is still proximate to an earlier point { B[counter,k]=y; // marking proximity (without swap....) if (counter+1<=nrows(Mult)) { e=e-Mult[counter+1,k]; } else { e=e-1; } } } if (counter==s) // given number of points determined { j=ncols(M)+1; i=nrows(M)+1; // leave procedure } } } } L=B,V; return(L); } //////////////////////////////////////////////////////////////////////////////// // compute "good" upper bound for needed number of help variables // static proc Determine_no_b(intmat U,matrix B) // U is assumed to be 3rd output of control_Matrix // B is assumed to be 1st output of inf_Tangents { int i,j,counter; for (j=1;j<=ncols(U);j++) { for (i=1;i<=nrows(U);i++) { if (U[i,j]>1) { if (B[i,j]<>x and B[i,j]<>y) { counter=counter+1; } } } } counter=counter+ncols(U); return(counter); } //////////////////////////////////////////////////////////////////////////////// // compute number of infinitely near free points corresponding to non-zero // entries in control_Matrix[1] (except first row) // static proc no_freePoints(intmat Mult,matrix B) // Mult is assumed to be 1st output of control_Matrix // U is assumed to be 3rd output of control_Matrix // B is assumed to be 1st output of inf_Tangents { int i,j,k,counter; for (j=1;j<=ncols(Mult);j++) { for (i=2;i<=nrows(Mult);i++) { if (Mult[i,j]>=1) { if (B[i-1,j]<>x and B[i-1,j]<>y) { counter=counter+1; } } } } return(counter); } /////////////////////////////////////////////////////////////////////////////// // COMPUTES string(minpoly) and substitutes the parameter by newParName static proc makeMinPolyString (string newParName) { int i; string parName = parstr(basering); int parNameSize = size(parName); string oldMinPolyStr = string (minpoly); int minPolySize = size(oldMinPolyStr); string newMinPolyStr = ""; for (i=1;i <= minPolySize; i++) { if (oldMinPolyStr[i,parNameSize] == parName) { newMinPolyStr = newMinPolyStr + newParName; i = i + parNameSize-1; } else { newMinPolyStr = newMinPolyStr + oldMinPolyStr[i]; } } return(newMinPolyStr); } /////////////////////////////////////////////////////////////////////////////// // // DEFINES: A new basering, "myRing", // with new names for the parameters and variables. // The new names for the parameters are a(1..k), // and t(1..s),x,y for the variables // The ring ordering is ordStr. // NOTE: This proc uses 'execute'. static proc createMyRing_new(poly p_F, string ordStr, string minPolyStr, int no_b) { def r_old = basering; int chara = char(basering); string charaStr; int i; string helpStr; int nDefParams = nvars(r_old)-2; ideal qIdeal = ideal(basering); if (npars(basering) == 0 and minPolyStr == "") { helpStr = "ring myRing1 =" + string(chara)+ ", (t(1..nDefParams), x, y),("+ ordStr +");"; execute(helpStr); } else { charaStr = charstr(basering); if (charaStr == string(chara) + "," + parstr(basering) or minPolyStr<>"") { if (minPolyStr<>"") { helpStr = "ring myRing1 = (" + string(chara) + ",a), (t(1..nDefParams), x, y),(" + ordStr + ");"; execute(helpStr); execute (minPolyStr); } else // no minpoly given { helpStr = "ring myRing1 = (" + string(chara) + ",a(1..npars(basering)) ), (t(1..nDefParams), x, y),(" + ordStr + ");"; execute(helpStr); } } else { // ground field is of type (p^k,a).... i = find (charaStr,","); helpStr = "ring myRing1 = (" + charaStr[1,i] + "a), (t(1..nDefParams), x, y),(" + ordStr + ");"; execute (helpStr); } } ideal qIdeal = fetch(r_old, qIdeal); poly p_F = fetch(r_old, p_F); // Extension by no_b auxiliary variables if (no_b>0) { if (npars(basering) == 0) { ordStr = "(dp("+string(no_b)+"),"+ordStr+")"; helpStr = "ring myRing =" + string(chara)+ ", (b(1..no_b), t(1..nDefParams), x, y)," + ordStr +";"; execute(helpStr); } else { charaStr = charstr(basering); if (charaStr == string(chara) + "," + parstr(basering)) { if (minpoly !=0) { ordStr = "(dp(" + string(no_b) + ")," + ordStr + ")"; minPolyStr = makeMinPolyString("a"); helpStr = "ring myRing = (" + string(chara) + ",a), (b(1..no_b), t(1..nDefParams), x, y)," + ordStr + ";"; execute(helpStr); helpStr = "minpoly =" + minPolyStr + ";"; execute (helpStr); } else // no minpoly given { ordStr = "(dp(" + string(no_b) + ")," + ordStr + ")"; helpStr = "ring myRing = (" + string(chara) + ",a(1..npars(basering)) ), (b(1..no_b), t(1..nDefParams), x, y)," + ordStr + ";"; execute(helpStr); } } else { i = find (charaStr,","); ordStr = "(dp(" + string(no_b) + ")," + ordStr + ")"; helpStr = "ring myRing = (" + charaStr[1,i] + "a), (b(1..no_b), t(1..nDefParams), x, y)," + ordStr + ";"; execute (helpStr); } } ideal qIdeal = imap(myRing1, qIdeal); if(qIdeal != 0) { def r_base = basering; kill myRing; qring myRing = std(qIdeal); } poly p_F = imap(myRing1, p_F); kill myRing1; } else { if(qIdeal != 0) { def r_base = basering; kill myRing1; qring myRing1 = std(qIdeal); } def myRing=myRing1; kill myRing1; } keepring(myRing); } //////////////////////////////////////////////////////////////////////////////// // returns list of coef, leadmonomial // static proc determine_coef (poly Fm) { def r_base = basering; // is assumed to be the result of createMyRing int chara = char(basering); string charaStr; int i; string minPolyStr = ""; string helpStr = ""; if (npars(basering) == 0) { helpStr = "ring myRing1 =" + string(chara)+ ", (y,x),ds;"; execute(helpStr); } else { charaStr = charstr(basering); if (charaStr == string(chara) + "," + parstr(basering)) { if (minpoly !=0) { minPolyStr = makeMinPolyString("a"); helpStr = "ring myRing1 = (" + string(chara) + ",a), (y,x),ds;"; execute(helpStr); helpStr = "minpoly =" + minPolyStr + ";"; execute (helpStr); } else // no minpoly given { helpStr = "ring myRing1 = (" + string(chara) + ",a(1..npars(basering)) ), (y,x),ds;"; execute(helpStr); } } else { i = find (charaStr,","); helpStr = " ring myRing1 = (" + charaStr[1,i] + "a), (y,x),ds;"; execute (helpStr); } } poly f=imap(r_base,Fm); poly g=leadmonom(f); setring r_base; poly g=imap(myRing1,g); kill myRing1; def M=coef(Fm,xy); for (i=1; i<=ncols(M); i++) { if (M[1,i]==g) { poly h=M[2,i]; // determine coefficient of leading monomial (in K[t]) i=ncols(M)+1; } } return(list(h,g)); } //////////////////////////////////////////////////////////////////////////////// // Defines a new ring without deformation-parameters. static proc createHNERing() { string str; string minpolyStr = string(minpoly); str = " ring HNERing = (" + charstr(basering) + "), (x,y), ls;"; execute (str); str = "minpoly ="+ minpolyStr+";"; execute(str); keepring(HNERing); } /////////////////////////////////////////////////////////////////////////////// // RETURNS: 1, if p_f = 0 or char(basering) divides the order of p_f // or p_f is not squarefree. // 0, otherwise static proc checkPoly (poly p_f) { int i_print = printlevel - voice + 3; int i_ord; if (p_f == 0) { print("Input is a 'deformation' of the zero polynomial!"); return(1); } i_ord = mindeg1(p_f); if (number(i_ord) == 0) { print("Characteristic of coefficient field " +"divides order of zero-fiber !"); return(1); } if (squarefree(p_f) != p_f) { print("Original polynomial (= zero-fiber) is not reduced!"); return(1); } return(0); } //////////////////////////////////////////////////////////////////////////////// static proc make_ring_small(ideal J) // returns varstr for new ring, the map and the number of vars { attrib(J,"isSB",1); int counter=0; ideal newmap; string newvar=""; for (int i=1; i<=nvars(basering); i++) { if (reduce(var(i),J)<>0) { newmap[i]=var(i); if (newvar=="") { newvar=newvar+string(var(i)); counter=counter+1; } else { newvar=newvar+","+string(var(i)); counter=counter+1; } } else { newmap[i]=0; } } list L=newvar,newmap,counter; attrib(J,"isSB",0); return(L); } /////////////////////////////////////////////////////////////////////////////// // The following procedure is called by esStratum (typ=0), resp. by // isEquising (typ=1) /////////////////////////////////////////////////////////////////////////////// static proc esComputation (int typ, poly p_F, list #) { // Initialize variables int branch=1; int blowup=1; int auxVar=1; int nVars; intvec upper_bound, upper_bound_old, fertig, soll; list blowup_string; int i_print= printlevel-voice+2; int no_b, number_of_branches, swapped; int i,j,k,m, counter, dummy; string helpStr = ""; string ordStr = ""; string MinPolyStr = ""; if (nvars(basering)<=2) { print("family is trivial (no deformation parameters)!"); if (typ==1) //isEquising { return(1); } else { return(list(ideal(0),0)); } } if (size(#)>0) { if (typeof(#[1])=="int") { def artin_bd=#[1]; // compute modulo maxideal(artin_bd) if (artin_bd <= 1) { print("Do you really want to compute over Basering/maxideal(" +string(artin_bd)+") ?"); print("No computation performed !"); if (typ==1) //isEquising { return(1); } else { return(list(ideal(0),int(1))); } } if (size(#)>1) { if (typeof(#[2])=="list") { def @L=#[2]; // is assumed to be the Hamburger-Noether matrix } } } if (typeof(#[1])=="list") { def @L=#[1]; // is assumed to be the Hamburger-Noether matrix } } def old_ring=basering; if (defined(HNEring)) { def @HNEring = HNEring; kill HNEring; } if(defined(@L)<=0) { createHNERing(); // Basering changed to HNERing (variables x,y, with ls ordering) k=nvars(old_ring); matrix Map_Phi[1][k]; Map_Phi[1,k-1]=x; Map_Phi[1,k]=y; map phi=old_ring,Map_Phi; poly f=phi(p_F); // Heuristics: if x,y are transversal parameters then computation of HNE // can be much faster when exchanging variables...! if (2*size(coeffs(f,x))0) { setring old_ring; return(0); } } F(1)=p_F; // and reduce the remaining terms in F(1): bl_Map(1)=maxideal(1); attrib(J,"isSB",1); bl_Map(1)=reduce(bl_Map(1),J); attrib(J,"isSB",0); phi=myRing,bl_Map(1); F(1)=phi(F(1)); // simplify F(1) attrib(J,"isSB",1); F(1)=reduce(F(1),J); attrib(J,"isSB",0); // now we compute the m-jet: Fm=m_Jet(F(1),m); G=1; counter=branch; k=upper_bound[branch]; F_save=F(1); // is truncated differently in the following loop while(counter<=k) { F(counter)=m_Jet(F_save,maxDeg[blowup,counter]); if (V[counter]==0) // 2nd ring variable is tangent to this branch { G=G*(y-(b(auxVar)+B[blowup,counter])*x)^(M[3][blowup,counter]); } else // 1st ring variable is tangent to this branch { G=G*(x-(b(auxVar)+B[blowup,counter])*y)^(M[3][blowup,counter]); F(counter)=swapXY(F(counter)); } bl_Map(counter)=maxideal(1); bl_Map(counter)[nvars(basering)]=xy+(b(auxVar)+B[blowup,counter])*x; auxVar=auxVar+1; upper_bound[counter]=counter+M[2][blowup+1,counter]-1; counter=counter+M[2][blowup+1,counter]; } list LeadDataFm=determine_coef(Fm); def LeadDataG=coef(G,xy); for (i=1; i<=ncols(LeadDataG); i++) { if (LeadDataG[1,i]==LeadDataFm[2]) { poly LeadG = LeadDataG[2,i]; // determine the coefficient of G i=ncols(LeadDataG)+1; } } G=LeadDataFm[1]*G-LeadG*Fm; // leading terms in y should cancel... coef_Mat = coef(G,xy); Jnew=coef_Mat[2,1..ncols(coef_Mat)]; if (defined(artin_bd)) // the artin_bd-th power of the maxideal of // deformation parameters can be cutted off { Jnew=jet(Jnew,artin_bd-1); } // simplification of Jnew Jnew=interred(Jnew); J=J,Jnew; if (typ==1) // isEquising { if(ideal(nselect(J,1,no_b))<>0) { setring old_ring; return(0); } } while (fertig<>soll and blowup1) { if (M[3][blowup-1,branch]==1 and ((B[blowup,branch]<>x and B[blowup,branch]<>y) or (blowup==nrows(M[3])) )) { fertig[branch]=1; dbprint(i_print,"// 1 branch finished"); } } if (branch<=upper_bound_old[branch] and fertig[branch]<>1) { for (i=branch;i>=1;i--) { if (M[1][blowup-1,i]<>0) { m=M[1][blowup-1,i]; // multiplicity before blowup i=0; } } // we blow up the branch and take the strict transform: attrib(J,"isSB",1); bl_Map(branch)=reduce(bl_Map(branch),J); attrib(J,"isSB",0); phi=myRing,bl_Map(branch); F(branch)=phi(F(branch))/x^m; // simplify F attrib(Jnew,"isSB",1); F(branch)=reduce(F(branch),Jnew); attrib(Jnew,"isSB",0); m=M[1][blowup,branch]; // multiplicity after blowup Fm=m_Jet(F(branch),m); // homogeneous part of lowest degree // we check for Fm=F[k]*...*F[k+s] where // // F[j]=(y-b'(j)*x)^m(j), respectively F[j]=(-b'(j)*y+x)^m(j) // // according to the entries m(j)= M[3][blowup,j] and // b'(j) mod m_A = B[blowup,j] // computed from the HNE of the special fibre of the family: G=1; counter=branch; k=upper_bound[branch]; F_save=F(branch); while(counter<=k) { F(counter)=m_Jet(F_save,maxDeg[blowup,counter]); if (B[blowup,counter]<>x and B[blowup,counter]<>y) { G=G*(y-(b(auxVar)+B[blowup,counter])*x)^(M[3][blowup,counter]); bl_Map(counter)=maxideal(1); bl_Map(counter)[nvars(basering)]= xy+(b(auxVar)+B[blowup,counter])*x; auxVar=auxVar+1; } else { if (B[blowup,counter]==x) { G=G*x^(M[3][blowup,counter]); // branch has tangent x !! F(counter)=swapXY(F(counter)); // will turn x to y for blow up bl_Map(counter)=maxideal(1); bl_Map(counter)[nvars(basering)]=xy; } else { G=G*y^(M[3][blowup,counter]); // tangent has to be y bl_Map(counter)=maxideal(1); bl_Map(counter)[nvars(basering)]=xy; } } upper_bound[counter]=counter+M[2][blowup+1,counter]-1; counter=counter+M[2][blowup+1,counter]; } G=determine_coef(Fm)[1]*G-Fm; // leading terms in y should cancel coef_Mat = coef(G,xy); Jnew=coef_Mat[2,1..ncols(coef_Mat)]; if (defined(artin_bd)) // the artin_bd-th power of the maxideal of // deformation parameters can be cutted off { Jnew=jet(Jnew,artin_bd-1); } // simplification of J Jnew=interred(Jnew); J=J,Jnew; if (typ==1) // isEquising { if(ideal(nselect(J,1,no_b))<>0) { setring old_ring; return(0); } } } } if (number_of_branches>=2) { J=interred(J); if (typ==1) // isEquising { if(ideal(nselect(J,1,no_b))<>0) { setring old_ring; return(0); } } } } dbprint(i_print,"// "); dbprint(i_print,"// Elimination starts:"); dbprint(i_print,"// -------------------"); poly gg; int b_left=no_b; for (i=1;i<=no_b; i++) { attrib(J,"isSB",1); gg=reduce(b(i),J); if (gg==0) { b_left = b_left-1; // another b(i) has to be 0 } J = subst(J, b(i), gg); attrib(J,"isSB",0); } J=simplify(J,10); if (typ==1) // isEquising { if(ideal(nselect(J,1,no_b))<>0) { setring old_ring; return(0); } } ideal J_no_b = nselect(J,1,no_b); if (size(J) > size(J_no_b)) { dbprint(i_print,"// std computation started"); // some b(i) didn't appear in linear conditions and have to be eliminated if (defined(artin_bd)) { // first we make the ring smaller (removing variables, which are // forced to 0 by J list LL=make_ring_small(J); ideal Shortmap=LL[2]; minPolyStr = ""; if (minpoly !=0) { minPolyStr = "minpoly = "+string(minpoly); } ordStr = "dp(" + string(b_left) + "),dp"; ideal qId = ideal(basering); helpStr = "ring Shortring = (" + charstr(basering) + "),("+ LL[1] +") , ("+ ordStr +");"; execute(helpStr); execute(minPolyStr); // ring has changed to "Shortring" ideal MM=maxideal(artin_bd); MM=subst(MM,x,0); MM=subst(MM,y,0); MM=simplify(MM,2); dbprint(i_print-1,"// maxideal("+string(artin_bd)+") has " +string(size(MM))+" elements"); dbprint(i_print-1,"//"); // we change to the qring mod m^artin_bd // first, we have to check if we were in a qring when starting ideal qId = imap(myRing, qId); if (qId == 0) { attrib(MM,"isSB",1); qring QQ=MM; } else { qId=qId,MM; qring QQ = std(qId); } ideal Shortmap=imap(myRing,Shortmap); map phiphi=myRing,Shortmap; ideal J=phiphi(J); J=std(J); J=nselect(J,1,no_b); setring myRing; // back to "myRing" J=nselect(J,1,no_b); Jnew=imap(QQ,J); J=J,Jnew; J=interred(J); } else { J=std(J); J=nselect(J,1,no_b); } } dbprint(i_print,"// finished"); dbprint(i_print,"// "); minPolyStr = ""; if (minpoly !=0) { minPolyStr = "minpoly = "+string(minpoly); } // remove the ring created by reddevelop kill HNEring; if (defined(@HNEring)) { def HNEring=@HNEring; export(HNEring); } if (typ==1) // isEquising { if(J<>0) { setring old_ring; return(0); } else { setring old_ring; return(1); } } setring old_ring; // we are back in the original ring if (npars(myRing)<>0) { ideal qIdeal = ideal(basering); helpStr = "ring ESSring = (" + string(char(basering))+ "," + parstr(myRing) + ") , ("+ varstr(basering)+") , ("+ ordstr(basering) +");"; execute(helpStr); execute(minPolyStr); // basering has changed to ESSring ideal qIdeal = fetch(old_ring, qIdeal); if(qIdeal != 0) { def r_base = basering; kill ESSring; qring ESSring = std(qIdeal); } kill qIdeal; ideal SSS; for (int ii=1;ii<=nvars(basering);ii++) { SSS[ii+no_b]=var(ii); } map phi=myRing,SSS; // b(i) variables are mapped to zero ideal ES=phi(J); kill phi; if (defined(p_F)<=0) { poly p_F=fetch(old_ring,p_F); export(p_F); } export(ES); setring old_ring; dbprint(i_print+1," // 'esStratum' created a list M of a ring and an integer. // To access the ideal defining the equisingular stratum, type: def ESSring = M[1]; setring ESSring; ES; "); return(list(ESSring,0)); } else { // no new ring definition necessary ideal SSS; for (int ii=1;ii<=nvars(basering);ii++) { SSS[ii+no_b]=var(ii); } map phi=myRing,SSS; // b(i) variables are mapped to zero ideal ES=phi(J); kill phi; setring old_ring; dbprint(i_print,"// output of 'esStratum' is list consisting of: // _[1] = ideal defining equisingular stratum // _[2] = 0"); return(list(ES,0)); } } //////////////////////////////////////////////////////////////////////////////// proc tau_es (poly f,list #) "USAGE: tau_es(f); f poly ASSUME: f is a reduced bivariate polynomial, the basering has precisely two variables, is local and no qring. RETURN: int, the codimension of the mu-const stratum in the semi-universal deformation base. NOTE: printlevel>=1 displays additional information. When called with any additional parameter, the computation of the Milnor number is avoided (no check for NND). SEE ALSO: esIdeal, tjurina, invariants EXAMPLE: example tau_es; shows an example. " { int i,j,k,s; int slope_x, slope_y, upper; int i_print = printlevel - voice + 3; string MinPolyStr; // some checks first if ( nvars(basering)<>2 ) { print("// basering has not the correct number (two) of variables !"); print("// computation stopped"); return(0); } if ( mult(std(1+var(1)+var(2))) <> 0) { print("// basering is not local !"); print("// computation stopped"); return(0); } if (mult(std(f))<=1) { // f is rigid return(0); } if ( deg(squarefree(f))!=deg(f) ) { print("// input polynomial was not reduced"); print("// try squarefree(f); first"); return(0); } def old_ring=basering; execute("ring @myRing=("+charstr(basering)+"),("+varstr(basering)+"),ds;"); poly f=imap(old_ring,f); ideal Jacobi_Id = jacob(f); // check for A_k singularity // ---------------------------------------- if (mult(std(f))==2) { dbprint(i_print-1,"// "); dbprint(i_print-1,"// polynomial defined A_k singularity"); dbprint(i_print-1,"// "); return( vdim(std(Jacobi_Id)) ); } // check for D_k singularity // ---------------------------------------- if (mult(std(f))==3 and size(factorize(jet(f,3))[1])>=3) { dbprint(i_print,"// "); dbprint(i_print,"// polynomial defined D_k singularity"); dbprint(i_print,"// "); ideal ES_Id = f, jacob(f); return( vdim(std(Jacobi_Id))); } if (size(#)==0) { // check if Newton polygon non-degenerated // ---------------------------------------- Jacobi_Id=std(Jacobi_Id); int mu = vdim(Jacobi_Id); poly f_tilde=f+var(1)^mu+var(2)^mu; //to obtain convenient Newton-polygon list NP=newtonpoly(f_tilde); dbprint(i_print-1,"// Newton polygon:"); dbprint(i_print-1,NP); dbprint(i_print-1,""); if(is_NND(f,mu,NP)) // f is Newton non-degenerate { upper=NP[1][2]; ideal ES_Id= x^k*y^upper; dbprint(i_print-1,"polynomial is Newton non-degenerated"); dbprint(i_print-1,""); k=0; for (i=1;i<=size(NP)-1;i++) { slope_x=NP[i+1][1]-NP[i][1]; slope_y=NP[i][2]-NP[i+1][2]; for (k=NP[i][1]+1; k<=NP[i+1][1]; k++) { while ( slope_x*upper + slope_y*k >= slope_x*NP[i][2] + slope_y*NP[i][1]) { upper=upper-1; } upper=upper+1; ES_Id=ES_Id, x^k*y^upper; } } ES_Id=std(ES_Id); dbprint(i_print-2,"ideal of monomials above Newton bd. is generated by:"); dbprint(i_print-2,ES_Id); ideal ESfix_Id = ES_Id, f, maxideal(1)*jacob(f); ES_Id = ES_Id, Jacobi_Id; ES_Id = std(ES_Id); dbprint(i_print-1,"// "); dbprint(i_print-1,"// Equisingularity ideal is computed!"); dbprint(i_print-1,""); return(vdim(ES_Id)); } else { dbprint(i_print-1,"polynomial is Newton degenerated !"); dbprint(i_print-1,""); } } // for Newton degenerate polynomials, we compute the HN-development, and // count the number of free points ..... if (defined(HNEring)) { def @HNEring = HNEring; kill HNEring; } dbprint(i_print-1,"// "); dbprint(i_print-1,"// Compute HN development"); dbprint(i_print-1,"// ----------------------"); i=printlevel; printlevel=printlevel-5; if (2*size(coeffs(f,x))2 ) { print("// basering has not the correct number (two) of variables !"); print("// computation stopped"); return(list(0,0)); } if ( mult(std(1+var(1)+var(2))) <> 0) { print("// basering is not local !"); print("// computation stopped"); return(list(0,0)); } if (mult(std(f))<=1) { // f is rigid return(list(ideal(1),ideal(1))); } if ( deg(squarefree(f))!=deg(f) ) { print("// input polynomial was not squarefree"); print("// try squarefree(f); first"); return(list(0,0)); } if (char(basering)<>0) { if (mult(std(f)) mod char(basering)==0) { print("// characteristic of ground field divides " + "multiplicity of polynomial !"); print("// computation stopped"); return(list(0,0)); } } // check for A_k singularity // ---------------------------------------- if (mult(std(f))==2) { dbprint(i_print,"// "); dbprint(i_print,"// polynomial defined A_k singularity"); dbprint(i_print,"// "); ideal ES_Id = f, jacob(f); ES_Id = interred(ES_Id); ideal ESfix_Id = f, maxideal(1)*jacob(f); ESfix_Id= interred(ESfix_Id); return( list(ES_Id,ESfix_Id) ); } // check for D_k singularity // ---------------------------------------- if (mult(std(f))==3 and size(factorize(jet(f,3))[1])>=3) { dbprint(i_print,"// "); dbprint(i_print,"// polynomial defined D_k singularity"); dbprint(i_print,"// "); ideal ES_Id = f, jacob(f); ES_Id = interred(ES_Id); ideal ESfix_Id = f, maxideal(1)*jacob(f); ESfix_Id= interred(ESfix_Id); return( list(ES_Id,ESfix_Id) ); } // check if Newton polygon non-degenerated // ---------------------------------------- int mu = milnor(f); poly f_tilde=f+var(1)^mu+var(2)^mu; //to obtain a convenient Newton-polygon list NP=newtonpoly(f_tilde); dbprint(i_print-1,"// Newton polygon:"); dbprint(i_print-1,NP); dbprint(i_print-1,""); if(is_NND(f,mu,NP)) // f is Newton non-degenerate { upper=NP[1][2]; ideal ES_Id= x^k*y^upper; dbprint(i_print,"polynomial is Newton non-degenerated"); dbprint(i_print,""); k=0; for (i=1;i<=size(NP)-1;i++) { slope_x=NP[i+1][1]-NP[i][1]; slope_y=NP[i][2]-NP[i+1][2]; for (k=NP[i][1]+1; k<=NP[i+1][1]; k++) { while ( slope_x*upper + slope_y*k >= slope_x*NP[i][2] + slope_y*NP[i][1]) { upper=upper-1; } upper=upper+1; ES_Id=ES_Id, x^k*y^upper; } } ES_Id=std(ES_Id); dbprint(i_print-1,"ideal of monomials above Newton bd. is generated by:"); dbprint(i_print-1,ES_Id); ideal ESfix_Id = ES_Id, f, maxideal(1)*jacob(f); ES_Id = ES_Id, f, jacob(f); dbprint(i_print,"// "); dbprint(i_print,"// equisingularity ideal is computed!"); return(list(ES_Id,ESfix_Id)); } else { dbprint(i_print,"polynomial is Newton degenerated !"); dbprint(i_print,""); } def old_ring=basering; dbprint(i_print,"// "); dbprint(i_print,"// versal deformation with triv. section"); dbprint(i_print,"// ====================================="); dbprint(i_print,"// "); ideal JJ=maxideal(1)*jacob(f); ideal kbase_versal=kbase(std(JJ)); s=size(kbase_versal); string ring_versal="ring @Px = ("+charstr(basering)+"),(t(1.."+string(s)+")," +varstr(basering)+"),(ds("+string(s)+")," +ordstr(basering)+");"; MinPolyStr = string(minpoly); execute(ring_versal); if (MinPolyStr<>"0") { MinPolyStr = "minpoly="+MinPolyStr; execute(MinPolyStr); } // basering has changed to @Px poly F=imap(old_ring,f); ideal kbase_versal=imap(old_ring,kbase_versal); for (i=1; i<=s; i++) { F=F+var(i)*kbase_versal[i]; } dbprint(i_print-1,F); dbprint(i_print-1,""); ideal ES_Id; dbprint(i_print,"// "); dbprint(i_print,"// Compute equisingular Stratum over Spec(C[t]/t^2)"); dbprint(i_print,"// ================================================"); dbprint(i_print,"// "); list M=esStratum(F,2); dbprint(i_print,"// finished"); dbprint(i_print,"// "); if (M[2]==1) // error occured during esStratum computation { print("Some error has occured during the computation"); return(list(0,0)); } if ( typeof(M[1])=="ideal" ) { int defpars = nvars(basering)-2; poly Fred; poly g; F=reduce(F,std(M[1])); for (i=1; i<=defpars; i++) { Fred=reduce(F,std(var(i))); g=subst(F-Fred,var(i),1); ES_Id=ES_Id, g; F=Fred; } setring old_ring; // back to original ring ideal ES_Id = imap(@Px,ES_Id); ES_Id = interred(ES_Id); ideal ESfix_Id = ES_Id, f, maxideal(1)*jacob(f); ES_Id = ES_Id, f, jacob(f); return(list(ES_Id,ESfix_Id)); } else { def AuxRing=M[1]; dbprint(i_print,"// "); dbprint(i_print,"// change ring to ESSring"); setring AuxRing; // contains p_F, ES int defpars = nvars(basering)-2; poly Fred; poly g; ideal ES_Id; p_F=reduce(p_F,std(ES)); for (i=1; i<=defpars; i++) { Fred=reduce(p_F,std(var(i))); g=subst(p_F-Fred,var(i),1); ES_Id=ES_Id, g; p_F=Fred; } dbprint(i_print,"// "); dbprint(i_print,"// back to the original ring"); setring old_ring; // back to original ring ideal ES_Id = imap(AuxRing,ES_Id); ES_Id = interred(ES_Id); kill @Px; kill AuxRing; ideal ESfix_Id = ES_Id, f, maxideal(1)*jacob(f); ES_Id = ES_Id, f, jacob(f); dbprint(i_print,"// "); dbprint(i_print,"// equisingularity ideal is computed!"); return(list(ES_Id,ESfix_Id)); } } example { "EXAMPLE:"; echo=2; ring r=0,(x,y),ds; poly f=x7+y7+(x-y)^2*x2y2; list K=esIdeal(f); option(redSB); // Wall's equisingularity ideal: std(K[1]); ring rr=0,(x,y),ds; poly f=x4+4x3y+6x2y2+4xy3+y4+2x2y15+4xy16+2y17+xy23+y24+y30+y31; list K=esIdeal(f); vdim(std(K[1])); // the latter should be equal to: tau_es(f); } /////////////////////////////////////////////////////////////////////////////// proc esStratum (poly p_F, list #) "USAGE: esStratum(F[,m,L]); F poly, m int, L list ASSUME: F defines a deformation of a reduced bivariate polynomial f and the characteristic of the basering does not divide mult(f). @* If nv is the number of variables of the basering, then the first nv-2 variables are the deformation parameters. @* If the basering is a qring, ideal(basering) must only depend on the deformation parameters. COMPUTE: equations for the stratum of equisingular deformations with fixed (trivial) section. RETURN: list l: either consisting of an ideal and an integer, where @format l[1]=ideal defining the equisingular stratum l[2]=1 if some error has occured, l[2]=0 otherwise; @end format or consisting of a ring and an integer, where @format l[1]=ESSring is a ring extension of basering containing the ideal ES (describing the ES-stratum) and the poly p_F=F, l[2]=1 if some error has occured, l[2]=0 otherwise. @end format NOTE: L is supposed to be the output of reddevelop (with the given ordering of the variables appearing in f). @* If m is given, the ES Stratum over A/maxideal(m) is computed. @* This procedure uses @code{execute} or calls a procedure using @code{execute}. printlevel>=2 displays additional information. SEE ALSO: esIdeal, isEquising KEYWORDS: equisingular stratum EXAMPLE: example esStratum; shows examples. " { list l=esComputation (0,p_F,#); return(l); } example { "EXAMPLE:"; echo=2; int p=printlevel; printlevel=1; ring r = 0,(a,b,c,d,e,f,g,x,y),ds; poly F = (x2+2xy+y2+x5)+ax+by+cx2+dxy+ey2+fx3+gx4; list M = esStratum(F); M[1]; printlevel=3; // displays additional information esStratum(F,2); // es stratum over Q[a,b,c,d,e,f,g] / ^2 ideal I = f-fa,e+b; qring q = std(I); poly F = imap(r,F); esStratum(F); printlevel=p; } /////////////////////////////////////////////////////////////////////////////// proc isEquising (poly p_F, list #) "USAGE: isEquising(F[,m,L]); F poly, m int, L list ASSUME: F defines a deformation of a reduced bivariate polynomial f and the characteristic of the basering does not divide mult(f). @* If nv is the number of variables of the basering, then the first nv-2 variables are the deformation parameters. @* If the basering is a qring, ideal(basering) must only depend on the deformation parameters. COMPUTE: tests if the given family is equisingular along the trivial section. RETURN: int: 1 if the family is equisingular, 0 otherwise. NOTE: L is supposed to be the output of reddevelop (with the given ordering of the variables appearing in f). @* If m is given, the family is considered over A/maxideal(m). @* This procedure uses @code{execute} or calls a procedure using @code{execute}. printlevel>=2 displays additional information. EXAMPLE: example isEquising; shows examples. " { int check=esComputation (1,p_F,#); return(check); } example { "EXAMPLE:"; echo=2; ring r = 0,(a,b,x,y),ds; poly F = (x2+2xy+y2+x5)+ay3+bx5; isEquising(F); ideal I = ideal(a); qring q = std(I); poly F = imap(r,F); isEquising(F); ring rr=0,(A,B,C,x,y),ls; poly f=x7+y7+(x-y)^2*x2y2; poly F=f+A*y*diff(f,x)+B*x*diff(f,x); isEquising(F); isEquising(F,2); // computation over Q[a,b] / ^2 } /* Examples: LIB "ES.lib"; ring r = 0,(x,y),ds; poly p1 = y^2+x^3; poly p2 = p1^2+x5y; poly p3 = p2^2+x^10*p1; poly p=p3^2+x^20*p2; p; versal(p); setring Px; poly F=Fs[1,1]; int t=timer; list M=esStratum(F); timer-t; //-> 3 LIB "ES.lib"; option(prot); printlevel=2; ring r=0,(x,y),ds; poly f=(x-yx+y2)^2-(y+x)^31; versal(f); setring Px; poly F=Fs[1,1]; int t=timer; list M=esStratum(F); timer-t; //-> 233 LIB "ES.lib"; printlevel=2; option(prot); timer=1; ring r=32003,(x,y),ds; poly f=(x4-y4)^2-x10; versal(f); setring Px; poly F=Fs[1,1]; int t=timer; list M=esStratum(F,3); timer-t; //-> 8 LIB "ES.lib"; printlevel=2; timer=1; ring rr=0,(x,y),ls; poly f=x7+y7+(x-y)^2*x2y2; list K=esIdeal(f); // tau_es vdim(std(K[1])); //-> 22 // tau_es_fix vdim(std(K[2])); //-> 24 LIB "ES.lib"; printlevel=2; timer=1; ring rr=0,(x,y),ls; poly f=x7+y7+(x-y)^2*x2y2+x2y4; // Newton non-deg. list K=esIdeal(f); // tau_es vdim(std(K[1])); //-> 21 // tau_es_fix vdim(std(K[2])); //-> 23 LIB "ES.lib"; ring r=0,(w,v),ds; poly f=w2-v199; list L=reddevelop(f); // change basering to HNEring poly f=fetch(r,f); versal(f); setring Px; list L=imap(HNEring,L); poly F=Fs[1,1]; list M=esStratum(F,2,L); LIB "ES.lib"; printlevel=2; timer=1; ring rr=0,(A,B,C,x,y),ls; poly f=x7+y7+(x-y)^2*x2y2; poly F=f+A*y*diff(f,x)+B*x*diff(f,x)+C*diff(f,y); list M=esStratum(F,6); std(M[1]); // standard basis of equisingularity ideal LIB "ES.lib"; printlevel=2; timer=1; ring rr=0,(x,y),ls; poly f=x20+y7+(x-y)^2*x2y2+x2y4; // Newton non-degenerated list K=esIdeal(f); ring rr=0,(x,y),ls; poly f=x6y-3x4y4-x4y5+3x2y7-x4y6+2x2y8-y10+2x2y9-y11+x2y10-y12-y13; list K=esIdeal(f); versal(f); setring Px; poly F=Fs[1,1]; list M=esStratum(F,2); LIB "ES.lib"; printlevel=2; ring rr=0,(x,y),ls; poly f=x6y-3x4y4-x4y5+3x2y7-x4y6+2x2y8-y10+2x2y9-y11+x2y10-y12-y13; list K=esIdeal(f); vdim(std(K[1])); //-> 51 tau_es(f); //-> 51 printlevel=3; f=f*(y-x2)*(y2-x3)*(x-y5); int t=timer; list L=esIdeal(f); vdim(std(L[1])); //-> 99 timer-t; //-> 42 t=timer; tau_es(f); //-> 99 timer-t; //-> 23 LIB "ES.lib"; printlevel=3; ring rr=0,(x,y),ds; poly f=x4+4x3y+6x2y2+4xy3+y4+2x2y15+4xy16+2y17+xy23+y24+y30+y31; list K=esIdeal(f); vdim(std(K[1])); //-> 68 tau_es(f); //-> 68 versal(f); setring Px; poly F=Fs[1,1]; int t=timer; list M=esStratum(F); timer-t; //-> 0 */