########################################################################################
#            Computation of Involutive Basis for Syzygy Module                         #
########################################################################################

##############################################
#            Involutive Divisions            #
##############################################
##############Janet Division#############
Janet:=proc(u,U,Vars)
local n,m,d,L,i,j,dd,V,v,Mult;
#option trace;
if degree(u)=0 then RETURN(Vars); fi;
n:=nops(Vars);
m:=ArrayNumElems(U);
d:=[seq(max(seq(degree(U[j],Vars[i]),j=1..m)),i=1..n)];
Mult:=NULL;
    if degree(u,Vars[1])=d[1] then
       Mult:=Mult,Vars[1];
    fi:
for j from 2 to n do
    	dd:=[seq(degree(u,Vars[l]),l=1..j-1)];
    	V:=NULL:
    	for v in U do
    	    if [seq(degree(v,Vars[l]),l=1..j-1)]=dd then
	       V:=V,v:
	    fi:
	od:
	if degree(u,Vars[j])=max(seq(degree(v,Vars[j]), v in [V])) then
	   Mult:=Mult,Vars[j];
	fi:
od:
RETURN([Mult]);
end:

########### Pommaret Division #############
LeftPommaret:=proc(u,U,Vars)
local N,Ind,i;
N:=NULL:
Ind:=indets(u):
for i from 1 to nops(Vars) while not (Vars[i] in Ind) do
    N:=N,Vars[i]:
od:
N:=N,Vars[i]:
RETURN([N]);
end:
###############################
RightPommaret:=proc(u,U,Vars)
local N,Ind,i;
N:=NULL:
Ind:=indets(u):
for i from  nops(Vars) by -1 to 1 while not (Vars[i] in Ind) do
    N:=N,Vars[i]:
od:
N:=N,Vars[i]:
RETURN([N]);
end:
###################################
#     Schreyer’s ordering test    #
###################################
with(Groebner):
Sher1:=proc(p,q)
#option trace;
local i,j,a,b;
global tord,LLMI;
 i := p[2];
 j := q[2];  
if i < j then 
     return(true);
   else 
    return(false); 
   fi;
end:
########################
with(Groebner):
Sher22:=proc(p,q)
#option trace;
local i,j,a,b;
global tord,LLMI;
 i := p[6][2];
 j := q[6][2];  
 a:=p[7][2];
 b:=q[7][2];
if a<>b then
  RETURN(TestOrder(a,b,tord));
  elif i < j then 
     return(true);
   elif i >= j then
    return(false); 
   fi;
end:
################TOP###################
ord:=proc(p,q)
#option trace;
local i,j,a,b;
global tord,LR,Vars;
if p[1]<>q[1] then
RETURN(TestOrder(p[1],q[1],tord));
else
    RETURN(evalb(p[2]<q[2]));
   fi;
end:
############################
#         Reduction        # 
############################                          
with(PolynomialIdeals):
############################
LM:=proc(f)
global tord;
if f<>0 then
  RETURN(LeadingMonomial(f,tord));
fi:
RETURN(0);
end:

#########################
 LT := proc (f)
 local A,B;
 global tord;
 if f=0 then return(0); fi;
 A,B:=LeadingTerm(p,trd);
  return(A*B);
end proc:
############################
LC:=proc(f)
global tord;
  RETURN(LeadingCoefficient(f,tord));
end: 
#############################
Crit1:=proc(p,gg,pp)
global c1;
local InvDiv,flag,g,l;
#option trace;
l:=p[2][2];
g:=gg[2][2];
    if l*g=p[7][2] then 
    c1:=c1+1; 
    RETURN(true);
fi;
RETURN(false);
end:
################ Normal Form Algorithm #################
with(Groebner):
Normalform:=proc(pp,TermOrder,InvDivision)
#option trace;
local LMF,Vars,p,r,flag,g,u,Mult,InvDiv,FF,j,i,qq,c1,m,M,A,B,QQQ,a,b,cd,LO;
global tord,T,Polys,ISR2,sy,C3,rdz;
LMF:=copy(Array([seq(v[7][2], v in T)])); 
Vars:=[op(TermOrder)];
p:=Polys[pp[1]]:                                 
r:=0;
FF:=copy(Array([seq(Polys[v[1]], v in T)]));
j:=0:m:=ArrayNumElems(FF);
qq:=NULL;
QQQ:=NULL;
A,B:=op(pp[7]);
while p<>0 do     
      j:=j+1:
      flag:=false:
      i:=1;
      while not flag and i<=m do
          g:=FF[i]:
      	  if divide(B,LMF[i],'u') then
	           Mult:=InvDivision(LMF[i],LMF,Vars);
	           if indets(u) subset {op(Mult)} then 
	           if j=1 then
	              cd:=[pp[6],[u*T[i][6][1],T[i][6][2]]];
	              M := sort(cd, proc (a, b) options operator, arrow; Sher1(a, b) end proc);
	              if M[-1]=pp[6] then
	               qq:=pp[6];  
	             fi;
	          fi; 
	     	   flag:=true:                   
               if j=1 then
                 if Crit1(pp,T[i],InvDivision) and pp[8][1]<>true then                 
                sy(ArrayNumElems(sy)+1):=[[Polys[T[i][2][1]]*Polys[pp[2][1]],Polys[pp[2][1]],pp[2][3]],[Polys[T[i][2][1]]*Polys[pp[2][1]],Polys[T[i][2][1]],T[i][2][3]]];
                   RETURN([0,[0,qq],[[Polys[pp[2][1]],T[i][2][1],i],[Polys[T[i][2][1]],pp[2][1],pp[2][1]]]]);  
                 fi:
               fi:
               QQQ:=QQQ,[A/T[i][7][1]*u,T[i][1],i];
		       p:=simplify(p-(A/T[i][7][1])*u*g);
		       A,B:=LeadingTerm(p,TermOrder);
	        else 
                 i:=i+1;
	        fi:
	     else 
            i:=i+1;
	  fi:
    od:
      if not flag  then 
            j:=j+1;
	        r:=simplify(r+A*B);
            p:=simplify(p-A*B);
            A,B:=LeadingTerm(p,TermOrder);
      fi:
od:
if r=0 then
sy(ArrayNumElems(sy)+1):=[seq([Polys[a[3]]*a[1],a[1],a[3]] , a in [QQQ]),[-Polys[pp[1]],pp[4],pp[8][3]]];
if pp[8][1]=true then      
   for t from 1 to ArrayNumElems(T) do
      ZS:=NULL;           
      for i from 1 to nops(T[t][8][2])-1 do 
          if pp[1]=T[t][8][2][i][2] then                                                                    
              ZS:=ZS,seq([tt[1]*T[t][8][2][i][1],tt[2],tt[3]], tt in [QQQ]);                                         
           else      
          ZS:=ZS,T[t][8][2][i];                        
          fi; 
        od;
       T[t][8][2]:= [ZS,T[t][8][2][-1]];
    od;  
 fi;
fi;

if r<>0 then
if pp[8][1]=true then      
   for t from 1 to ArrayNumElems(T) do
      ZS:=NULL;           
      for i from 1 to nops(T[t][8][2])-1 do 
          if pp[1]=T[t][8][2][i][2] then                     
              ZS:=ZS,seq([tt[1]*T[t][8][2][i][1],tt[2],tt[3]], tt in [QQQ]),[T[t][8][2][i][1],p[1],ArrayNumElems(T)+1];            
           else      
          ZS:=ZS,T[t][8][2][i];                        
          fi; 
        od;
       T[t][8][2]:= [ZS,T[t][8][2][-1]];
    od;  
 fi;
 fi;        
if r=0  then 
  rdz:=rdz+1: 
