with(Groebner):with(PolynomialTools):with(PolynomialIdeals):
################ AutoReduced Algorithm #################
AutoReduction2:=proc(f,F,TermOrder)
#option trace;
local LMF,Vars,p,r,flag,g,u,Mult,InvDiv;
LMF:=[seq(LM(v), v in F)];
Vars:=[op(TermOrder)];
p:=f:
r:=0;
FF:=F:
while p<>0 do
      flag:=false:
      i:=1;
      while not flag and i<= nops(FF) do
          g:=FF[i]:
      	  if divide(LeadingMonomial(p,TermOrder),LeadingMonomial(g,TermOrder),'u') then
	     	flag:=true:
		p:=simplify(p-(LeadingCoefficient(p,TermOrder)/LeadingCoefficient(g,TermOrder))*u*g);
	     else i:=i+1;
	     fi:
       od:
       if not flag  then 
	  r:=simplify(r+LeadingCoefficient(p,TermOrder)*LeadingMonomial(p,TermOrder));
          p:=simplify(p-LeadingCoefficient(p,TermOrder)*LeadingMonomial(p,TermOrder));
       fi:
od:
RETURN(r);
end:

################ AutoReduce #################
AutoReduce:=proc(F,TermOrder)
local InvDiv,H,i;
#option trace;
H:=F;
for i from 1 to nops(F) do
    H[i]:=AutoReduction2(H[i],[op({op(H)} minus {H[i],0})],TermOrder);
od:
RETURN([op({op(H)} minus {0})]);
end:

##########################
#           G2v          #
##         |  |         ##
###        |  |        ###
####      _|  |_      ####
#####     \    /     #####
######     \  /     ######
#######     \/     #######
########          ########
##########################
LM:=proc(f)
global Tord;
if f<>0 then
    RETURN(LeadingMonomial(f,Tord));
fi:
RETURN(0);
end:
#################
LC:=proc(f)
global Tord;
RETURN(LeadingCoefficient(f,Tord));
end:   
#############################
#############################
tidy:=proc(p,q)
global Arxiv,Tord;
if Arxiv[p][4][2]<>Arxiv[q][4][2] then
   RETURN(evalb(Arxiv[p][4][2]>Arxiv[q][4][2]));
elif Arxiv[p][4][1]<>Arxiv[q][4][1] then
          RETURN(TestOrder(Arxiv[p][4][1],Arxiv[q][4][1],Tord));
else 
     RETURN(TestOrder(Arxiv[p][2],Arxiv[q][2],Tord));
fi:
end:
########
tidy2:=proc(p,q)
global Arxiv,Tord;
if p[2]<>q[2] then
   RETURN(evalb(p[2]>q[2]));
fi;
     RETURN(TestOrder(p[1],q[1],Tord));
end:
#######
tidy3:=proc(p,q)
global Arxiv,Tord;
if degree(Arxiv[p][4][1]*Arxiv[Arxiv[p][4][2]][3])<>degree(Arxiv[q][4][1]*Arxiv[Arxiv[q][4][2]][3]) then
   RETURN(evalb(degree(Arxiv[p][4][1]*Arxiv[Arxiv[p][4][2]][3])<degree(Arxiv[q][4][1]*Arxiv[Arxiv[q][4][2]][3])));
fi;
     RETURN(tidy(p,q));
end:
#####################
#####################
Normalform:=proc(pp)
#option trace;
local LMF,p,r,flag,g,u,v,Mult,InvDiv,LMFF,lmp,lcp;
global Vars,Arxiv,ArxivLM,Arxiv0,TT,Grob,L,NumOfSup,NumOfZero,Deg,NumOfF5,NumIs,LMG,Tord,cc1;
p:=Arxiv[pp][3]:
r:=0;
FF:=[seq(Arxiv[v][3], v in TT)];
LMFF:=[seq(Arxiv[v][2], v in TT)];
j:=0:
while p<>0 do
      lmp:=LeadingMonomial(p,Tord);
      lcp:=LeadingCoefficient(p,Tord):
      j:=j+1:
      flag:=false:
      i:=1;
      while not flag and i<= nops(FF) do
          g:=FF[i]: 
      	  if divide(lmp,LMFF[i],'u') and tidy2([lmp/LMFF[i]*Arxiv[TT[i]][4][1],Arxiv[TT[i]][4][2]],Arxiv[pp][4]) then 
		    if j=1 then 
                      if Arxiv[pp][4]=[simplify(u*Arxiv[TT[i]][4][1]),Arxiv[TT[i]][4][2]] then 
			    i:=i+1:
                      else
                              flag:=true:
                              p:=simplify(p-(lcp/LeadingCoefficient(g,Tord))*u*g);
                      fi:
                else
                      flag:=true:
                      p:=simplify(p-(lcp/LeadingCoefficient(g,Tord))*u*g);
                fi:
	   else 
                i:=i+1;
	   fi:
      od:
      if not flag  then 
            j:=j+1;
	    r:=simplify(r+lcp*lmp);
            p:=simplify(p-lcp*lmp);
      fi:
od:
RETURN(r);
end:
##############################
Invg2v:=proc(Ideal0,tord)
global Vars,Arxiv,ArxivLM,Arxiv0,TT,Grob,L,NumOfSup,NumOfZero,Deg,NumOfF5,NumIs,LMG,Tord,cc1,cc2,cc3:
local i,t1,t2,b1,b2,Ideal,temp,n,A,u1,u2:
#option trace;
t1,b1:=kernelopts(cputime,bytesused):
NumOfSup:=0:
NumOfZero:=0:
NumOfF5:=0:
NumIs:=0:
cc1:=0:
cc2:=0:
cc3:=0:
Deg:=0:
Tord:=tord:
Vars:=[op(Tord)];
L:=AutoReduce(Ideal0,Tord,InvDivision):
L:=sort(L, (a,b) -> TestOrder(b,a,Tord));
n:=nops(L):
Arxiv:=Array(1..n,i->[i,LM(L[i]),L[i],[1,i],i,{}]);
ArxivLM:=Array(1..n,i->[LM(L[i])]);
Arxiv0:=Array(1..n,i->[]);
TT:=[n];
Q:=[seq(i,i=1..n-1)];
Q:=sort(Q,tidy);
while  nops(Q)<>0 do
        p:=Q[1];
        Q:=Q[2..-1];  
	Deg:=max(Deg,degree(Arxiv[p][3]));
	h:=Normalform(p); 
	if h=0 then NumOfZero:=NumOfZero+1: fi:
    	if h=0 and Arxiv[p][4][2]>1 then 
	   Arxiv0[Arxiv[p][4][2]]:=[op(Arxiv0[Arxiv[p][4][2]]),Arxiv[p][4][1]]:
	fi:

       if h<>0 and h<>c1 then 
            lmp:=LeadingMonomial(h,Tord);
            flag:=false:
            i:=1;
            while not flag and i<= nops(TT) do
               g:=TT[i]: 
      	       if divide(lmp,Arxiv[g][2],'u') then 
                   if [u*Arxiv[g][4][1],u*Arxiv[g][4][2]]=Arxiv[p][4] then
	     	      flag:=true:
                      NumOfSup:=NumOfSup+1:
                      h:=c1:
                   else
                      i:=i+1:
                   fi:
               else
                   i:=i+1:
               fi:
             od:
        fi:
	if h<>0 and h<>c1 then
	   Arxiv[p]:=[p,LM(h),h,Arxiv(p)[4],p,{}]:
	   ArxivLM[Arxiv[p][4][2]]:=[op(ArxivLM[Arxiv[p][4][2]]),LM(h)]:
       	   for i from 1 to nops(TT) do 
       	         q:=TT[i];	 
                 A:=lcm(Arxiv[p][2],Arxiv[q][2]):
		 #if A<>Arxiv[p][2]*Arxiv[q][2] then 
		    u1:=A/Arxiv[p][2];
		    u2:=A/Arxiv[q][2]; 
		    if tidy2([u1*Arxiv[p][4][1],Arxiv[p][4][2]],[u2*Arxiv[q][4][1],Arxiv[q][4][2]]) then 
                      if not IdealMembership(u2*Arxiv[q][4][1], <seq(op(ArxivLM[u]) ,u=Arxiv[q][4][2]+1..ArrayNumElems(ArxivLM))>)  then
                           n:=n+1:
		      	   Arxiv(n):=[n,u2*Arxiv[q][2],expand(u2*Arxiv[q][3]),[u2*Arxiv[q][4][1],Arxiv[q][4][2]],Arxiv[q][5],{}]:
	              	   Q:=[op(Q),n];
                      else
                           NumOfF5:=NumOfF5+1:
                      fi:
		    else      
                      if not IdealMembership(u1*Arxiv[p][4][1], <seq(op(ArxivLM[u]) ,u=Arxiv[p][4][2]+1..ArrayNumElems(ArxivLM))>) then
                          n:=n+1:
		          Arxiv(n):=[n,u1*Arxiv[p][2],u1*Arxiv[p][3],[u1*Arxiv[p][4][1],Arxiv[p][4][2]],Arxiv[p][5],{}]:
	                  Q:=[op(Q),n];
                      else
                          NumOfF5:=NumOfF5+1:
                      fi:
		    fi:
           od:
	   TT:=[op(TT),p];
    	fi;
        Q:=sort(Q,tidy);    
	if Q<>[] then
	    Q2:=Q[1];
	    for i from 2 to nops(Q) do
		if Arxiv[Q[i-1]][4]<>Arxiv[Q[i]][4] then
		    Q2:=Q2,Q[i];
		fi:
	    od:
	    NumOfSup:=NumOfSup+(nops(Q)-nops([Q2]));
	    Q:=[Q2]:
	fi:

od:
t2,b2:=kernelopts(cputime,bytesused):
G:=[seq(Arxiv[p][3],p in TT)];
#appendto(testInc);
printf("%-1s %1s %1s %1s         : %3a %3a\n",The, cpu, time, is,t2-t1,(sec)):
printf("%-1s %1s %1s         : %5a %3a\n",The,used,memory,b2-b1,(bytes)):
printf("%-1s %1s %1s %1s       : %5a\n",Number, of, zero, reduction,NumOfZero):
printf("%-1s %1s %1s           : %5g\n",Num,of,F5criteria,NumOfF5):
printf("%-1s %1s %1s %1s       : %5a\n",Num, of, superReduction, is,NumOfSup):
printf("%-1s %1s %1s        : %3a\n",Num,of,criteria,[cc1,cc2,cc3]):
printf("%-1s %1s %1s %1s       : %5a\n",The, max, degree, is,Deg):
printf("%-1s %1s %1s             : %5a\n",Num,of,poly,nops(G)):
print("IsGrobner",IsGrobner(Ideal0,G,Tord));
#appendto(terminal):
RETURN();
end:
############################
IsGrobner:=proc(A,H,T)
#option trace;
local s,j,S,p,F,L,LL;
    F:=H;
    while member(0, F, 'p') do
        F:=subsop(p=NULL,F);
        unassign('p');
    od;
    L:=LeadingMonomial(F, T):
    LL:=LeadingMonomial(Basis(A,T),T):
    if evalb(LeadingMonomial(<op(L)>, T)<>LeadingMonomial(<op(LL)>, T)) then 
        RETURN(false);
    fi:
    for s from 1 to nops(A) do
        if evalb(Reduce(A[s],F,T)<>0) then
            RETURN(false);
        fi;
    od;
    RETURN(true);
end:

print("####################Example 1");

F:=[8*x^2-2*x*y-6*x*z+3*x+3*y^2-7*y*z+10*y+10*z^2-8*z-4,10*x^2-2*x*y+6*x*z-6*x+9*y^2-y*z-4*y-2*z^2+5*z-9 ,5*x^2+8*x*y+4*x*z+8*x+9*y^2-6*y*z+2*y-z^2-7*z+5]:
Invg2v(F,tdeg(x,y,z),J):


print("####################Example 2");

F:=[ttt+vvv-aaa,xxx+yyy+zzz+ttt-uuu-www-aaa,xxx*zzz+yyy*zzz+xxx*ttt+zzz*ttt-uuu*www-uuu*aaa-www*aaa]:
Invg2v(F,tdeg(xxx,yyy,zzz,ttt,uuu,vvv,www,aaa),J):


print("####################Example 3");

F:=[a*y^2+b*x^3+c, 2*a*y, 3*b*x^2]:
Invg2v(F,tdeg(x,y,a,b,c),J):


print("####################Example 4");

F:=[x^2+b*y^2+2*c*x*y+2*d*x+2*e*y+f, x+c*y+d, b*y+c*x+e]:
Invg2v(F,tdeg(x,y,a,b,c,d,e,f),J):


print("####################Example 5");

F:=[a*x^3*y+c*x*y^2, x^2*y+3*d*y, c*x^2+b*x*y]:
Invg2v(F,tdeg(x,y,a,b,c,d),J):


print("####################Example 6");

F:=[a*x^2*y+b*x+y^3, a*x^2*y+b*x*y, y^2+b*x^2*y+c*x*y]:
Invg2v(F,tdeg(x,y,a,b,c),J):



print("####################HAAS 3");

F := [x^6+a*y^3-y, y^6+b*x^3-x, 36*y^5*x^5-9*x^2*b*a*y^2+3*a*y^2+3*b*x^2-1]:
Invg2v(F,tdeg(x,y,a,b),J):




print("####################Liu");

F:=[y*z-y*t0-x*h+a*h, z*t0-z*x-y*h+a*h, t0*x-y*t0-z*h+a*h, x*y-z*x-t0*h+a*h]:
Invg2v(F,tdeg(z, y, x, t0, a, h),J):


print("####################Cyclic 5");

F:=[a*b*c*d*e-1, a*b*c*d + a*b*c*e + a*b*d*e + a*c*d*e +
    b*c*d*e, a*b*c + a*b*e + a*d*e + b*c*d + c*d*e, a*b + a*e + b*c +
    c*d + d*e, a + b + c + d + e]:
Invg2v(F,tdeg(a,b,c,d,e),J):


print("####################Noon");

F := [10*x1^2*x4+10*x2^2*x4+10*x3^2*x4-11*x4*h^2+10*h^3, 10*x1^2*x3+10*x2^2*x3+10*x3*x4^2-11*x3*h^2+10*h^3, 10*x1*x2^2+10*x1*x3^2+10*x1*x4^2-11*x1*h^2+10*h^3, 10*x1^2*x2+10*x2*x3^2+10*x2*x4^2-11*x2*h^2+10*h^3]:
Invg2v(F, tdeg(x1, x2, x3, x4, h), J):


print("####################Weispfenning94");

F := [y^4+x*y^2*z+x^2*h^2-2*x*y*h^2+y^2*h^2+z^2*h^2, x*y^4+y*z^4-2*x^2*y*h^2-3*h^5, -y^2*x^3+x*y*z^3+y^4*h+x*y^2*z*h-2*x*y*h^3]:
Invg2v(F, tdeg(x,y,z,h), J):



print("####################Katsura5");

F := [2*ax^2+2*ay^2+2*az^2+2*at^2+2*au^2+av^2-av, ax*ay+ay*az+2*az*at+2*at*au+2*au*av-au, 2*ax*az+2*ay*at+2*az*au+au^2+2*at*av-at, 2*ax*at+2*ay*au+2*at*au+2*az*av-az, at^2+2*ax*av+2*ay*av+2*az*av-ay, 2*ax+2*ay+2*az+2*at+2*au+av-1]:
Invg2v(F, tdeg(ax, ay, az, at, au, av), J):




print("####################Lichtblau");

F := [x-110*t^2+495*t^3-1320*t^4+2772*t^5-5082*t^6+7590*t^7-8085*t^8+5555*t^9-2189*t^10+374*t^11, y-22*t+110*t^2-330*t^3+1848*t^5-3696*t^6+3300*t^7-1650*t^8+550*t^9-88*t^10-22*t^11]:
Invg2v(F,tdeg(x,y,t), J):



print("####################Cyclic6");


    F:=[a*b*c*d*e*f -1,
    a*b*c*d*e + a*b*c*d*f + a*b*c*e*f + a*b*d*e*f + a*c*d*e*f + b*c*d*e*f,
    a*b*c*d + a*b*c*f + a*b*e*f + a*d*e*f + b*c*d*e + c*d*e*f,
    a*b*c + a*b*f + a*e*f + b*c*d + c*d*e + d*e*f,
    a*b + a*f + b*c + c*d + d*e + e*f,
    a + b + c + d + e + f]:
Invg2v(F,tdeg(a,b,c,d,e,f),J):


print("####################Katsura6");

F:=[
	    1*x1+2*x2+2*x3+2*x4+2*x5+2*x6+2*x7-1,
	    2*x4*x3+2*x5*x2+2*x6*x1+2*x7*x2-1*x6,
	    1*x3^2+2*x4*x2+2*x5*x1+2*x6*x2+2*x7*x3-1*x5,
	    2*x3*x2+2*x4*x1+2*x5*x2+2*x6*x3+2*x7*x4-1*x4,
            1*x2^2+2*x3*x1+2*x4*x2+2*x5*x3+2*x6*x4+2*x7*x5-1*x3,
	    2*x2*x1+2*x3*x2+2*x4*x3+2*x5*x4+2*x6*x5+2*x7*x6-1*x2,
            1*x1^2+2*x2^2+2*x3^2+2*x4^2+2*x5^2+2*x6^2+2*x7^2-1*x1
	   ]:
Invg2v(F, tdeg(x1, x2, x3, x4,x5,x6,x7), J):



print("####################Eco7");

F := [(x1+x1*x2+x2*x3+x3*x4+x4*x5+x5*x6)*x7-1, (x2+x1*x3+x2*x4+x3*x5+x4*x6)*x7-2, (x3+x1*x4+x2*x5+x3*x6)*x7-3, (x4+x1*x5+x2*x6)*x7-4, (x5+x1*x6)*x7-5, x6*x7-6, x1+x2+x3+x4+x5+x6+1]:
Invg2v(F, tdeg(x1, x2, x3, x4,x5,x6,x7), J):



print("####################Sturmfels and Eisenbud");

F := [ss*uu+bb*vv, tt*uu+bb*ww, tt*vv+ss*ww, ss*xx+bb*yy, tt*xx+bb*zz, tt*yy+ss*zz, vv*xx+uu*yy, ww*xx+uu*zz, ww*yy+vv*zz]:
Invg2v(F, tdeg(bb, xx, yy, zz, ss, tt, uu, vv, ww), J):