LIB "matrix.lib";
LIB "sing.lib";
LIB "primdec.lib";




/////////////Compute Codimension and Milnor number /////////

proc coDimiL(ideal f)
{
  module D=jacob(f);
  module TK=maxideal(1)*D+f*freemodule(2);   //you forgot maxideal(1)
  module STK=std(TK);
  int c=vdim(STK)-2;                                          // -2 since it is dim( maxideal(1)*freemodule(2)/TK)
  int mi=vdim(std(f))-1;
  list L=c,mi;
  return(L);
 }
/////////////Compute Sigma //////////////////////////////////

proc siGmaa(ideal f)
{
  int k=1;
  module D=jacob(f);
  module TK=maxideal(1)*D+f*freemodule(2);
  module STK=std(TK);
  int t = size(reduce(maxideal(k)*freemodule(2),STK));
  while(t!=0)
    {
      k=k+1;
      t= size(reduce(maxideal(k)*freemodule(2),STK));
      
       } 
     return(k);
    }








/////////////Compute Fold number //////////////////////////////////

proc dFoldn(ideal I)
 {
  int k, h;
  k=1;
 ideal I1=I+maxideal(k);
 ideal I2=I+maxideal(k+1);
 h=vdim(std(I2))-vdim(std(I1));
 while (h==2)
     {  k=k+1;
        I1=I+maxideal(k);
        I2=I+maxideal(k+1); 
        h=vdim(std(I2))-vdim(std(I1));       
        }  
        return(k-1);
        }








/////////////Compute Boardman Symbol //////////////////////////
proc Br1Symbol(ideal I)
    {
      int k;
      list L;
      matrix M=jacob(I);
     ideal T=minor(M,2);
      ideal J=std(I+T);
      ideal J1=std(reduce(I,std(J)));
      poly f=reduce(1,std(J));
      if(f!=0 && J1==0)
       {
        ideal K= minor(M,1);
        ideal K1=std(I+K);
        ideal K2=std(reduce(J,std(K1)));
        poly g=reduce(1,std(K1));
        if(g!=0 && K2==0)
         {
           return(list(2,K1));
           }
         }
       else
         { return(list(0,0));}
           } 
         
     proc Br2Symbol(ideal I)
      {
      intvec u,v;
      list L;
      L=Br1Symbol(I);
      if(L[1]!=0)
       {
        matrix N=jacob(L[2]);
        ideal P=minor(N,2);
        ideal A=std(L[2]+P);
        ideal A1=std(reduce(L[2],std(A)));
        poly h=reduce(1,std(A));
        if(h!=0 && A1==0)
         {
           ideal B=std(minor(N,1));
           ideal B1=std(L[2]+B);
           ideal B2=std(reduce(A,std(B1)));
           poly h1=reduce(1,std(B1));
           if(h1==0 && B2==0)
            {
           return(list(2,1));
           }
       else
         { return(list(0,0));}
           } 
         }
       }

/////////////Decide m>l or m<=l //////////////////////////
proc valueK(ideal f)
{
  matrix M=jacob(f);
  poly q=det(M);
  ideal J=radical(q);
  poly g=J[1];
  int k=milnor(g);
    return(k);
      }

///////////// Compute Alpha //////////////////////////
proc vAlpha(int a,int b)
{
  if(a<=b)
   {
    return(a);
     }
  if(b<=a)
   {
    return(b);
     }
      }

///////////// reduce to x2 //////////////////////////


proc red2Jet(ideal I)
 {
   ideal J;
   poly f=I[1];
   poly g=jet(f,2);
   list L=factorize(g);
   poly f1=L[1][2];
   poly g1=subst(f1,y,0);
   number b1=leadcoef(g1);
   poly g2=subst(f1,x,0);
   number b2=leadcoef(g2);
   if(b1!=0)
    {
       map phi=basering,(x-g2)/b1,y;
       J=phi(I);
      }
  else
   {
      map shi=basering,y,x/b2;
      J=shi(I);
      }
     return(J);   
     }



proc red2jet2(ideal I,int k)
 {
   poly f=I[1];
   poly h=f-x2;
   map phi;
   while(ord(h/x)<k && h!=0)
    {
      phi=basering,x-1/2*h/x,y;
      I=phi(I);
      return(I);
       }  
   }