fi:        
RETURN([r,[[LeadingTerm(r,tord)],qq],[QQQ]]);
end:
#####################################################
#            Involutive Basis                       #
#####################################################
InvBasis:=proc(Ideal,Tord,InvDivision)
#option trace;
global tord,T,Q,Polys,Vars,c1,c2,rdz,LMF,critelim, reduct,isr,deg,LLMI, mmtelim,sy;
local firsttime, firstbytes, cov, F, n, syz,  p, FLAGE, j, hh, h, k, rpp, t, i, SS, S, s, TT, jj, fl, c, q, NM, y, Q2, Q3, ii0;
F:=Array(1..nops(Ideal),Ideal);
F:=copy(sort(F,(a,b)->TestOrder(a,b,Tord)));
n:=ArrayNumElems(F):
syz:=Array([]);sy:=Array([]);
LLMI:=Array(1..n,[seq([LeadingTerm(i,tord)], i in F)]);
deg:=0;
Polys := copy(Array(1 .. n, proc (i) options operator, arrow; F[i] end proc));
T:=copy(Array(1..1,[[1,[1,LLMI[1][2],1],{},1,1,[1,1],LLMI[1],[false,[[1,1,1],[-Polys[1],1,1]]]]]));
Q:=copy(Array(1..ArrayNumElems(F)-1,[seq([i,[i,LLMI[i][2],i],{},1,i,[1,i],LLMI[i],[false,1,i]],i=2..n)]));
Q := sort(Q, proc (a, b) options operator, arrow; Sher22(a, b) end proc);
while ArrayNumElems(Q) <> 0 do 
  p := Q[1];                                           
  deg := max(deg, degree(Polys[p[1]])); 
  Q := copy(Array(1 .. ArrayNumElems(Q)-1, proc(t) options operator, arrow; Q[t+1] end proc)); 
         FLAGE := false;             
       if FLAGE = false then
           hh:=Normalform(p,tord,InvDivision);
           h:=hh[1]: 
           k := ArrayNumElems(T)+1;
          rpp:= Polys[p[1]];
 
          if nops(hh[2])>1 then                     
          if hh[2][2][1]<>1 then            
            syz(ArrayNumElems(syz)+1):=hh[2][2]; 
         fi;  
         fi;                
        if h=0 and p[1]=p[2][1] and p[8][1]<>true  then
	    SS:=Array([]);
            S:=copy(Q):
             for s in S do 
	           if s[2][1]<>p[1] then	                   
	       	     SS(ArrayNumElems(SS)+1):=s;
		    else
		     mmtelim :=mmtelim+1:                                                              
                  fi;
              od:
             Q:=copy(SS);
         elif h<>0 and hh[2][1][2]<>p[7][2] then 
              k := ArrayNumElems(T)+1;
	          Polys[p[1]]:=h; 
	          T(k):=[p[1],[p[1],hh[2][1][2],k],{},1,p[1],[1,p[1]],hh[2][1],[false, [[1,p[1],k],seq(t, t in  hh[3]),[-rpp,p[4],p[8][3]]] ]];
              TT:=Array([]);            
                for k from 1 to ArrayNumElems(T) do
	               if divide(T[k][7][2], hh[2][1][2]) and hh[2][1][2] <> T[k][7][2] then
	      	           Q(ArrayNumElems(Q)+1):=[seq(T[k][i], i=1..7),[true,T[i][8][2][-1][2],T[i][8][2][-1][3]]];
	      	           for jj from 1 to ArrayNumElems(T) do
	      	                   T[jj][2][3]:=T[jj][2][3]-1;
                               fl:=[seq( T[jj][8][2][g][3] , g=1..nops( T[jj][8][2]))];
	      	             for c from 1 to nops(fl) do
	      	    	       if  fl[c] >= ArrayNumElems(TT)+1  then	      	          
	      	                 T[jj][8][2][c][3]:=T[jj][8][2][c][3]-1;	      	          
	      	               fi;	      	   
	      	             od;	      	            
	      	         od;	      	           
                       else 
		           TT(ArrayNumElems(TT)+1):=T[k]; 
                       fi:                
                od:               
                 T:=copy(TT);
                   for i from 1 to ArrayNumElems(T) do
       	               q:=T[i];
       	               NM:={op(Vars)} minus {op(InvDivision(q[7][2],Array([seq(w[7][2], w in T)]),Vars))};
	                       for y in NM minus q[3] do
		                   n:=n+1: 
	                           Polys(n):=expand(y*Polys[q[1]]);	                       
	                           Q(ArrayNumElems(Q)+1):=[n,q[2],{},y,q[1],[expand(y*q[6][1]),q[6][2]],[q[7][1],y*q[7][2]],[false,y,q[8][2][1][3]]];	                            
	   	                   od:
                            T[i][3]:=q[3] union NM;
                      od:
           elif h<>0 and hh[2][1][2]=p[7][2] then
	           Polys[p[1]]:=h;
	          
	           T(k):=[seq(p[i],i=1..7),[false, [[1,p[1],k],seq(t , t in hh[3]),[-rpp,p[4],p[8][3]]]]];
	             
                    for i from 1 to ArrayNumElems(T) do
       	                 q:=T[i];
                 	 NM:={op(Vars)} minus {op(InvDivision(q[7][2],Array([seq(w[7][2], w in T)]),Vars))};
	                 for y in NM minus q[3] do 
		             n:=n+1: 
	                    Polys(n):=expand(y*Polys[q[1]]);
	                    Q(ArrayNumElems(Q)+1):=[n,q[2],{},y,q[1],[expand(y*q[6][1]),q[6][2]],[q[7][1],y*q[7][2]],[false,y,q[8][2][1][3]]];                      
	   	         od:
                         T[i][3]:=q[3] union NM;
                       od; 
              fi;  
         Q:=convert(Q,list); 
         Q := sort(Q, proc (a, b) options operator, arrow; Sher22(a, b) end proc);      
         Q:=copy(Array(1..nops(Q),Q));              
   fi;   
od:
SY:=Array(1..ArrayNumElems(T), [seq(tt[8][2] , tt in T)] );
SY1:=Array([]);
for i in SY do 
   if i[-1][2]<>1 then
      SY1(ArrayNumElems(SY1)+1):= [seq([Polys[i[ii][2]]*i[ii][1],i[ii][1],i[ii][3]],ii=1..nops(i)-1), i[-1]];                                                           
   fi;
 od;
SYY:=Array(1..ArrayNumElems(sy)+ArrayNumElems(SY1) ,[seq(a , a in SY1),seq( b ,b in sy)]);                                                                                                     
TH:=Array(1..ArrayNumElems(T),[seq([1,tt[7],1], tt in T)]):
THH:=Array(1..ArrayNumElems(T),[seq(tt[7][2], tt in T)]):
T:=Array(1..ArrayNumElems(T),[seq([Polys[T[tt][1]],[1,tt]], tt=1..ArrayNumElems(T))]);
RETURN([T,TH,THH,SYY]);
end:
######################################################################################
######################################################################################





########################################################################################
#                        Pommaret Basis Computation                                    #
########################################################################################


##############################################
#            Involutive Divisions            #
##############################################

##############Janet Division#############
Janet:=proc(u,U,Vars)
local n,m,d,L,i,j,dd,V,v,Mult;
#option trace;
if degree(u)=0 then RETURN(Vars); fi;
n:=nops(Vars);
m:=ArrayNumElems(U);
d:=[seq(max(seq(degree(U[j],Vars[i]),j=1..m)),i=1..n)];
Mult:=NULL;
    if degree(u,Vars[1])=d[1] then
       Mult:=Mult,Vars[1];
    fi:
for j from 2 to n do
    	dd:=[seq(degree(u,Vars[l]),l=1..j-1)];
    	V:=NULL:
    	for v in U do
    	    if [seq(degree(v,Vars[l]),l=1..j-1)]=dd then
	       V:=V,v:
	    fi:
	od:
	if degree(u,Vars[j])=max(seq(degree(v,Vars[j]), v in [V])) then
	   Mult:=Mult,Vars[j];
	fi:
od:
RETURN([Mult]);
end:

########### Pommaret Division #############
LeftPommaret:=proc(u,U,Vars)
local N,Ind,i;
N:=NULL:
Ind:=indets(u):
for i from 1 to nops(Vars) while not (Vars[i] in Ind) do
    N:=N,Vars[i]:
