LIB"classifyMapGerms.lib";
LIB"sing.lib";

////////////////////////////////////////////////////////////////////////////////
proc class(ideal I)
{
   //=========     classification of simple (K^2,0)--->(K^3,0), ===========
   //=========             I ideal defining this map            ===========
   //=========          assume that R has 2 variables           ===========
   //=========  return the normal form or -1 if not simple      ===========

   def R=basering;
   if(nvars(R)>2){return(-1);}
   if(ncols(I)>3){return(-1);}

   int i,mu;
   poly p;
   map phi;

   //===============  compute the multiplicity mu =========================
   //========= and keep only the cases with multiplicity 1,2,3  ===========
   mu=vdim(std(I));   //multiplicity
   if((mu>3)||(mu<0)){return(-1);}  //not simple
   if(mu==1){return(ideal(var(1),var(2)));}
   if(size(I)<3){return(-1);}       //not simple

   //======= finds an element of order 1 in I and change it to I[1] =======
   for(i=1;i<=ncols(I);i++)
   {
      if(ord(I[i])==1){break;}
   }
   if(i>ncols(I)){return(-1);} //not rank 1
   p=I[i];
   I[i]=I[1];
   I[1]=p;     
   
   //=====            computes a bound L[1] of the determinacy        =====
   //=====                 and the A^e-codimension  L[2]              =====
   list L=coDimMap(I); 

   //== maps I to <x,b,c> and kills the pure powers of x in I[2] and I[3]== 
   I=norMap(I,L[1]);   
   I[2]=y*(I[2]/y);    
   I[3]=y*(I[3]/y);

   //= computes the coefficient matrix of the quadratic part of I[2],I[3] =
   ideal J=I[2],I[3];
   J=jet(J,2);         
   if(size(J)==0){return(-1);} 
   matrix M=coef(J,var(2)); 
   M=subst(M,var(1),1); 
 
   //= check if the determinant is different from zero (the case x,y^2,xy)= 
   if((ncols(M)==2)&&(nrows(M)==3))
   {
      M=submat(M,2..3,1..2);
      poly d=det(M);
      if(d!=0){return(ideal(var(1),var(2)^2,var(1)*var(2)));}
   }
   
   //==== Now the determinant is zero and we change I to jet(I[2],2)!=0 ===
   //====              and jet(I[3],2)=0                                ===
   if(J[1]==0){p=I[2];I[2]=I[3];I[3]=p;}
   I[2]=simplify(I[2],1);
   I[3]=simplify(I[3],1);
   I[3]=I[3]-I[2];
   p=jet(I[2],2);

   //=========== the case  x,xy+y^(3k-1), k+2 the codimension    ==========
   if(p==var(1)*var(2))   //the case <x,xy+hot,c>
   {
      I[2]=var(1)*var(2)+var(2)^(3*L[2]-1);
      I[3]=var(2)^3;
      return(I);
   }

   //========= jet(I[2],2)=xy+cy^2 and c!=0                      ==========
   //========= make a transformation xy+cy^2 ----> y^2           ==========
   if(leadmonom(p)==var(1)*var(2)) //the case <x,xy+s*y^2+hot,c>
   {
      number b=leadcoef(I[2][2]);
      phi=R,var(1),var(2)-1/(2*b)*var(1);
      I=phi(I);
      I[2]=y*(I[2]/y);   //kill pure powers of x
      I[3]=y*(I[3]/y);
      I=jet(I,L[1]);
      I[2]=simplify(I[2],1);
   }

   //========  now jet(I[2],2)=y^2+hot and transform it to y^2   ==========
   while(I[2]!=var(2)^2)     //the case <x,y^2+hot,c>
   {
      p=(I[2]-var(2)^2)/y;
      phi=R,var(1),var(2)-1/2*p;
      I=phi(I);
      I[2]=y*(I[2]/y);   //kill pure powers of x
      I[3]=y*(I[3]/y);
      I=jet(I,L[1]);
   }

   //========  now I[2]=y^2 and we look for the normal form of I[3] ========
   I[3]=clean(I[3]);     //brings I[3] to the form yp(x,y^2)
   I[3]=simplify(I[3],1);
   I[3]=findType(I[3],L[2]+2); //finds the type of I[3]
   if(I[3]==-1){return(-1);}  //not simple
   return(I);
}

////////////////////////////////////////////////////////////////////////////////
proc norMap(ideal I, int bound)
{
// returns (x,b,c) A-equivalent to I modulo <x,y>^bound+1
// assumes that ord I[1]=1

   I=jet(I,bound);
   def R=basering;   //assume that the ordering is ds or Ds
   map phi;
   I[1]=simplify(I[1],1);
   if(lead(I[1])!=var(1))
   {
      phi=R,var(2),var(1);
      I=phi(I);
   }
   phi=R,2*var(1)-I[1],var(2);
   I=phi(I);
   I=jet(I,bound);
   while(size(I[1])>1)
   {
      phi=R,2*var(1)-I[1],var(2);
      I=phi(I);
      I=jet(I,bound);
   }
   return(I);
}
////////////////////////////////////////////////////////////////////////////////
proc clean(poly p)
{
   // removes the monomials  x^iy^2j from p
   intvec v;
   while(1)
   {
      if(p==0){return(p);}
      v=leadexp(p);
      if((v[2] mod 2)==0)
      {
         p=p-lead(p);
      }
      else
      {
         break
      }
   }
   return(lead(p)+clean(p-lead(p)));
}
////////////////////////////////////////////////////////////////////////////////
proc findType(poly p,int k)
{
//k is the codimension
   p=p/var(2);
   if(ord(p)==2)
   {
      if(milnor(p)==k-2){return(var(2)^3+var(1)^(k-1)*var(2));}
      return(var(1)^2*var(2)+var(2)^(2*k-3));
   }
   if(ord(p)==3)
   {
      if(leadmonom(p)==var(1)^3){return(var(1)^3*var(2)+var(2)^5);}
      return(var(1)*var(2)^3+var(1)^(k-2)*var(2));
   }
   return(-1); 
}
////////////////////////////////////////////////////////////////////////////////
proc mix(ideal I,ideal J)
{
//to create examples
   if(size(I)<=1){return(I);}
   def R=basering;
   map phi=R,J;
   I=phi(I);
   poly p;
   p=I[1];
   I[1]=I[2];
   I[2]=p;
   I[2]=I[2]+I[1];
   return(I);
}
///////////////////////////      examples      /////////////////////////////////
ring R=0,(x,y),ds;
ideal J=x+2y+7xy,3x+y+xy2;

ideal I=x,y;
class(I);
I=mix(I,J);
class(I);

ideal I=x,y2,xy;
class(I);
I=mix(I,J);
class(I);

ideal I=x,y2,y3+x7y;
class(I);
I=mix(I,J);
class(I);

ideal I=x,y2,x2y+y5;
class(I);
I=mix(I,J);
class(I);

ideal I=x,y2,xy3+x4y;
class(I);
I=mix(I,J);
class(I);

ideal I=x,y2,x3y+y5;
class(I);
I=mix(I,J);
class(I);

ideal I=x,xy+y8,y3;
class(I);
I=mix(I,J);
class(I);

ideal I=y;
class(I);
I=mix(I,J);
class(I);

ideal I=x,xy+y8,y5;
class(I);
I=mix(I,J);
class(I);

ideal I=x,y2,y5;
class(I);
I=mix(I,J);
class(I);