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
*/