od:
N:=N,Vars[i]:
RETURN([N]);
end:
    	   

RightPommaret:=proc(u,U,Vars)
local N,Ind,i;
N:=NULL:
Ind:=indets(u):
for i from  nops(Vars) by -1 to 1 while not (Vars[i] in Ind) do
    N:=N,Vars[i]:
od:
N:=N,Vars[i]:
RETURN([N]);
end:
###############################
#    Leading Monomial         #
###############################                            
with(PolynomialIdeals):
#################
with(Groebner):
LM:=proc(f)
global tord;
if f<>0 then
  RETURN(LeadingMonomial(f,tord));
fi:
RETURN(0);
end:

##################
 LT := proc (f)
 local A,B;
 global tord;
 if f=0 then return(0); fi;
 A,B:=LeadingTerm(p,trd);
  return(A*B);
end proc:
#################
LC:=proc(f)
global tord;
    RETURN(LeadingCoefficient(f,tord));
end: 
#################################
#     Schreyer’s ordering test  #
#################################
with(Groebner):
Sher1:=proc(p,q)
#option trace;
local i,j,a,b;
global tord,LLMI;
 i := p[2];
 j := q[2];  
if i < j then 
     return(true);
   else 
    return(false); 
   fi;
end:

########################
with(Groebner):
Sher22:=proc(p,q)
#option trace;
local i,j,a,b;
global tord,LLMI;
 i := p[6][2];
 j := q[6][2];  
 a:=p[7][2];
 b:=q[7][2];
if a<>b then
  RETURN(TestOrder(a,b,tord));
  elif i < j then 
     return(true);
   elif i >= j then
    return(false); 
   fi;
end:ord:=proc(p,q)
#option trace;##TOP3
local i,j,a,b;
global tord,LR,Vars;
if p[1]<>q[1] then
RETURN(TestOrder(p[1],q[1],tord));
else
    RETURN(evalb(p[3]<q[3]));
   fi;
end:
#########################
 LTSY := proc (f)
 #option trace;
 local A,B;
 global tord;
 FF:=[seq([expand(f[i][1]),f[i][2],f[i][3]], i=1..nops(f))];
 A:= sort(FF, proc (a, b) options operator, arrow; ord(a, b) end proc)[-1];
  return([LM(A[2]),A[3]]);
end proc:


############################
#         Reduction        # 
############################
Crit1:=proc(p,gg,pp)
global c1;
local InvDiv,flag,g,l;
#option trace;
l:=p[2][2];
g:=gg[2][2];
    if l*g=p[7][2] then 
    c1:=c1+1; 
    RETURN(true);
fi;
RETURN(false);
end:
############################
Crit2:=proc(p,gg,pp)
global c2;
local g,l,lc;
#option trace;
l:=p[2][2];
g:=gg[2][2];
lc:=lcm(l,g);
if  divide(p[7][2],lc) and p[7][2]<>lc  then
    c2:=c2+1; 
    RETURN(true);
fi;
RETURN(false);
end:

###############################
#   Normal Form Algorithm     #
###############################
with(Groebner):
Normalform:=proc(pp,TermOrder,InvDivision)
#option trace;
local LMF,Vars,p,r,flag,g,u,Mult,InvDiv,FF,j,i,qq,c1,m,M,A,B,QQQ,a,b,cd,LO;
global tord,T,Polys,ISR2,sy,C3,rdz;
LMF:=copy(Array([seq(v[7][2], v in T)])); 
Vars:=[op(TermOrder)];
p:=Polys[pp[1]]:                                 
r:=0;
FF:=copy(Array([seq(Polys[v[1]], v in T)]));
j:=0:m:=ArrayNumElems(FF);
qq:=NULL;
QQQ:=NULL;
A,B:=op(pp[7]);
while p<>0 do     
      j:=j+1:
      flag:=false:
      i:=1;
      while not flag and i<=m do
          g:=FF[i]:
      	  if divide(B,LMF[i],'u') then
	           Mult:=InvDivision(LMF[i],LMF,Vars);
	           if indets(u) subset {op(Mult)} then 
	           if j=1 then
	              cd:=[pp[6],[u*T[i][6][1],T[i][6][2]]];
	              M := sort(cd, proc (a, b) options operator, arrow; Sher1(a, b) end proc);
	              if M[-1]=pp[6] then
	               qq:=pp[6];  
	             fi;
	          fi; 
	     	   flag:=true:                   
               if j=1 then
                 if Crit1(pp,T[i],InvDivision) and pp[8][1]<>true then                 
                sy(ArrayNumElems(sy)+1):=[[Polys[T[i][2][1]]*Polys[pp[2][1]],Polys[pp[2][1]],pp[2][3]],[Polys[T[i][2][1]]*Polys[pp[2][1]],Polys[T[i][2][1]],T[i][2][3]]];
                   RETURN([0,[0,qq],[[Polys[pp[2][1]],T[i][2][1],i],[Polys[T[i][2][1]],pp[2][1],pp[2][1]]]]);  
                 fi:
               fi:
               QQQ:=QQQ,[A/T[i][7][1]*u,T[i][1],i];
		       p:=simplify(p-(A/T[i][7][1])*u*g);
		       A,B:=LeadingTerm(p,TermOrder);
	        else 
                 i:=i+1;
	        fi:
	     else 
            i:=i+1;
	  fi:
    od:
      if not flag  then 
            j:=j+1;
	        r:=simplify(r+A*B);
            p:=simplify(p-A*B);
            A,B:=LeadingTerm(p,TermOrder);
      fi:
od:
if r=0 then
sy(ArrayNumElems(sy)+1):=[seq([Polys[a[3]]*a[1],a[1],a[3]] , a in [QQQ]),[-Polys[pp[1]],pp[4],pp[8][3]]];
if pp[8][1]=true then      
   for t from 1 to ArrayNumElems(T) do
      ZS:=NULL;           
      for i from 1 to nops(T[t][8][2])-1 do 
          if pp[1]=T[t][8][2][i][2] then                                                                    
              ZS:=ZS,seq([tt[1]*T[t][8][2][i][1],tt[2],tt[3]], tt in [QQQ]);                                         
           else      
          ZS:=ZS,T[t][8][2][i];                        
          fi; 
        od;
       T[t][8][2]:= [ZS,T[t][8][2][-1]];
    od;  
 fi;
fi;

if r<>0 then
if pp[8][1]=true then      
   for t from 1 to ArrayNumElems(T) do
      ZS:=NULL;           
      for i from 1 to nops(T[t][8][2])-1 do 
          if pp[1]=T[t][8][2][i][2] then                     
              ZS:=ZS,seq([tt[1]*T[t][8][2][i][1],tt[2],tt[3]], tt in [QQQ]),[T[t][8][2][i][1],p[1],ArrayNumElems(T)+1];            
           else      
          ZS:=ZS,T[t][8][2][i];                        
          fi; 
        od;
       T[t][8][2]:= [ZS,T[t][8][2][-1]];
    od;  
 fi;
 fi;        
if r=0  then 
  rdz:=rdz+1: 
fi:        
RETURN([r,[[LeadingTerm(r,tord)],qq],[QQQ]]);
end:
################  Head AutoReduce #################
HAutoReduction:=proc(f,F,TermOrder,InvDivision)
#option trace;
local LMF,Vars,p,r,flag,g,u,Mult,InvDiv,i,rr,LMFF,A,B;
global tord;
LMF:=copy(Array([seq(LeadingMonomial(v,TermOrder), v in F)]));
Vars:=[op(TermOrder)];
p:=LM(f):
rr:=ArrayNumElems(F);
LMFF:=copy(Array([seq(v, v in LMF),LM(f)]));
A,B:=LeadingTerm(p,TermOrder);
FLAG := false;
while FLAG = false and p<>0 do
      flag:=false:i:=1;
      FLAG := true;
      while not flag and i<=rr do
          g:=F[i]:
      	  if divide(B,LMF[i],'u') then
	         Mult:=InvDivision(LMF[i],LMFF,Vars);
	         if indets(u) subset {op(Mult)} then
	     	 flag:=true:
	     	 FLAG := false;
	     	 p:=simplify(p-(A/LeadingCoefficient(g,TermOrder))*u*g);
	     	 A,B:=LeadingTerm(p,TermOrder);
	      else i:=i+1;
	     fi:
	   else i:=i+1;
	   fi:
     od:
od:
RETURN(p);
end:                
###############################
HAutoReduce:=proc(F,TermOrder,InvDivision)
local InvDiv,H,i;
#option trace;
H:=copy(F);
for i from 1 to ArrayNumElems(F) do
    H(i):=HAutoReduction(H[i],Array([op({seq(h , h in H)} minus {H[i],0})]),TermOrder,InvDivision);
od:
RETURN(Array([op({seq(h, h in H)} minus {0})]));
end:                 

################ Normal Form Algorithm2  ###############
with(Groebner):
Normalform2:=proc(pp,TermOrder,InvDivision)
#option trace;
local LMF,Vars,p,r,flag,g,u,Mult,InvDiv,FF,j,i,qq,c1,m,M,A,B,QQQ,a,b,cd,LO;
global tord,T,Polys,ISR2,C3,rdz;
LMF:=copy(Array([seq(v[7][2], v in T)])); 
Vars:=[op(TermOrder)];
p:=Polys[pp[1]]:                                 
r:=0;
FF:=copy(Array([seq(Polys[v[1]], v in T)]));
j:=0:m:=ArrayNumElems(FF);
qq:=NULL;
QQQ:=NULL;
A,B:=op(pp[7]);
while p<>0 do     
      j:=j+1:
      flag:=false:
      i:=1;
      while not flag and i<=m do
          g:=FF[i]:
      	  if divide(B,LMF[i],'u') then
	           Mult:=InvDivision(LMF[i],LMF,Vars);
	           if indets(u) subset {op(Mult)} then
	           if j=1 and pp[6][2]=T[i][6][2] and pp[6][1]=u*T[i][6][1] then
	              C3:=C3+1;
	              RETURN([0,[0,[1,0]]]);
	           fi;
	           if j=1 then
	              cd:=[pp[6],[u*T[i][6][1],T[i][6][2]]];
	              M := sort(cd, proc (a, b) options operator, arrow; Sher1(a, b) end proc);
	              if M[-1]=pp[6] then
	               qq:=pp[6];
	             fi;
	          fi; 
	     	   flag:=true:                   
               if j=1 then
                 if Crit1(pp,T[i],InvDivision) then 
                   RETURN([0,[0,[1,0]],[[Polys[pp[2][1]],T[i][2][1],i],[Polys[T[i][2][1]],pp[2][1],pp[2][1]]]]);  
                 fi:
               fi:
             if j=1 then
                 if Crit2(pp,T[i],InvDivision) then 
                     RETURN([0,[0,[1,0]]]);
                 fi:
               fi:
		       p:=simplify(p-(A/T[i][7][1])*u*g);
		       A,B:=LeadingTerm(p,TermOrder);
	        else 
                 i:=i+1;
	        fi:
	     else 
            i:=i+1;
	  fi:
    od:
      if not flag  then 
            j:=j+1;
	        r:=simplify(r+A*B);
            p:=simplify(p-A*B);
            A,B:=LeadingTerm(p,TermOrder);
      fi:
od:
if r=0  then 
  rdz:=rdz+1: 
 fi:
RETURN([r,[[LeadingTerm(r,tord)],qq]]);
end:
###################################################
#                Involutive Basis                 #
###################################################
InvBasis:=proc(Ideal,Tord,InvDivision)
#option trace;
global tord,T,Q,Polys,Vars,c1,c2,rdz,LMF,critelim, reduct,isr,deg,LLMI, mmtelim,sy;
local firsttime, firstbytes, cov, F, n, syz,  p, FLAGE, j, hh, h, k, rpp, t, i, SS, S, s, TT, jj, fl, c, q, NM, y, Q2, Q3, ii0;
F:=Array(1..nops(Ideal),Ideal);
F:=copy(sort(F,(a,b)->TestOrder(a,b,Tord)));
n:=ArrayNumElems(F):
syz:=Array([]);sy:=Array([]);
LLMI:=Array(1..n,[seq([LeadingTerm(i,tord)], i in F)]);
deg:=0;
Polys := copy(Array(1 .. n, proc (i) options operator, arrow; F[i] end proc));
T:=copy(Array(1..1,[[1,[1,LLMI[1][2],1],{},1,1,[1,1],LLMI[1],[false,[[1,1,1],[-Polys[1],1,1]]]]]));
Q:=copy(Array(1..ArrayNumElems(F)-1,[seq([i,[i,LLMI[i][2],i],{},1,i,[1,i],LLMI[i],[false,1,i]],i=2..n)]));
Q := sort(Q, proc (a, b) options operator, arrow; Sher22(a, b) end proc);
while ArrayNumElems(Q) <> 0 do 
  p := Q[1];                                           
  deg := max(deg, degree(Polys[p[1]])); 
  Q := copy(Array(1 .. ArrayNumElems(Q)-1, proc(t) options operator, arrow; Q[t+1] end proc)); 
         FLAGE := false;             
       if FLAGE = false then
           hh:=Normalform(p,tord,InvDivision);
           h:=hh[1]: 
           k := ArrayNumElems(T)+1;
          rpp:= Polys[p[1]];
 
          if nops(hh[2])>1 then                     
          if hh[2][2][1]<>1 then            
            syz(ArrayNumElems(syz)+1):=hh[2][2]; 
         fi;  
         fi;                
        if h=0 and p[1]=p[2][1] and p[8][1]<>true  then
	    SS:=Array([]);
            S:=copy(Q):
             for s in S do 
	           if s[2][1]<>p[1] then	                   
	       	     SS(ArrayNumElems(SS)+1):=s;
		    else
		     mmtelim :=mmtelim+1:                                                              
                  fi;
              od:
             Q:=copy(SS);
         elif h<>0 and hh[2][1][2]<>p[7][2] then 
              k := ArrayNumElems(T)+1;
	          Polys[p[1]]:=h; 
	          T(k):=[p[1],[p[1],hh[2][1][2],k],{},1,p[1],[1,p[1]],hh[2][1],[false, [[1,p[1],k],seq(t, t in  hh[3]),[-rpp,p[4],p[8][3]]] ]];
              TT:=Array([]);            
                for k from 1 to ArrayNumElems(T) do
	               if divide(T[k][7][2], hh[2][1][2]) and hh[2][1][2] <> T[k][7][2] then
	      	           Q(ArrayNumElems(Q)+1):=[seq(T[k][i], i=1..7),[true,T[i][8][2][-1][2],T[i][8][2][-1][3]]];
	      	           for jj from 1 to ArrayNumElems(T) do
	      	                   T[jj][2][3]:=T[jj][2][3]-1;
                               fl:=[seq( T[jj][8][2][g][3] , g=1..nops( T[jj][8][2]))];
	      	             for c from 1 to nops(fl) do
	      	    	       if  fl[c] >= ArrayNumElems(TT)+1  then	      	          
	      	                 T[jj][8][2][c][3]:=T[jj][8][2][c][3]-1;	      	          
	      	               fi;	      	   
	      	             od;	      	            
	      	         od;	      	           
                       else 
		           TT(ArrayNumElems(TT)+1):=T[k]; 
                       fi:                
                od:               
                 T:=copy(TT);
                   for i from 1 to ArrayNumElems(T) do
       	               q:=T[i];
       	               NM:={op(Vars)} minus {op(InvDivision(q[7][2],Array([seq(w[7][2], w in T)]),Vars))};
	                       for y in NM minus q[3] do
		                   n:=n+1: 
	                           Polys(n):=expand(y*Polys[q[1]]);	                       
	                           Q(ArrayNumElems(Q)+1):=[n,q[2],{},y,q[1],[expand(y*q[6][1]),q[6][2]],[q[7][1],y*q[7][2]],[false,y,q[8][2][1][3]]];	                            
	   	                   od:
                            T[i][3]:=q[3] union NM;
                      od:
           elif h<>0 and hh[2][1][2]=p[7][2] then
	           Polys[p[1]]:=h;
	          
	           T(k):=[seq(p[i],i=1..7),[false, [[1,p[1],k],seq(t , t in hh[3]),[-rpp,p[4],p[8][3]]]]];
	             
                    for i from 1 to ArrayNumElems(T) do
       	                 q:=T[i];
                 	 NM:={op(Vars)} minus {op(InvDivision(q[7][2],Array([seq(w[7][2], w in T)]),Vars))};
	                 for y in NM minus q[3] do 
		             n:=n+1: 
	                    Polys(n):=expand(y*Polys[q[1]]);
	                    Q(ArrayNumElems(Q)+1):=[n,q[2],{},y,q[1],[expand(y*q[6][1]),q[6][2]],[q[7][1],y*q[7][2]],[false,y,q[8][2][1][3]]];                      
	   	         od:
                         T[i][3]:=q[3] union NM;
                       od; 
              fi;  
         Q:=convert(Q,list); 
         Q := sort(Q, proc (a, b) options operator, arrow; Sher22(a, b) end proc);      
         Q:=copy(Array(1..nops(Q),Q));              
   fi;   