///////////// Value of lambda for F and tilde G ///////////////////

proc valuelambda(ideal I,int k,int l)
{
   poly f;
   ideal J;
   J=red2jet2(I,k);
   int c=k+l;
   int d=l+k div 2;
        f=jet(J[2],d);
        module f1=[f];
        module f2=[J[1]];
        module f3=jet(system("reduce_bound",f1,f2,c),c);
        matrix f4=f3;
        f=jet(f4[1,1],d); 
        f=f-x*(f/x);
    if(ord(f)>d)
     { return(0);}
    else
    { 
       return(leadcoef(f));
    }
}












///////////// Main Procedure //////////////////////////



proc contactMapgerms(ideal I)
{
  int c,mu,d,sig,de,a,b,e;
  number lmd;
  ideal J1;
  ideal J2;
  a=ord(I[1]);
  b=ord(I[2]);
  e=vAlpha(a,b);
  list B= Br2Symbol(I);
  list L=coDimiL(I);
  c=L[1];
  mu=L[2];
  d= dFoldn(I);
  de=valueK(I);
  sig= siGmaa(I);
    if(B[1]==2&&B[2]==1)
     {
       if(de==1)
        {
////////////////////C-Type /////////////////////////////////

          if(c==3*d+1)
           {
       return("f is of type C_"+string(d+1)+"");
             }
////////////////////E-Type /////////////////////////////////

          if(c<3*d+1)
           {
     return("f is of type E_"+string(c-2*d)+","+string(d+1)+"");
             }
            }
        if(de>1)
         {
////////////////////D-Type /////////////////////////////////
           if(mu==2*sig-1 && mu==c-d)
            {
     return("f is of type D_"+string(d)+","+string(sig)+"");
             }
////////////////////G_k,l,m-Type /////////////////////////////////
           if(mu==2*sig-1 && mu!=c-d)
            {
return("f is of type G_"+string(c-d-sig+1)+","+string(d)+" ,"+string(sig)+"");
             }
////////////////////F-Type /////////////////////////////////
           if(mu!=2*sig-1 && c-mu+d!=sig)
            {
              if(mu==sig+2*d-1)
              {
     return("f is of type F_"+string(sig)+","+string(d)+"");
                }
////////////////////G^i_k/2,l,n-Type /////////////////////////////////
           if(mu<sig+2*d-1)           
            {
              int k= c-d-sig+1;
               if(k mod 2==0)
                {
return("f is of type G^i_"+string((c-d-sig+1)*1/2)+","+string(d)+" ,"+string(sig)+"");
                }
               }
              }
////////////////////F_k,l or tildG_k/2,l -Type ////////////////
           if(c-mu+d==sig)           
            {
              int p=mu-2*d+1;
              int alf=vAlpha(p-1,d);
              if(c==mu+alf && p mod 2!=0)
               {
     return("f is of type F_"+string(p)+","+string(d)+"");
                }
               if(c==mu+alf && p mod 2==0)
               {
                 J1=red2Jet(I);
                 J2=red2jet2(J1,p);
                 lmd=valuelambda(J2,p,d);
                  if(lmd==0)
                   {
           return("f is of type F_"+string(p)+","+string(d)+"");
                }
                  if(lmd!=0)
                   {
      return("f is of type tildG _"+string(p div 2)+","+string(d)+"");

                   }
                }
               
////////////////////G_k,l,m or tildG_k/2,l,n -Type ////////////////

              if(c!=mu+alf && p mod 2!=0)
               {
return("f is of type G_"+string(p)+","+string(d)+" ,"+string(sig)+"");

                }
               if(c!=mu+alf && p mod 2==0)
               {
                 J1=red2Jet(I);
                 J2=red2jet2(J1,p);
                 lmd=valuelambda(J2,p,d);
                  if(lmd==0)
                   {
return("f is of type G_"+string(p)+","+string(d)+" ,"+string(sig)+"");

                }
                  if(lmd!=0)
                   {
return("f is of type tildG_"+string(p div 2)+","+string(d)+" ,"+string(sig)+"");

                   }
                }
               }


              }
             }
            }