//////// Procedure to classify the right unimodal and bimodal Hypersurface singularities in positive characteristic by invariants///////////////////
////////// Authors: Amir Shehzad, Muhammad Ahsan Binyamin, Junaid Alam Khan, Umer Shuaib, Khawar Mehmood/////////////
LIB"general.lib";
LIB "hnoether.lib";
proc uniBimodalReq(poly f)
{ /// main loop-start
int p=char(basering); // charterstic of the basering
int m=mult(std(f)); // multiplicity
int u=vdim(std(jacob(f))); // milnornumber
int T=vdim(std(jacob(f)+f));// Tjurina number
int d1=1;
int d2=2;
int d3=3;
int mu1=1;
int mu2=2;
int mu3=3;
int mu4=4;
////////////////////// Non-isolated Case//////////////////////
if(m==-1)
{
return("f is a non-isolated singularity");
}
////////////////////// Smooth Case////////////////////////////
if(m==0)
{
return("f is smooth");
}
if(m==1)
{
//================= the A_1 ===================================
return("A1=xy");
}
///////////////////////// multiplicity 2/// A-type/////////////////////
if(m==2)
{ //// multiplicity 2 loop start
//////////right simple////////////////
if (u>1 && u<=p-2)
{
classifyReq(f);
}
/////////////unimodal of type A ////////////////
if (u>=p && u<=2*p-2)
{
poly g=red2Jet2(f);
list L=red2NF1(g);
return("f is unimodal of type A"+string(u)+": x2+"+string(leadcoef(L[2]))+"y "+string(p)+" + y"+string(u+1)+" ");
}
//////////////bimodal of type A ////////////////
if (u>=2*p && u<=3*p-2)
{
poly g=red2Jet2(f);
list L=red2NF1(g);
list Q =red2NF11(L[1],p);
return("f is bimodal of type A"+string(u)+": x2+"+string(leadcoef(L[2]))+"y "+string(p)+" + y "+string(2*p)+"+ y"+string(u+1)+" ");
}
return(0);
} //// multiplicity 2 loop end
///////////////////////// multiplicity 3///////D,J,E type/////////////
if(m==3)
{ //// multiplicity 3 loop start
poly g=jet(f,3);
///////////// jet(f,3)=x2y//// D type////////////////
list F12=sortFacMult(g,d1,mu2);
if (size(F12)==1)
////right simple/////////////(jet(f,3)=x2y)///
{
if(u>=4 && u
p+1 && u<=2*p-1)
{
poly g=req3Jet2(f);
list L=red2NF(g);
number d=leadcoef(L[2]); // alpha;
return("f is unimodal of type D"+string(u)+":x2y+"+string(d)+"y"+string(p)+"+y"+string(u-1)+"");
}
///////////bimodal//////D-type///////(jet(f,3)=x2y)///
if (u==2*p)
{
poly g=req3Jet2(f);
list L=red2NF(g);
number d=leadcoef(L[2]); // alpha;
return("f is bimodal of type D"+string(u)+":x2y+"+string(d)+"y"+string(p)+"+y"+string(2*p-1)+"");
}
if (u>2*p+1 && u<=3*p-1)
{
poly g=req3Jet2(f);
list L=red2NF(g);
number d=leadcoef(L[2]); // alpha
list P=red2NF2(g,p);
number e=leadcoef(P[2]); // beta
return("f is bimodal of type D"+string(u)+":x2y+"+string(d)+"y"+string(p)+"+"+string(e)+"y"+string(2*p)+"+y"+string(u-1)+"");
}
return(0);
}
////////////////////jet(f,3)=x3 ///////// Case of J & E-type///////////////////
list F13=sortFacMult(g,d1,mu3);
if (size(F13)==1)
{
int u2=Blow2up(f); // 2nd blow+milnor number
int u3=Blow3up(f); // 3rd blow+milnor number
poly f1=stricTrans(f);
poly h=jet(f1,3);
////right simple////////////////
if(u>=6 && u<=8)
{
classifyReq(f);
}
///////////unimodal///J & E/////////(jet(f,3)=x3)///////
if (u==10)
{
return("f is unimodal of type J2,0");
}
if (u==11)
{
return("f is unimodal of type J2,1");
}
if (u>=12 && u<=14)
{
if (u2==u-11)
{
return("f is unimodal of type J2,"+string(u-10)+"");
}
if (u2==u-12)
{
return("f is unimodal of type E"+string(u)+"");
}
}
if (u==15 && p>11)
{
return("f is unimodal of type J2,5");
}
//////// unimodal & bimodal/// J & E type// milnor number 16 OR 17///////
list BF111 = sortFacMult(h,d1,mu1); // degree-multiplicity-number of factors required
list BF121 = sortFacMult(h,d1,mu2); // degree-multiplicity-number of factors required
list BF131 = sortFacMult(h,d1,mu3); // degree-multiplicity-number of factors required
if (u==16 || u==17) ///////// milnor number 16 0R 17////J type
{
if(u
p+4 && u<2*p+4 && size(BF111)==1 && size(BF121)==1)
{
return("f is bimodal of type J2,"+string(u-10)+"");
}
if(u
p+4 && u<2*p+4 && size(BF111)==1 && size(BF121)==1)
{
return("f is bimodal of type J2,"+string(u-10)+"");
}
if(u
20) ///////// milnor number u>20////J type
{
if(u
p+4 && u<2*p+4 && size(BF111)==1 && size(BF121)==1)
{
return("f is bimodal of type J2,"+string(u-10)+"");
}
if(u
9 && u
=17 && u<2*p+6 && size(B2H121)==1 )
{
return("f is bimodal of type W#1,"+string(u-15)+"");
}
if(u==18 && h2==0 && size(B3H131)==1)
{
return("f is bimodal of type W18");
}
if (u>=17 && u
17)
{
return("f is bimodal of type W1,"+string(u-15)+"");
}
}
return(0);
}
///////////////////////jet(f,4)=one linear factor of multiplicity 3 and one linear factor of multiplicity 1 =x3y //////////////////////////////
list F131=sortFacMult(g,d1,mu3);
if (size(F131)==1) // degree-multiplicity-number of factors required
{
int u2z=Blow2up1(f);
if (u>=11 && u<=13)
{
return("f is unimodal of type Z"+string(u)+"");
}
if (u>=17 && u<=19 && u2z==u-17)
{
return("f is bimodal of type Z"+string(u)+"");
}
if (u==15)
{
return("f is bimodal of type Z1,0");
}
if (u>=16 && u<=p+8 && u2z==u-16)
{
return("f is bimodal of type Z1,"+string(u-15)+"");
}
return(0);
}
}//// multiplicity 4 loop end
} /// main loop end
//////////////========== LIB"classifyRightsimple.lib"======= Procedure to classify the right simple singularities in positive characteristic===////////
///// Authors: Deeba Afzal, Muhammad Ahsan Binyamin and Faira Kanwal Janjua//////
//////==This procedure based on the algorithm given in the article “On the Classification of Simple Singularities in Positive Characteristic ”, /////
///////==An. St. Univ. Ovidius Constanta, Vol- 22(2), 2014, 5-19, Article Authors: Deeba Afzal, Muhammad Ahsan Binyamin and Faira Kanwal Janjua===///////////
LIB "classify.lib";
LIB "presolve.lib";
proc classifyReq(poly f)
{
def R=basering;
int n=nvars(R);
ring S=char(basering),(x(1..n)),(c,ds);
poly f=fetch(R,f);
int p=char(basering);
ideal J=std(jacob(f));
int k=system("HC")+1;
int m=vdim(J);
int c=corank(f);
if(m2)
//================= corank > 2 not simple =====================
{
setring R;
return("f is not simple");
}
if(c==1)
{
setring R;
if(m<=p-2)
{
//================= the A_m singularities ==================
return("A"+string(m)+"=x"+string(m+1)+"+y2");
}
return("f is not simple");
}
//==================== the case of corank 2 ===================
//==================== the splitting lemma ===================
poly g =MorseR(jet(f,k),k,c)[1];
ideal I=findvars(g);
int i;
intvec v;
for(i=1;i<=n;i++)
{
if(var(i)==I[1])
{
v[1]=i;
}
if(var(i)==I[2])
{
v[2]=i;
}
}
ring T=p,(x,y),ds;
ideal J;
J[v[1]]=x;
J[v[2]]=y;
map phi=S,J;
poly f=phi(g);
poly g=jet(f,3);
if(g==0)
{
setring R;
return("f is not simple");
}
list L=factorize(g);
int b=size(L[1]);
if(b==2)
{
if(L[2][2]==3)
{
//================= the E_m singularities ==================
if(m >= 9)
{
setring R;
return("f is not simple");
}
setring R;
if(m==6 && p!=3)
{
return("E6=x3+y4");
}
if(m==7 && p!=3)
{
return("E7=x3+xy3");
}
if(m==8 && p!=3 && p!=5)
{
return("E8=x3+y5");
}
}
else
{
//===================== the D_4 ===============================
// In this case j^3(f) is irreducible over the base ring and therefore
// splits into 3 factors over the algebraic closure of the field
setring R;
return("D4=x2y+y3");
}
}
setring R;
//================= the D_m singularities ======================
if(b>=3 && 4<=m && m<=p)
{
return("D"+string(m)+"=x2y+y"+string(m-1)+"");
}
//================= the remaining cases =======================
return("f is not simple");
}
////////////////////////////////// procedures from classify.lib ///////////////
proc MorseR(poly f, int K, int corank)
{
poly fc, f2, a, P, Beta, fi;
ideal Jfx, B;
int n, i, j, k, Rang, Done;
matrix Mat;
map Id, Psi, Phi, PhiG;
intvec Abb, RFlg;
list v;
fi = f;
n = nvars(basering);
def ring_top=basering;
for( j=1; j<=n ; j++)
{
Abb[j] = 0;
}
RFlg = GetRfR(fi, n);
PhiG=ring_top,maxideal(1);
//----------------- find quadratic term, if there is only one -------
B = maxideal(1);
if(corank == (n-1))
{
Done = 0;
f2 = jet(fi, 2);
j = 1;
Jfx = f2, diff(f2, x(j));
while(j<=n && (diff(f2, x(j))==0))
{
j++;
Jfx = f2, diff(f2, x(j));
}
Mat = matrix(syz(Jfx));
Beta = 2*Mat[2,1]/Mat[1,1];
for( j=1; j<=n ; j++)
{
f2 = CoeffR(Beta, x(RFlg[j]), x(RFlg[j]));
if(f2!=0)
{
k = RFlg[j];
break;
}
}
for( j=1; j<=n ; j=j+1)
{
f2 = CoeffR(Beta, x(j), x(j));
if(j == k)
{
B[rvar(x(j))] = (2*f2*x(j)-Beta) / number(f2);
}
}
Phi =ring_top,B;
fi = Phi(fi);
PhiG = Phi(PhiG);
}
//------------------------ compute spliting lemma -------------------
fc = fi;
i = 1;
while( i <= n)
{
Phi=ring_top,maxideal(1);
j = i + 1;
f2 = jet(fc,2);
if( (f2 - subst(f2, x(RFlg[i]), 0)) == 0 )
{
Abb[RFlg[i]] = 1;
}
if( (f2 - subst(f2, x(RFlg[i]), 0)) != 0 )
{
while( (j<=n) || (i==n) )
{
a = CoeffR(jet(fc,2), x(RFlg[i]), x(RFlg[i])^2);
if( (a != 0) || (i==n) )
{
break;
}
B = maxideal(1);
for( k=1; k<=n ; k=k+1)
{
if(k==RFlg[j])
{
B[rvar(x(k))] = x(k) + x(RFlg[i]);
}
}
Phi = ring_top,B;
fc = Phi(fi);
j++;
}
PhiG = Phi(PhiG);
if( (j<=n) || (i==n))
{
P = CoeffR(fc, x(RFlg[i]), x(RFlg[i]));
if(P != 0) { P = P / number (2*a);
B = maxideal(1);
for(k=1;k<=n;k=k+1)
{
if(k==RFlg[i])
{
B[rvar(x(k))]=x(k)-P;
}
}
Phi =ring_top,B;
fi = Phi(fc);
PhiG = Phi(PhiG);
fc = jet(fi,K);
P = CoeffR(fc, x(RFlg[i]), x(RFlg[i]));
if( P != 0)
{
fi = fc;
continue;
}
}
}
}
fi = fc; i++;
}
//--------------------------- collect results --------------
Rang = 0;
B = maxideal(1);
for(i=1;i<=n;i++)
{
if(Abb[i]!=1)
{
Rang ++;
B[rvar(x(i))]=0;}} Phi = ring_top,B;
PhiG = Phi(PhiG);
fi = Phi(fi);
v = fi, PhiG;
return(v);
}
proc GetRfR(poly fi, int n)
{
//---------------------------- initialisation --------------
int j, k, l1, l1w;
matrix Koef;
intvec RFlg;
RFlg[n] = 0;
intvec Ha = RFlg;
for( j=1; j<=n ; j=j+1)
{
Koef = coef(fi, x(j));
Ha[j] = ncols(Koef);
if(CoeffR(fi, x(j),0) == 0)
{
Ha[j] = Ha[j] + 1;
}
}
for( j=n; j>0 ; j=j-1)
{
l1 = 0;
l1w = 0;
for(k=1;k<=n;k=k+1)
{
if(Ha[k]>l1w)
{
l1=k;l1w=Ha[k];
}
}
RFlg[j] = l1;
Ha[l1] = 0;
}
return(RFlg);
}
//////////////////////////////////////////////////////////////
proc CoeffR(poly f, list #)
{
poly a, term;
int n, i;
matrix K;
n = nvars(basering);
i = 1;
term = #[2];
K = coef(f, #[1]);
while((i<=ncols(K)) && (K[1,i] != term))
{
i++;
if(i>ncols(K))
{
break;
}
}
if(i<=ncols(K))
{
a = K[2,i];
}
if(i>ncols(K))
{
a = 0;
}
return(a);
}
/////////////////////////////// Support Programs//////////////////////////////////////////
proc sortFacMult(poly f,int d,int mu)/// factors of given degree and multiplicity
{ list L=factorize(f,2);
int n=size(L[1]);
int i,j;
poly g;
list Q;
for(i=1;i<=n;i++)
{ g=L[1][i];
if(deg(g)==d)
{ if (L[2][i]==mu)
{ j=j+1;
Q[j]=g;
}
}
}
return(Q);
}
///////////////////////////////////////
proc red3Jet(poly f)
{
def R=basering;
poly g1=jet(f,3);
list L1=factorize(g1);
poly h=(1/L1[1][1])*f;
poly g=jet(h,3);
list L=factorize(g);
poly h1=subst(L[1][2],y,0);
number a=leadcoef(h1);
poly h2=subst(L[1][2],x,0);
number b=leadcoef(h2);
if(a!=0)
{
map phi=R,1/a*x-b/a*y,y;
poly f1=phi(h);
}
else
{
map shi=R,y,1/b*x;
poly f1=shi(h);
}
return(f1);
}
/////////// milnor number///2nd blow up////////multiplicity 3////J & E type///
proc Blow2up(poly f)
{
int bl_2;
poly f1=red3Jet(f);
poly g=subst(f1,x,xy);
poly g1=g/y3;
poly h=subst(g1,x,xy);
poly h1=h/y3;
bl_2=vdim(std(jacob(h1)));
return(bl_2);
}
////////// milnor number///// 3rd blow up////////multiplicity 3////J & E type///
proc Blow3up(poly f)
{
int bl_3;
poly f1=red3Jet(f);
poly g=subst(f1,x,xy);
poly g1=g/y3;
poly h=subst(g1,x,xy);
poly h1=h/y3;
poly t=subst(h1,x,xy);
poly t1=t/y3;
bl_3=vdim(std(jacob(t1)));
return(bl_3);
}
//////////////////////////////////
proc red3Jet1(poly f)
{
poly g=jet(f,3);
list L=factorize(g);
poly f1=L[1][2];
poly f2=L[1][3];
poly g1=subst(f1,y,0);
number b1=leadcoef(g1);
poly g2=subst(f1,x,0);
number b2=leadcoef(g2);
poly h1=subst(f2,y,0);
number c1=leadcoef(h1);
poly h2=subst(f2,x,0);
number c2=leadcoef(h2);
if(b1!=0)
{
poly h=subst(f,x,(x-g2)/b1);
}
else
{
poly h=subst(f,x,(x-h2)/c1);
}
return(h);
}
////////////////////////////////////////////////////
proc req3Jet1(poly f)
{
def R=basering;
poly T;
map phi=R,y,x;
poly h=red3Jet1(f);
poly h1=jet(h,3);
list M=factorize(h1);
if(subst(M[1][2],x,0)==0)
{
poly f1=M[1][3];
}
else
{
poly f1=M[1][2];
}
poly g1=subst(f1,y,0);
number b1=leadcoef(g1);
poly g2=subst(f1,x,0);
number b2=leadcoef(g2);
poly g=subst(h,y,(y-g1)/b2);
poly g3=jet(g,3);
poly g4=g3/y;
poly g5=subst(g4,y,0);
if(g5!=0)
{
return(g);
}
else
{
T=phi(g);
}
return(T);
}
//////////////////////////////////
proc red4Jet1(poly f)
{
poly g=jet(f,4);
list L=factorize(g);
poly f1=L[1][2];
poly f2=L[1][3];
poly g1=subst(f1,y,0);
number b1=leadcoef(g1);
poly g2=subst(f1,x,0);
number b2=leadcoef(g2);
poly h1=subst(f2,y,0);
number c1=leadcoef(h1);
poly h2=subst(f2,x,0);
number c2=leadcoef(h2);
if(b1!=0)
{
poly h=subst(f,x,(x-g2)/b1);
}
else
{
poly h=subst(f,x,(x-h2)/c1);
}
return(h);
}
//////////////////////////////////////////////////
proc req4Jet1(poly f)
{
def R=basering;
poly T;
map phi=R,y,x;
poly h=red4Jet1(f);
poly h1=jet(h,4);
list M=factorize(h1);
if(subst(M[1][2],x,0)==0)
{
poly f1=M[1][3];
}
else
{
poly f1=M[1][2];
}
poly g1=subst(f1,y,0);
number b1=leadcoef(g1);
poly g2=subst(f1,x,0);
number b2=leadcoef(g2);
poly g=subst(h,y,(y-g1)/b2);
poly g3=jet(g,4);
poly g4=g3/y;
poly g5=subst(g4,y,0);
if(g5!=0)
{
return(g);
}
else
{
T=phi(g);
}
return(T);
}
//////////////////////////milnor number/// 2nd blowup////multiplicity 4////Z-type////////
proc Blow2up1(poly f)
{
int a;
int b;
poly f1=req4Jet1(f);
poly g=subst(f1,x,xy);
poly g1=g/y4;
list L=factorize(jet(g1,3));
if(size(L[1])==2)
{
poly f2=red3Jet(g1);
poly h=subst(f2,x,xy);
poly h1=h/y3;
a=vdim(std(jacob(h1)));
return(a);
}
if(size(L[1])==3)
{
poly f3=req3Jet1(g1);
poly h2=subst(f3,x,xy);
poly h3=h2/y3;
b=vdim(std(jacob(h3)));
return(b);
}
}
///////////////// 1st Strict Transform/////multiplicity 3///J type/////
proc stricTrans(poly f)
{
int a;
poly f1=red3Jet(f);
poly g=subst(f1,x,xy);
poly g1=g/y3;
return(g1);
}
///////////////////////////////////////////////////////
proc sortsmall(int i, int j)
{ list K;
if(i<=j)
{
K[1]=i;
K[2]=j;
}
if(j 3/////////alpha coefficient
proc red2NF(poly f)
{
if(f==x2*y){return(0);}
int k=4;
poly l,m,p;
while(1)
{
while(jet(f-x*y2,k)==0)
{
k++;
}
poly f1=jet(f,k)-x^2*y;
poly f2=(f1-subst(f1,x,0)- subst(f1,y,0))/xy;
poly f3= subst(f2,x,0);
poly g1=(f1- subst(f1,x,0))/x;
poly g2=(g1- subst(g1,x,0))/x;
m=subst(f,x,x-1/2*f3);
p=subst(m,y,y-g2);
if(leadcoef(jet(p,k)-jet(p,k-1))!=0)
{
return(leadcoef(jet(p,k)-jet(p,k-1)));
}
k=k+1;
}
}
//////////////////////////////
/////////////Compute modulus values/////////beta coefficient
proc red2NF(poly f)
{
if(f==x2*y){return(0);}
int k=4;
poly l,m,p;
while(1)
{
while(jet(f-x2*y,k)==0)
{
k++;
}
poly f1=jet(f,k)-x^2*y;
poly f2=(f1-subst(f1,x,0)- subst(f1,y,0))/xy;
poly f3= subst(f2,x,0);
poly g1=(f1- subst(f1,x,0))/x;
poly g2=(g1- subst(g1,x,0))/x;
f=subst(f,x,x-1/2*f3);
f=subst(f,y,y-g2);
if(leadcoef(jet(f,k)-jet(f,k-1))!=0)
{
return(list(f,leadcoef(jet(f,k)-jet(f,k-1))));
}
k=k+1;
}
}
proc red2NF2(poly f,int t)
{
if(f== jet(f,t)){return(0);}
int k=t+1;
poly l,m,p;
while(1)
{
while(jet(f- jet(f,t),k)==0)
{
k++;
}
poly f1=jet(f,k)- jet(f,t);
poly f2=(f1-subst(f1,x,0)- subst(f1,y,0))/xy;
poly f3= subst(f2,x,0);
poly g1=(f1- subst(f1,x,0))/x;
poly g2=(g1- subst(g1,x,0))/x;
f=subst(f,x,x-1/2*f3);
f=subst(f,y,y-g2);
if(leadcoef(jet(f,k)-jet(f,k-1))!=0)
{
return(list(jet(f,k),leadcoef(jet(f,k)-jet(f,k-1))));
}
k=k+1;
}
}
/////////////////////////////////////////////////
/////////////////////////Reduce Fourth jet int0 x2y2///////////////////
proc red4Jet(poly f)
{
poly g=jet(f,4);
list L=factorize(g);
poly f1=L[1][2];
poly f2=L[1][3];
poly g1=subst(f1,y,0);
number b1=leadcoef(g1);
poly g2=subst(f1,x,0);
number b2=leadcoef(g2);
poly g3=subst(f2,y,0);
number b3=leadcoef(g3);
poly g4=subst(f2,x,0);
number b4=leadcoef(g4);
if(b1!=0)
{
poly h=subst(f,x,x-(b2/b1)*y);
}
else
{
poly h=subst(f,x,x-(b4/b3)*y);
}
return(h);
}
proc req4Jet(poly f)
{
poly h=red4Jet(f);
poly h1=jet(h,4);
list M=factorize(h1);
if(subst(M[1][2],x,0)==0)
{
poly f1=M[1][3];
}
else
{
poly f1=M[1][2];
}
poly g1=subst(f1,y,0);
number b1=leadcoef(g1);
poly g2=subst(f1,x,0);
number b2=leadcoef(g2);
poly g=subst(h,y,y-(b1/b2)*x);
return(g);
}
/////////////////////////////