od:
SY:=Array(1..ArrayNumElems(T), [seq(tt[8][2] , tt in T)] );
SY1:=Array([]);
for i in SY do 
   if i[-1][2]<>1 then
      SY1(ArrayNumElems(SY1)+1):= [seq([Polys[i[ii][2]]*i[ii][1],i[ii][1],i[ii][3]],ii=1..nops(i)-1), i[-1]];                                                           
   fi;
 od;
SYY:=Array(1..ArrayNumElems(sy)+ArrayNumElems(SY1) ,[seq(a , a in SY1),seq( b ,b in sy)]);                                                                                                     
TH:=Array(1..ArrayNumElems(T),[seq([1,tt[7],1], tt in T)]):
THH:=Array(1..ArrayNumElems(T),[seq(tt[7][2], tt in T)]):
T:=Array(1..ArrayNumElems(T),[seq([Polys[T[tt][1]],[1,tt]], tt=1..ArrayNumElems(T))]);
RETURN([T,TH,THH,SYY]);
end:
################  Involutive Basis 2 ########################
InvBasis2:=proc(Ideal,Tord,InvDivision,SYZ)
#option trace;
global tord,T,Q,Polys,Vars,c1,c2,rdz,LMF,critelim, reduct,isr,deg,LLMI,syzelim,HE, mmtelim;
local firsttime, firstbytes, cov, F, n, syz, sy, p, FLAGE, j, hh, h, k, rpp, t, i, SS, S, s, TT, jj, fl, c, q, NM, y, Q2, Q3, ii0;
F:=Ideal;
F:=copy(sort(F,(a,b)->TestOrder(a[1],b[1],tord)));
n:=ArrayNumElems(F):
syz:=Array([]);sy:=Array([]);
LLMI:=[seq([LeadingTerm(i[1],tord)], i in F)];
deg:=0;
Polys := copy(Array(1 .. n, proc (i) options operator, arrow; F[i][1] end proc));
T:=copy(Array(1..1,[[1,[1,LLMI[1][2]],{},1,1,[1,1],LLMI[1],F[1][2]]]));
Q:=copy(Array(1..ArrayNumElems(F)-1,[seq([i,[i,LLMI[i][2]],{},1,i,[1,i],LLMI[i],F[i][2]],i=2..n)]));
Q := sort(Q, proc (a, b) options operator, arrow; Sher22(a, b) end proc);
while ArrayNumElems(Q) <> 0 do 
  p := Q[1];     
  Q := copy(Array(1 .. ArrayNumElems(Q)-1, proc(t) options operator, arrow; Q[t+1] end proc));                                        
  deg := max(deg, degree(Polys[p[1]])); 
  flage:=false;  
    for j from 1 to ArrayNumElems(SYZ) while flage = false  do
     if  p[8][2]=SYZ[j][2] and  divide(p[8][1],SYZ[j][1]) = true then
        flage := true;     
        syzelim := syzelim+1; 
     fi;
   end do;            
  if flage=false then              
        FLAGE3:=false;
         if ArrayNumElems(syz) <> 0 then
             for j from 1 to ArrayNumElems(syz) while FLAGE3 = false  do
                 if  p[6][2]=syz[j][2] and  divide(p[6][1],syz[j][1],'uu') = true and uu<>1 then 
                     FLAGE3 := true;  
                     mmtelim := mmtelim+1;    # print(p[7]); 
                  fi;
             end do;
           end if;                                                                  
     if  FLAGE3=false then
              flage1:=false;
     if HFI(T,deg,Tord,Vars,Janet)=MMM(deg) then
        HE:=HE+1;                
        flage1:=true:                  
     fi;               
       if flage1 = false then
           hh:=Normalform2(p,tord,InvDivision);  
           h:=hh[1]: 
           k := ArrayNumElems(T)+1;
           rpp:= Polys[p[1]];         
          if nops(hh[2])>1 then                     
          if hh[2][2][1]<>1 then           
            syz(ArrayNumElems(syz)+1):=hh[2][2]; 
         fi;  
         fi;                 

        if h=0 and p[1]=p[2][1] then
	    SS:=Array([]);
            S:=copy(Q):
             for s in S do 
	           if s[2][2]*s[2][1]<>p[1] then	                   
	       	     SS(ArrayNumElems(SS)+1):=s;
		    else
		      isr:=isr+1:
                  fi;
              od:
             Q:=copy(SS);
         elif h<>0 and hh[2][1][2]<>p[7][2] then 
              k := ArrayNumElems(T)+1; 
	          Polys[p[1]]:=h; 
	          T(k):=[p[1],[p[1],hh[2][1][2]],{},1,p[1],[1,p[1]],hh[2][1],p[8]];
              TT:=Array([]);            
                for k from 1 to ArrayNumElems(T) do
	               if divide(T[k][7][2], hh[2][1][2]) and hh[2][1][2] <> T[k][7][2] then
	      	            Q(ArrayNumElems(Q)+1):=T[k];	      	     	      	           
                       else 
		           TT(ArrayNumElems(TT)+1):=T[k]; 
                       fi:                
                od:               
                 T:=copy(TT);
                   for i from 1 to ArrayNumElems(T) do
       	               q:=T[i];
       	               NM:={op(Vars)} minus {op(InvDivision(q[7][2],Array([seq(w[7][2], w in T)]),Vars))};
	                       for y in NM minus q[3] do
		                   n:=n+1: 
	                           Polys(n):=expand(y*Polys[q[1]]);	                       
	                           Q(ArrayNumElems(Q)+1):=[n,q[2],{},y,q[1],[expand(y*q[6][1]),q[6][2]],[q[7][1],y*q[7][2]],[expand(y*q[8][1]),q[8][2]]];	                            
	   	                   od:
                            T[i][3]:=q[3] union NM;
                      od:
           elif h<>0 and hh[2][1][2]=p[7][2] then
	           Polys[p[1]]:=h;
	           T(k):=p;	             
                    for i from 1 to ArrayNumElems(T) do
       	                 q:=T[i];
                 	 NM:={op(Vars)} minus {op(InvDivision(q[7][2],Array([seq(w[7][2], w in T)]),Vars))};
	                 for y in NM minus q[3] do 
		            n:=n+1: 
	                    Polys(n):=expand(y*Polys[q[1]]);
	                    Q(ArrayNumElems(Q)+1):=[n,q[2],{},y,q[1],[expand(y*q[6][1]),q[6][2]],[q[7][1],y*q[7][2]],[expand(y*q[8][1]),q[8][2]]];                      
	   	         od:
                         T[i][3]:=q[3] union NM;
             od; 
         fi;  
         Q:=convert(Q,list); 
         Q := sort(Q, proc (a, b) options operator, arrow; Sher22(a, b) end proc);      
         Q:=copy(Array(1..nops(Q),Q));              
   fi; 
   fi;
   fi;
