LIB "classifyMapGerms.lib";


///////////////////////////////////////////////////////////////////////////////
proc classifyPlaneMaps1(ideal I)
"USAGE:   classifyPlaneMaps1(I);  I=an ideal with 2 generators in a polynomial
         ring with 2 variables and local ordering defining a map germ C^2 to C^2
COMPUTE:  The normal form of the germ in Riegers classification if it is
          simple
RETURN:   normal form of I, of type ideal
NOTE:     If I has codimension greater than 6 it returns (0,0)
EXAMPLE:  example classifyPlaneMaps1; shows an example"

{
   if(size(jet(I,1))==0){return(ideal(0,0));}
   list L=coDimMap(I,1);
   int determinacy=L[1];
   int c=L[2];
   poly g=I[2];
   if(determinacy<0){return(ideal(0,0));}
   I=normalMap(I,determinacy);
   int m =vdim(std(I));
   int cu=vdim(std(diff(I[2],y),diff(diff(I[2],y),y)));

   if((cu==-1)&&(m>=2)){return(ideal(0,0));}
   if(m==1){return(maxideal(1));}

   if(m==2){return(ideal(var(1),var(2)^2));}
   
   if(m==7){return(ideal(var(1),var(1)*var(2)+var(2)^7+var(2)^9+var(2)^10+var(2)^11));}
      
      if(m==3)
      {
        if(cu==1)
        {
           return(ideal(var(1),var(1)*var(2)+var(2)^3));
        }
        if(cu==2)
        {
           return(ideal(var(1),var(2)^3+var(1)^2*var(2)));
        }
       if(cu==3)
        {
           return(ideal(var(1),var(2)^3+var(1)^3*var(2)));
         }
       if(cu==4)
        {
           return(ideal(var(1),var(2)^3+var(1)^4*var(2)));
          }
      if(cu==5)
        {
           return(ideal(var(1),var(2)^3+var(1)^5*var(2)));
        }
      
      
      
      
      return(ideal(0,0));
   }
   if(m==4)
   {
      if(cu==2){return(ideal(var(1),var(1)*var(2)+var(2)^4));}
      if(cu==6){return(ideal(var(1),var(1)^3*var(2)+var(1)^2*var(2)^2+var(2)^4+var(1)^3*var(2)^2));}
      if(cu==3)
      {
        if(c==4)
        
          {
         return(ideal(var(1),var(1)*var(2)^2+var(2)^4+var(2)^5));
          }
        
        if(c==5)
        {

         return(ideal(var(1),var(1)*var(2)^2+var(2)^4+var(2)^7));
        }
        if(c==6)
        {

         return(ideal(var(1),var(1)*var(2)^2+var(2)^4+var(2)^9));
        }
        }

      if(cu==4)
      {
        if(c==5)
        {

         return(ideal(var(1),var(1)^2*var(2)+var(2)^4+var(2)^5));
        }
        if(c==6)
        {

         return(ideal(var(1),var(1)^2*var(2)+var(2)^4));
        }
        }
        
        return(ideal(0,0));
       }
      if(m==5)
      {
       if((cu==3) && (c==4))
         {

         return(ideal(var(1),var(1)*var(2)+var(2)^5+var(2)^7));
         }
        if((cu==3) && (c==5))
        {

         return(ideal(var(1),var(1)*var(2)+var(2)^5));
        }
        
        if(cu==4)
      {
        if(c==5)
        {

         return(ideal(var(1),var(1)*var(2)^2+var(2)^5+var(2)^6));
        }
        if(c==6)
        {

         return(ideal(var(1),var(1)*var(2)^2+var(2)^5+var(2)^9));
        }
        }
        
        if(cu==6){return(ideal(var(1),var(1)^2*var(2)+var(1)*var(2)^3+var(2)^5+var(2)^6+var(2)^7));}
        
       
       }
      if(m==6)
      {
       if(cu==4)
      {
        if(c==6)
         {
            g=jet(g,9);
            I[2]=modulus(g);
            return(I);
         }

        }

       if(cu==5){return(ideal(var(1),var(1)*var(2)^2+var(2)^6+var(2)^7+var(2)^9));}
         }
        return(ideal(0,0));
       }

///////////////////////////////////////////////////////////////////////////////
 proc normalMap(ideal I, int bound)
"USAGE:    normalMap(I);  I=f1,f2 in K[x,y]
 RETURN:   (0,0) if jet(I,1)=0
           (x,h) A-equivalent to I modulo <x,y>^bound+1
           if ord(h)=1 then h=y
           if ord(h)=2 then h=xy+y^m(I)+h.o.t.if m(I)>=3 or (x,y^2)
           if ord(h)=3 milnor(diff(h,y))=1 and m(I)>3 then h=xy2+y^m(I)+h.o.t.
           if ord(h)=4 and m(I)=4 then h=y^4+a*x^3y+b*x^2y^2+h.o.t.
           m(I)=vdim(std(I))
 NOTE: the procedure is designed for simple and unimodular maps
