LIB"sagbi.lib";

///////////// Procedure based on the main theorem/////

proc comp( ideal G, ideal C, intmat m)

{ int n=nvars(basering);
  def bsr=basering;
 ideal I=G;
 ideal J=C;
 ideal L,S;
 L=lead(J);
 int i,j,p,k;
 p=IP(L);

 if (p==0)

    {    
 
      ERROR("composition is not iterated powering");
    } 

ideal q=std(L);

  if (size(q)< n)

   
    {    
 
      ERROR("composition is not iterated powering");
    } 

intmat W[n][n]= mak(J);
intmat O[n][n]= m*W;

execute("ring R1=("+charstr(bsr)+"),("+varstr(bsr)+"),(M(O)) ;") ;

ideal H=imap(bsr,I);

 H=sagbi(H,0);

setring bsr;

map phi= R1,C;


ideal GC=phi(H);

 
return(GC);

}



////////////////////////////////////

///// Procedure for the checking of Iterative Power//////


proc IP( ideal I)

{ int n=nvars(basering);
  
 ideal m=maxideal(1);
 ideal J=lead(I);
  int i,j,p,k;
  p=1;
  poly f;
  list K;
 for(i=1;i<=size(I);i++)
    {  
       f=J[i];
    
      for(k=1;k<=n;k++)

      { 
     
       if ( leadmonom(f)/leadmonom(m[k])!=0)
          
       {
             for(j=k+1;j<=n;j++)
            
             {
               
                if(leadmonom(f)/leadmonom(m[j])!=0) 
                {
                  p=0;
                  break;    
                }
                
              }
          
        }
      }

    }
return(p);

}

/////////////////////////////////////////////

//////////// Input is the "composition ideal", Output is the matrix (useful in the computation)


proc mak( ideal I)

{ int n=nvars(basering);
  
 ideal J=I;
  int i,j,p,k;
 intvec v,w,u;
 p=n*n;
 intmat m[n][n];
 w=leadexp(J[1]);
 for(i=2;i<=n;i++)
    {  
       v=leadexp(J[i]);

       w=w,v;
       v=0;  
       
    }

m=w[1..p];
m=transpose(m);
return(m);

}


///////////////////////////////////////////////////