od:
THH:=Array(1..ArrayNumElems(T),[seq(T[tt][7][2], tt=1..ArrayNumElems(T))]):
T:=Array(1..ArrayNumElems(T),[seq([Polys[T[tt][1]],[1,T[tt][8][2]]], tt=1..ArrayNumElems(T))]);
RETURN([T,THH]);
end:
#############################
#       Hilbert function    #   
#############################
 binom:=proc(n,r) 
  if n>=r and r>=0 then
     return(binomial(n,r));
  else
    return(0);
    fi;
 end:
  
########################
HFI:=proc(F,s,T,V,invdiv)
local n,G,K,H;
#option trace;
n:=nops(V);  
G:=copy(Array([seq(g[-2][2],g in F)]));
K:= [seq(nops({op(invdiv(g,G,V))}),g in G)]; 
H:=add(i , i in [seq(binom(s-degree(G[i])+K[i]-1,s-degree(G[i])),i=1..ArrayNumElems(G))]); 
return(H);  
end: 
################################
 #          TEST               #
################################
 TEST := proc (u, F, V,Tord)
local M,VV, k,i,Ind,L,lu;
#option trace;
L:=F;
lu:=LM(u);
M[1] := {op(Janet(lu, L,V))};
M[2] := {op(RightPommaret(lu, L, V))};
if nops(M[1]) = nops(M[2]) then
 return (true);
 else 
   VV:= sort([op(`minus`(M[1], M[2]))], proc (a, b) options operator, arrow; TestOrder(b, a, Tord) end proc);  
	Ind:=indets(u):
	for i from  nops(V) to 1 by -1 while not (V[i] in Ind) do
	od: 
    return (false, VV[1],V[i]);   
end if;
end proc:
###########################
GTEST := proc (F, V,Tord)
local u, A; 
for u in F do 
 A := TEST(u, F,V, Tord);
  if nops([A]) <> 1 then 
  RETURN(A);
  end if;
end do;
RETURN(true);
end proc:

############################## 
#       Pommaert Basis       #
############################## 
POM := proc (F, Tord,InDiv) 
local G, L, chen, A, u, LL, LLL,i, B, K,V,uu,m,secondtime,firsttime,secondbytes,firstbytes,sj,LR,ltsj,sjj,KK,LH,LS,LLH; 
global tord,T,Q,Polys,Vars,c1,c2,c3,deg,rdz, reduct,ltshh, syzelim, mmtelim,MMM,HE,isr,ISR2,C3;
#option trace;
firsttime, firstbytes := kernelopts(cputime, bytesused);
HE:=0;
deg:=0:
ISR2:=0;
reduct := 0;
c1:=0:c2:=0:rdz:=0;C3:=0;
mmtelim := 0;
syzelim := 0;
isr:=0:
tord:=Tord;
Vars:=[seq([op(Tord)][i], i = 1 .. nops([op(Tord)]))];
LH:=copy(InvBasis(F, Tord,Janet));
L :=copy(LH[1]);               
LLL:=copy(LH[2]);
SYZ:=LH[4];
MMM:=(d->HFI(LLL,d,Tord,Vars,Janet));
chen := NULL;
LR:=copy(LH[3]);
A := GTEST(LR,Vars,Tord);
if A=true then  
LL:=L;
fi; 
while A <> true do
   	  uu := rand(-2 .. 2)(); 
      while uu = 0 do 
         uu := rand(-2 .. 2)(); 
      end do;
      m := A[3] = uu*A[2]+A[3];
      SYZZ:=expand(subs(m, SYZ));
      LS:=expand(subs(m, L));
      TSYZ:=copy(Array(1..ArrayNumElems(SYZ),[seq(LTSY(f), f in SYZZ)]));  
      Sk:={seq(l , l in TSYZ)}; # print(Sk);
      LTSYZ:=copy(Array(1..nops(Sk),[op(Sk)]));
      LLH:=copy(InvBasis2(LS,Tord,Janet,LTSYZ));   
      LL:=copy(LLH[1]);                             
      LR:=copy(LLH[2]):
      LR:=HAutoReduce(LR,Tord,RightPommaret);
      B := GTEST(LR, Vars,Tord);
      if B <> A then        
         L:=LS;
         A := B;
         SYZ:=SYZZ;
         chen := chen, m;
      end if;
  end do;

 secondtime, secondbytes := kernelopts(cputime, bytesused); 
printf("%-1s %1s %1s %1s:         %3a %3a\n", The, cpu, time, is, secondtime-firsttime, sec);
printf("%-1s %1s %1s:         %3a %3a\n", The, used, memory, secondbytes-firstbytes, bytes);
printf("%-1s %1s %1s             : %3g\n",Num,of,Rds,rdz):
printf("%-1s %1s %1s        : %3a\n",Num,of,criteria,[c1,c2,C3]):
printf("%-1s %1s %1s        : %3a\n",Max,of,degree,deg):
printf("%-1s %1s %1s        : %3a\n",Num,of,linearchange , nops([chen])):
printf("%-1s %1s %1s %1s %1s %1s:    %3a\n", The, number, of, MMT, elimination, is, mmtelim);
printf("%-1s %1s %1s %1s %1s %1s:    %3a\n", The, number, of, syzygy , elimination, is, syzelim);
printf("%-1s %1s %1s %1s %1s %1s:    %3a\n", The, number, of, HilbertDriven, elimination, is, HE);
return ([seq(LM(l[1]), l in LL)],chen);
end:	
##############################################
##############################################









##############################################
#  Involutive Variant of the GVW Algorithm   #
##############################################

##############################################
#            Involutive Divisions            #
##############################################
###########Janet Division#######
Janet:=proc(u,U,Vars)
local n,m,d,L,i,j,dd,V,v,Mult;
#option trace;
if degree(u)=0 then RETURN(Vars); fi;
n:=nops(Vars);
m:=ArrayNumElems(U);
d:=[seq(max(seq(degree(U[j],Vars[i]),j=1..m)),i=1..n)];
Mult:=NULL;
    if degree(u,Vars[1])=d[1] then
       Mult:=Mult,Vars[1];
    fi:
for j from 2 to n do
    	dd:=[seq(degree(u,Vars[l]),l=1..j-1)];
    	V:=NULL:
    	for v in U do
    	    if [seq(degree(v,Vars[l]),l=1..j-1)]=dd then
	       V:=V,v:
	    fi:
	od:
	if degree(u,Vars[j])=max(seq(degree(v,Vars[j]), v in [V])) then
	   Mult:=Mult,Vars[j];
	fi:
od:
RETURN([Mult]);
end:

########### Pommaret Division #############
LeftPommaret:=proc(u,U,Vars)
local N,Ind,i;
N:=NULL:
Ind:=indets(u):
for i from 1 to nops(Vars) while not (Vars[i] in Ind) do
    N:=N,Vars[i]:
od:
N:=N,Vars[i]:
RETURN([N]);
end:
    	   

RightPommaret:=proc(u,U,Vars)
local N,Ind,i;
N:=NULL:
Ind:=indets(u):
for i from  nops(Vars) by -1 to 1 while not (Vars[i] in Ind) do
    N:=N,Vars[i]:
od:
N:=N,Vars[i]:
RETURN([N]);
end:
#################################
#     Schreyer’s ordering test   #
#################################
with(Groebner):
Sher22:=proc(p,q)
#option trace;##TOP3
local i,j,a,b;
global tord,LLMI;
 i := p[2];
 j := q[2];  
 a:=p[1]*LLMI[i][2];
 b:=q[1]*LLMI[j][2];