"
{
   I=jet(I,bound);
   //the trivial cases
   if(size(jet(I,1))==0){return(ideal(0,0));}
   if(vdim(std(jet(I,1)))==1){return(maxideal(1));}
   if(vdim(std(I))==2){I=var(1),var(2)^2;return(I);}

   def R=basering;
   map phi;
   poly p;
   int i;

   //transformation to (x,h)
   if(jet(I[1],1)==0){p=I[1];I[1]=I[2];I[2]=p;}
   I[1]=simplify(I[1],1);
   if(lead(I[1])==var(2)){phi=R,var(2),var(1);}
   else{phi=R,var(1)-I[1][2],var(2);}
   I=phi(I);
   while(I[1]!=var(1))
   {
      phi=R,2*var(1)-I[1],var(2);
      I=jet(phi(I),bound);
   }
   while((leadmonom(I[2])==var(1))||(leadmonom(I[2])==var(1)^2)||(leadmonom(I[2]
   )==var(1)^3)||(leadmonom(I[2])==var(1)^4))
   {
      I[2]=I[2]-lead(I[2]);
   }
   I[2]=simplify(I[2],1);
   if(deg(lead(I[2]))>=5){return(I);}

   //special treatment of h
   //the case ord(h)=4
   if(deg(lead(I[2]))==4)
   {
      if(vdim(std(I))>4){return(I);}
      def S=changevar("var(2),var(1)",R);
      setring S;
      map psi=R,var(2),var(1);
      ideal I=psi(I);
      I[2]=simplify(I[2],1);
      if(leadmonom(I[2][2])==var(1)^3*var(2))
      {
         map sigma=S,var(1)-1/4*leadcoef(I[2][2])*var(2),var(2);
         I=sigma(I);
         setring R;
         map lambda=S,var(2),var(1);
         I=lambda(I);
         if(leadmonom(I[2])==var(1)^4){I[2]=I[2]-lead(I[2]);}
         return(I);
      }
      else
      {
         setring R;
         return(I);
      }
   }

   //the case ord(h)=3
   if(deg(lead(I[2]))==3)
   {
      phi=R,var(1),var(2)-leadcoef(I[2])/(2*leadcoef(I[2][2]))*var(1);
      I=phi(I);
   }
   while((leadmonom(I[2])==var(1))||(leadmonom(I[2])==var(1)^2)||(leadmonom(I[2]
   )==var(1)^3)||(leadmonom(I[2])==var(1)^4))
   {
      I[2]=I[2]-lead(I[2]);
   }
   I[2]=simplify(I[2],1);
   if(leadmonom(I[2])==var(1)*var(2)^2)
   {
      while(i<bound)
      {
         i++;
         p=jet(I[2]/(var(1)*var(2))-var(2),i);
         phi=R,var(1),var(2)-1/2*p;
         I=jet(phi(I),bound);
      }
   }

   //the case ord(h)=2
   if(leadmonom(I[2])==var(1)*var(2))
   {
      while(i<bound)
      {
         i++;
         p=jet(I[2]/var(1)-var(2),i);
         phi=R,var(1),var(2)-p;
         I=jet(phi(I),bound);
      }
   }

   //all cases
   p=lead(I[2])+reduc(I[2]-lead(I[2]),I,bound);
   I[2]=p;
   if((leadmonom(I[2])==var(1)*var(2))||(leadmonom(I[2])==var(1)*var(2)^2))
   {
      number a=leadcoef(I[2][2]);
      I=subst(I,x,a*x);
      I=1/a*I;
   }
   return(I);
}

///////////////////////////////////////////////////////////////////////////////
  proc reduc(poly f,ideal I, int b)
"USAGE:   reduc(f, I, b)
COMPUTE:  the normal form of f with respect to the algebra generated by I
"
{
   list L=1;
   map psi;
   f=jet(f,b);
   if(f==0){return(f);}
   while((f!=0) && (L[1]!=0))
   {
     L= algebra_containment(lead(f),lead(I),1);
     if (L[1]==1)
     {
        def S= L[2];
        psi= S,maxideal(1),I;
        f=jet(f-psi(check),b);
        kill S;
     }
   }
   return (lead(f)+reduc(f-lead(f),I,b));
}

///////////////////////////////////////////////////////////////////////////////
 proc modulus(poly g)
// g=xy+y6+h.o.t. or xy2+y6+h.o.t.
{
   number a,b,c;
   matrix m=coef(g,var(2));
   if(m[1,1]==var(2)^9)
   {
      c=leadcoef(m[2,1]);
      if(m[1,2]==var(2)^8)
      {
         b=leadcoef(m[2,2]);
         if(m[1,3]==var(2)^7)
         {
            a=leadcoef(m[2,3]);
         }
      }
      else
      {
         if(m[1,2]==var(2)^7)
         {
           a=leadcoef(m[2,2]);
         }
      }
   }
   else
   {
      if(m[1,1]==var(2)^8)
      {
         b=leadcoef(m[2,1]);
         if(m[1,2]==var(2)^7)
         {
            a=leadcoef(m[2,2]);
         }
      }
      else
      {
         if(m[1,1]==var(2)^7)
         {
           a=leadcoef(m[2,1]);
         }
      }
   }
   if(leadmonom(g)==var(1)*var(2))
   {
      if(a!=0)
      {
         g=var(1)*var(2)+var(2)^6+(5*b-3*a^2)/(5*a^2)*var(2)^8+(25*c+14*a^3-35*a*b)/(25*a^3)*var(2)^9;
      }
   }
   if(leadmonom(g)==var(1)*var(2)^2)
   {
      if(a!=0)
      {
         g=var(1)*var(2)^2+var(2)^6+var(2)^7+(4*c*a-5*b)/(4*a^2)*var(2)^9;
      }
   }
   return(g);
}