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);
        }


//============================Boardman Symbol==================================
    

proc deltaI(ideal I, int s)

{
   
  int n=nvars(basering);
 
  if(n-s+1>0)
   {
  ideal J=minor(jacob(I),n-s+1);
 
  return(I+J);
   } 
  else
   {
     return(0);
    }
  }




   //==============================================================

proc extensionI(ideal I)

{
   
  int i;
   
  ideal J,JS,K;
   
  K=I;
   
  while(1)
   
  {
      
    i++;
      
    J=deltaI(I,i);
      
       JS=std(J);
      
       if(JS[1]==1||JS[1]==0)
       {
      return(list(K,i-1));
      }
      
      
    K=J;
   
   }

}

   //==============================================================

proc Brsymbol(ideal I)

{
   
  list L;
   
  ideal IS=std(I);
   
  intvec v;
   
  int q;
   
   
  while(1)
   
  {
     
    q++;
     
    L=extensionI(I);
     
    v[q]=L[2];
     
    I=L[1];
     
    if(size(reduce(I,IS))==0)
     
    {
        
      break;
     
       }
    
    IS=std(I);
   
   }
   
    return(v);

}



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



proc KSimPlePlaneGerms(ideal I)
{
  int c,mu,d,s,de,a,b,e;
  
  
  list L=coDimiL(I);
  c=L[1];
  mu=L[2];
   d= dFoldn(I);
  s= siGmaa(I);
    intvec v=Brsymbol(I);
   if(c==mu && v[1]==1 )
        {
        return("f is of type A_"+string(c)+"");
        }
       if(c==mu+1 && v[1]==2 && v[2]==0)
        {
        return("f is of type B_"+string(c- s)+","+string(s)+"");
        }
        
        
if(v[1]==2&&v[2]==1)
   {
////////////////////C-Type /////////////////////////////////      
          if(c==3*s-2 && mu==2*s-1 && s< 5)
           {
       return("f is of type C_"+string(s)+"");
             }

          if(c==mu+2 && c==2*s+1 && d>2)
           {
     return("f is of type E_"+string(3)+","+string(s)+"");
             }
          if(c==mu+2 && c==2*s)
           {
     return("f is of type F_"+string(3)+","+string(s-2)+"");
             }
         if(c==mu+2 && c==s+5 && d==2)
           {
     return("f is of type F_"+string(s)+","+string(2)+"");
             }
            }
else{return(0);}
            }