if p[1]<>q[1] then
RETURN(TestOrder(p[1],q[1],tord));
else
    RETURN(evalb(i<j));
   fi;
end:
#################################
ord:=proc(p,q)
#option trace;##TOP3
local i,j,a,b;
global tord,LR,Vars;
if p[1]<>q[1] then
RETURN(TestOrder(p[1],q[1],tord));
else
    RETURN(evalb(p[2]<q[2]));
   fi;
end:
#########################
 LTSY := proc (f)
# option trace;
 local A,B;
 global tord;
 A:= sort(f, proc (a, b) options operator, arrow; ord(a, b) end proc)[-1];
  return([LM(A[1]),A[2]]);
end proc:
####################################
#         Reduction  CRITERIA      # 
####################################
Crit1:=proc(p,gg)
global c1;
local InvDiv,flag,g,l;
#option trace;
l:=p[2][2];
g:=gg[2][2];
    if l*g=p[7][2] then 
    c1:=c1+1; 
    RETURN(true);
fi;
RETURN(false);
end:
###########################
Crit2:=proc(p,gg)
global c2;
local g,l,lc;
#option trace;
l:=p[2][2];
g:=gg[2][2];
lc:=lcm(l,g);
if  divide(p[7][2],lc) and p[7][2]<>lc  then
    c2:=c2+1; 
    RETURN(true);
fi;
RETURN(false);
end:
############################################## 
#           Normal Form Algorithm            #
##############################################
with(Groebner):
Normalform:=proc(pp,TermOrder,InvDivision)
#option trace;
local LMF,Vars,p,r,flag,g,u,Mult,InvDiv,FF,j,i,qq,c1,m,M,A,B,QQQ,a,b,cd,LO;
global tord,T,Polys,ISR2,sy,C3,rdz;
LMF:=copy(Array([seq(v[7][2], v in T)])); 
Vars:=[op(TermOrder)];
p:=Polys[pp[1]]:                                 
r:=0;
FF:=copy(Array([seq(Polys[v[1]], v in T)]));
j:=0:m:=ArrayNumElems(FF);
qq:=NULL;
QQQ:=NULL;
A,B:=op(pp[7]);
while p<>0 do     
      j:=j+1:
      flag:=false:
      i:=1;
      while not flag and i<=m do
          g:=FF[i]:
      	  if divide(B,LMF[i],'u') then
	           Mult:=InvDivision(LMF[i],LMF,Vars);
	           if indets(u) subset {op(Mult)} then 
	           if j=1 then
	              cd:=[pp[6],[u*T[i][6][1],T[i][6][2]]];
	              M := sort(cd, proc (a, b) options operator, arrow; Sher1(a, b) end proc);
	              if M[-1]=pp[6] then
	               qq:=pp[6];  
	             fi;
	          fi; 
	     	   flag:=true:                   
               if j=1 then
                 if Crit1(pp,T[i],InvDivision) and pp[8][1]<>true then                 
                sy(ArrayNumElems(sy)+1):=[[Polys[T[i][2][1]]*Polys[pp[2][1]],Polys[pp[2][1]],pp[2][3]],[Polys[T[i][2][1]]*Polys[pp[2][1]],Polys[T[i][2][1]],T[i][2][3]]];
                   RETURN([0,[0,qq],[[Polys[pp[2][1]],T[i][2][1],i],[Polys[T[i][2][1]],pp[2][1],pp[2][1]]]]);  
                 fi:
               fi:
               QQQ:=QQQ,[A/T[i][7][1]*u,T[i][1],i];
		       p:=simplify(p-(A/T[i][7][1])*u*g);
		       A,B:=LeadingTerm(p,TermOrder);
	        else 
                 i:=i+1;
	        fi:
	     else 
            i:=i+1;
	  fi:
    od:
      if not flag  then 
            j:=j+1;
	        r:=simplify(r+A*B);
            p:=simplify(p-A*B);
            A,B:=LeadingTerm(p,TermOrder);
      fi:
od:
if r=0 then
sy(ArrayNumElems(sy)+1):=[seq([Polys[a[3]]*a[1],a[1],a[3]] , a in [QQQ]),[-Polys[pp[1]],pp[4],pp[8][3]]];
if pp[8][1]=true then      
   for t from 1 to ArrayNumElems(T) do
      ZS:=NULL;           
      for i from 1 to nops(T[t][8][2])-1 do 
          if pp[1]=T[t][8][2][i][2] then                                                                    
              ZS:=ZS,seq([tt[1]*T[t][8][2][i][1],tt[2],tt[3]], tt in [QQQ]);                                         
           else      
          ZS:=ZS,T[t][8][2][i];                        
          fi; 
        od;
       T[t][8][2]:= [ZS,T[t][8][2][-1]];
    od;  
 fi;
fi;

if r<>0 then
if pp[8][1]=true then      
   for t from 1 to ArrayNumElems(T) do
      ZS:=NULL;           
      for i from 1 to nops(T[t][8][2])-1 do 
          if pp[1]=T[t][8][2][i][2] then                     
              ZS:=ZS,seq([tt[1]*T[t][8][2][i][1],tt[2],tt[3]], tt in [QQQ]),[T[t][8][2][i][1],p[1],ArrayNumElems(T)+1];            
           else      
          ZS:=ZS,T[t][8][2][i];                        
          fi; 
        od;
       T[t][8][2]:= [ZS,T[t][8][2][-1]];
    od;  
 fi;
 fi;        
if r=0  then 
  rdz:=rdz+1: 
fi:        
RETURN([r,[[LeadingTerm(r,tord)],qq],[QQQ]]);
end:
#########################################
#      LEADIN MONOMIAL                  #
#########################################
with(PolynomialIdeals):
#################
LM:=proc(f)
global tord;
if f<>0 then
  RETURN(LeadingMonomial(f,tord));
fi:
RETURN(0);
end:

##################
 LT := proc (f)
 local A,B;
 global tord;
 if f=0 then return(0); fi;
 A,B:=LeadingTerm(p,tord);
  return(A*B);
end proc:
#################
LC:=proc(f)
global tord;
    RETURN(LeadingCoefficient(f,tord));
end: 
############################################
#         Cover function                   #
############################################
Cover := proc (p, q,m) 
global tord;
local v,u,t;
#option trace:
u := p[1]; 
v := q[1];  
if not evalb(u[2] = v[2]) then 
   RETURN(false): 
 elif divide(u[1], v[1],'k') then 
   t := k*q[2]; 
   if TestOrder(t, p[2],tord) and t <> p[2] then 
       RETURN(true):
   end if:
end if; 
RETURN(false):
end proc:
#####################################################
#             INVOLUTIVE BASIS                      #
#####################################################
InvBasis:=proc(Ideal,Tord,InvDivision)
#option trace;
global tord,T,Q,Polys,Vars,c1,c2,rdz,LMF,critelim, reduct,isr,deg,fff,LLMI, mmtelim,sy;
local firsttime, firstbytes, cov, F, n, syz,  p, FLAGE, j, hh, h, k, rpp, t, i, SS, S, s, TT, jj, fl, c, q, NM, y, Q2, Q3, ii0;
F:=Array(1..nops(Ideal),Ideal);
F:=copy(sort(F,(a,b)->TestOrder(a,b,Tord)));
n:=ArrayNumElems(F):
syz:=Array([]);sy:=Array([]);
LLMI:=Array(1..n,[seq([LeadingTerm(i,tord)], i in F)]);
deg:=0;
Polys := copy(Array(1 .. n, proc (i) options operator, arrow; F[i] end proc));
T:=copy(Array(1..1,[[1,[1,LLMI[1][2],1],{},1,1,[1,1],LLMI[1],[false,[[1,1,1],[-Polys[1],1,1]]]]]));
Q:=copy(Array(1..ArrayNumElems(F)-1,[seq([i,[i,LLMI[i][2],i],{},1,i,[1,i],LLMI[i],[false,1,i]],i=2..n)]));
Q := sort(Q, proc (a, b) options operator, arrow; Sher22(a, b) end proc);
while ArrayNumElems(Q) <> 0 do 
  p := Q[1];                                           
  deg := max(deg, degree(Polys[p[1]])); 
  Q := copy(Array(1 .. ArrayNumElems(Q)-1, proc(t) options operator, arrow; Q[t+1] end proc)); 
         FLAGE := false;             
       if FLAGE = false then
           hh:=Normalform(p,tord,InvDivision);
           h:=hh[1]: 
           k := ArrayNumElems(T)+1;
          rpp:= Polys[p[1]];
 
          if nops(hh[2])>1 then                     
          if hh[2][2][1]<>1 then            
            syz(ArrayNumElems(syz)+1):=hh[2][2]; 
         fi;  
         fi;                

        fff:=0;
        if h=0 and p[1]=p[2][1] and p[8][1]<>true  then
	    SS:=Array([]);
            S:=copy(Q):
             for s in S do 
	           if s[2][1]<>p[1] then	                   
	       	     SS(ArrayNumElems(SS)+1):=s;
		    else
		     mmtelim :=mmtelim+1:                                                              
                  fi;
              od:
             Q:=copy(SS);
         elif h<>0 and hh[2][1][2]<>p[7][2] then 
              k := ArrayNumElems(T)+1;
	          Polys[p[1]]:=h; 
	          T(k):=[p[1],[p[1],hh[2][1][2],k],{},1,p[1],[1,p[1]],hh[2][1],[false, [[1,p[1],k],seq(t, t in  hh[3]),[-rpp,p[4],p[8][3]]] ]];
              TT:=Array([]);            
                for k from 1 to ArrayNumElems(T) do
	               if divide(T[k][7][2], hh[2][1][2]) and hh[2][1][2] <> T[k][7][2] then
	      	           Q(ArrayNumElems(Q)+1):=[seq(T[k][i], i=1..7),[true,T[i][8][2][-1][2],T[i][8][2][-1][3]]];
	      	           for jj from 1 to ArrayNumElems(T) do
	      	                   T[jj][2][3]:=T[jj][2][3]-1;
                               fl:=[seq( T[jj][8][2][g][3] , g=1..nops( T[jj][8][2]))];
	      	             for c from 1 to nops(fl) do
	      	    	       if  fl[c] >= ArrayNumElems(TT)+1  then	      	          
	      	                 T[jj][8][2][c][3]:=T[jj][8][2][c][3]-1;	      	          
	      	               fi;	      	   
	      	             od;	      	            
	      	         od;	      	           
                       else 
		           TT(ArrayNumElems(TT)+1):=T[k]; 
                       fi:                
                od:               
                 T:=copy(TT);
                   for i from 1 to ArrayNumElems(T) do
       	               q:=T[i];
       	               NM:={op(Vars)} minus {op(InvDivision(q[7][2],Array([seq(w[7][2], w in T)]),Vars))};
	                       for y in NM minus q[3] do
		                   n:=n+1: 
	                           Polys(n):=expand(y*Polys[q[1]]);	                       
	                           Q(ArrayNumElems(Q)+1):=[n,q[2],{},y,q[1],[expand(y*q[6][1]),q[6][2]],[q[7][1],y*q[7][2]],[false,y,q[8][2][1][3]]];	                            
	   	                   od:
                            T[i][3]:=q[3] union NM;
                      od:
           elif h<>0 and hh[2][1][2]=p[7][2] then
	           Polys[p[1]]:=h;
	          
	           T(k):=[seq(p[i],i=1..7),[false, [[1,p[1],k],seq(t , t in hh[3]),[-rpp,p[4],p[8][3]]]]];
	             
                    for i from 1 to ArrayNumElems(T) do
       	                 q:=T[i];
                 	 NM:={op(Vars)} minus {op(InvDivision(q[7][2],Array([seq(w[7][2], w in T)]),Vars))};
	                 for y in NM minus q[3] do 
		             n:=n+1: 
	                    Polys(n):=expand(y*Polys[q[1]]);
	                    Q(ArrayNumElems(Q)+1):=[n,q[2],{},y,q[1],[expand(y*q[6][1]),q[6][2]],[q[7][1],y*q[7][2]],[false,y,q[8][2][1][3]]];                      
	   	         od:
                         T[i][3]:=q[3] union NM;
                       od; 
              fi;  
         Q:=convert(Q,list); 
         Q := sort(Q, proc (a, b) options operator, arrow; Sher22(a, b) end proc);      
         Q:=copy(Array(1..nops(Q),Q));              
   fi;   
od:fff:=0;
SY:=Array(1..ArrayNumElems(T), [seq(tt[8][2] , tt in T)] );
SY1:=Array([]);
for i in SY do 
   if i[-1][2]<>1 then
      SY1(ArrayNumElems(SY1)+1):= [seq([Polys[i[ii][2]]*i[ii][1],i[ii][1],i[ii][3]],ii=1..nops(i)-1), i[-1]];                                                           
   fi;
 od;
SYY:=Array(1..ArrayNumElems(sy)+ArrayNumElems(SY1) ,[seq(a , a in SY1),seq( b ,b in sy)]);                                                                                                     
TH:=Array(1..ArrayNumElems(T),[seq([1,tt[7],1], tt in T)]):
THH:=Array(1..ArrayNumElems(T),[seq(tt[7][2], tt in T)]):
T:=Array(1..ArrayNumElems(T),[seq([Polys[T[tt][1]],[1,tt]], tt=1..ArrayNumElems(T))]);
RETURN([T,TH,THH,SYY]);
end:
###########  TEST IS INVOLUTIVE #################
with(Groebner):
Normalform1:=proc(pp,T,TermOrder,InvDivision)
#option trace;
local LMF,Vars,p,r,flag,g,u,Mult,InvDiv,FF,j,i,qq,c1,LO,Ma,b,cd,QQQ,m,A,B,M;
global tord,Polys,ISR2,fff;
LO:=[9];
LMF:=copy(Array([op({seq(LeadingMonomial(v,TermOrder), v in T)} minus {0})]));
Vars:=[op(TermOrder)];
p:=pp:                                 
r:=0;
FF:=copy(Array([op({seq(v, v in T)} minus {0})]  ));
j:=0:
m:=ArrayNumElems(FF);
qq:=NULL;
A,B:=LeadingTerm(p,TermOrder);
while p<>0 do     
      j:=j+1:
      flag:=false:
      i:=1;
      while not flag and i<=m do 
          g:=FF[i]:
      	  if divide(B,LM(g),'u') then
	         Mult:=InvDivision(LM(g),LMF,Vars);
	         if indets(u) subset {op(Mult)} then
	           
	     	  flag:=true:    
		         p:=simplify(p-(A/LeadingCoefficient(g,TermOrder))*u*g);
		         A,B:=LeadingTerm(p,TermOrder);
	           else 
                 i:=i+1;
	         fi:
	      else 
            i:=i+1;
	     fi:
      od:
      if not flag  then 
            j:=j+1;
	    r:=simplify(r+A*B);
            p:=simplify(p-A*B);
            A,B:=LeadingTerm(p,TermOrder);
      fi:
od:
RETURN(r);
end:
####################################
IsInvolutive:=proc(F,L,TermOrder,InvDivision)
local InvDiv,Vars,LMF,f,NM,u,G,GG;
#option trace;
Vars:=[op(TermOrder)];
G:=Basis(convert(F,list), TermOrder);GG:=Basis(convert(L,list), TermOrder);if G<>GG then
return(false);fi;
LMF:=Array(1..ArrayNumElems(F),[seq(LeadingMonomial(f,TermOrder), f in F)]);
for f in F do
    NM:={op(Vars)} minus {op(InvDivision(LeadingMonomial(f,TermOrder),LMF,Vars))};
    for u in NM do
    	if Normalform1(u*f,F,TermOrder,InvDivision)<>0 then
	RETURN(false);
	fi:
    od:
od:
RETURN(true);
end: