###################
###################This is the Maple implementation of our new algorithm for computing staggered linear bases.
###################
with(Groebner):
with(PolynomialIdeals):
LM := proc (f, tord)
options operator, arrow; LeadingMonomial(f, tord) end proc:

LT := proc (f, tord) options operator, arrow; LeadingMonomial(f, tord)*LeadingCoefficient(f, tord) end proc:

SOrt := proc (p1, p2, tord)
global G;
local u,v;
u := p1[1];
v:=p2[1];
RETURN(TestOrder(u, v, tord));
end proc:

Stag:=proc(FF,tord)
#option trace;
local firsttime, firstbytes,count1,count2,count3,secondtime, secondbytes,time,memory,PolList,CritPair,F,g,GG,EB,LCMList,L,TripSet,i,j,k,A,h,t;
       firsttime, firstbytes := kernelopts(cputime, bytesused);
       count1:=0: #the variable "count" counts the nubmer of reductions to zero.
       count2:=0: #the variable "count" counts the nubmer of reductions to zero.
       count3:=0: #the variable "count" pairs removed by Staggered criterion.
       F := [seq([FF[j], {seq(LM(FF[k], tord),k=1..j-1)}, LM(FF[j], tord),j], j = 1 .. nops(FF))]:
       PolList:=Array(1..nops(F),F):
       CritPair:=Array([]):
       for i from 2 to nops(F) do
                LCMList:=Array([seq(lcm(PolList[j][3],PolList[i][3]),j=1..i-1)]):
	        L:=NULL:
	     	for j from 1 to ArrayNumElems(LCMList) do
    	     	    #if not(LCMList[j]/PolList[j][3] in <op(PolList[j][2])>) then
		    if not(LCMList[j]/PolList[j][3] in <op(PolList[j][2])>) and not(LCMList[j]/PolList[i][3] in <op(PolList[i][2])>) then
		       L:=L,LCMList[j]:
		    else
		       count3:=count3+1:
		    fi:
 	        od:
	        L:=Generators(LeadingMonomial(<L>, tord));
                count2:=count2+i-1-nops(L):
                for j in L do
                    member(j,LCMList,'zz');
                    CritPair:=Array(1..ArrayNumElems(CritPair)+1,CritPair):
                    CritPair[-1]:=[j,PolList[zz],PolList[i]]:
                od:
	        PolList[i][2]:=PolList[i][2] union {seq(ss/PolList[i][3],ss=L)};
       od:
       while ArrayNumElems(CritPair)<>0 do
             CritPair:=sort(CritPair, proc (a, b) options operator, arrow; SOrt(a, b, tord) end proc):
             g:=CritPair[1]:
             CritPair:=Array(1..ArrayNumElems(CritPair)-1,i->CritPair[i+1]):
	     #PolList[g[3][4]][2]:=PolList[g[3][4]][2] union {g[1]/g[3][3]};
             h:=NormalForm(SPolynomial(g[2][1],g[3][1],tord),[seq(t[1],t=PolList)],tord);
             if h<>0 then
                       PolList:=Array(1..ArrayNumElems(PolList)+1,PolList):
                       PolList[-1]:=[h,{seq(PolList[j][3],j=1..ArrayNumElems(PolList)-1)},LM(h,tord),ArrayNumElems(PolList)]:
		       LCMList:=Array([seq(lcm(PolList[-1][3],PolList[j][3]),j=1..ArrayNumElems(PolList)-1)]):
		       L:=NULL:
	     	       for j from 1 to ArrayNumElems(LCMList) do
    	     	 	     	 #if not(LCMList[j]/PolList[j][3] in <op(PolList[j][2])>) then
				 if not(LCMList[j]/PolList[j][3] in <op(PolList[j][2])>) and not(LCMList[j]/PolList[-1][3] in <op(PolList[-1][2])>) then
		    	     	    L:=L,LCMList[j]:
		 	     	 else
				    count3:=count3+1:
		 	     	 fi:
 	     	       od:
	     	       L:=Generators(LeadingMonomial(<L>, tord));
                       count2:=count2+ArrayNumElems(PolList)-1-nops(L):
                       for j in L do
                       	   member(j,LCMList,'zz');
                    	   CritPair:=Array(1..ArrayNumElems(CritPair)+1,CritPair):
                    	   CritPair[-1]:=[j,PolList[zz],PolList[-1]]:
                       od:
		       PolList[-1][2]:=PolList[-1][2]  union {seq(ss/PolList[-1][3],ss=L)};
             else
                       count1:=count1+1:
             fi:
       od:
       Digits:=5:
       secondtime, secondbytes := kernelopts(cputime, bytesused);
       time := evalhf((secondtime-firsttime));
       memory := evalhf((evalf((secondbytes-firstbytes)/10^6)));
       GG:=[seq(PolList[i][1],i=1..ArrayNumElems(PolList))]:
       EB:=[seq(PolList[i][1..2],i=1..ArrayNumElems(PolList))];
       lprint("The required time is", time);
       lprint("The used memory", memory);
       lprint("The number of reductions to zero", count1);
       lprint("The number of removed pairs by quotient criterion", count2);
       lprint("The number of removed pairs by staggered criterion", count3);
       lprint("The number of Groebner basis elements", nops(GG));
       #lprint("Elimination basis", EB);
       #lprint("Is Grobner?", IsGrobner(FF,GG,tord));
       #lprint("Is Staggrred?", staggered(EB,tord));
       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(Basis(LeadingMonomial(F, T),T),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:
##############################
staggered:=proc(F,tord)
#option trace:
local count,n,i,j,Lcm,m1,m2,T,t,flag,l;
n:=nops(F):
count:=0:
for i from 1 to n do
    for j from i+1 to n do
        count:=count+1:
        Lcm:=lcm(LM(F[i][1],tord),LM(F[j][1],tord)):
        m1:=(Lcm)/(LM(F[i][1],tord)):
        m2:=(Lcm)/(LM(F[j][1],tord)):
        if not (m1 in LeadingMonomial(<op(F[i][2])>, tord)) and  not (m2 in LeadingMonomial(<op(F[j][2])>, tord)) then
           print("B is not an independent set and not Staggered basis",count);
           RETURN();
        fi:
    od:
od:
#print("B is independent set",count);
for i from 1 to n do
    T:=F[i][2]:
    while op(T)<> NULL do
    t:=T[1]:
    T:=T minus {t}:
    flag:=false:
    l:=1:
    while flag=false and l<= n do
          if divide(t*LM(F[i][1],tord), LM(F[l][1],tord), 'q') then
             if (q in LeadingMonomial(<op(F[l][2])>, tord)) then
                flag:=true:
             fi:
          fi:
          l:=l+1:
    od:
    if flag=false then
       print("B is not a generating set and not Staggered");
       RETURN(0):
    fi:
od:
od:
#print("B is a generating set");
RETURN(true):
end proc:
###################
###################This is the Maple implementation of Moller-Mora-Traverso to compute staggered linear bases.
###################
with(Groebner):
with(PolynomialIdeals):
LM := proc (f, tord)
options operator, arrow; LeadingMonomial(f, tord) end proc:

LT := proc (f, tord) options operator, arrow; LeadingMonomial(f, tord)*LeadingCoefficient(f, tord) end proc:

SOrt := proc (p1, p2, tord)
global G;
local u,v;
u := p1[1];
v:=p2[1];
RETURN(TestOrder(u, v, tord));
end proc:

Stag:=proc(FF,tord)
#option trace;
local firsttime, firstbytes,count1,count2,count3,secondtime, secondbytes,time,memory,PolList,CritPair,F,g,GG,EB,LCMList,L,TripSet,i,j,k,A,h,t;
       firsttime, firstbytes := kernelopts(cputime, bytesused);
       count1:=0: #the variable "count" counts the nubmer of reductions to zero.
       count2:=0: #the variable "count" counts the nubmer of reductions to zero.
       F := [seq([FF[j], {seq(LM(FF[k], tord),k=1..j-1)}, LM(FF[j], tord),j], j = 1 .. nops(FF))]:
       PolList:=Array(1..nops(F),F):
       CritPair:=[seq(seq([lcm(PolList[i][3],PolList[j][3]),PolList[i],PolList[j]],i=1..j-1),j=1..nops(FF))]:
       CritPair:=Array(1..nops(CritPair),CritPair):
       while ArrayNumElems(CritPair)<>0 do
             CritPair:=sort(CritPair, proc (a, b) options operator, arrow; SOrt(a, b, tord) end proc):
             g:=CritPair[1]:
             CritPair:=Array(1..ArrayNumElems(CritPair)-1,i->CritPair[i+1]):
                if not(g[1]/g[2][3] in <op(g[2][2])>) and not(g[1]/g[3][3] in <op(g[3][2])>) then
                    PolList[g[3][4]][2]:=PolList[g[3][4]][2] union {g[1]/g[3][3]};
                    h:=NormalForm(SPolynomial(g[2][1],g[3][1],tord),[seq(t[1],t=PolList)],tord);
                    if h<>0 then
                       PolList:=Array(1..ArrayNumElems(PolList)+1,PolList):
                       PolList[-1]:=[h,{seq(PolList[i][3],i=1..ArrayNumElems(PolList)-1)},LM(h,tord),ArrayNumElems(PolList)]:
                       CritPair:=Array(1..ArrayNumElems(CritPair)+ArrayNumElems(PolList)-1,CritPair):
                       for i from 1 to ArrayNumElems(PolList)-1 do
                                CritPair[ArrayNumElems(CritPair)-ArrayNumElems(PolList)+i+1]:=[lcm(PolList[i][3],PolList[-1][3]),PolList[i],PolList[-1]]:
                        od:
                    else
                       count1:=count1+1:
                    fi:
                else
                    count2:=count2+1:
                fi:
       od:
       secondtime, secondbytes := kernelopts(cputime, bytesused);
       time := evalf(secondtime-firsttime);
       memory := evalf((secondbytes-firstbytes)/10^6);
       GG:=[seq(PolList[i][1],i=1..ArrayNumElems(PolList))]:
       EB:=[seq(PolList[i][1..2],i=1..ArrayNumElems(PolList))];
       lprint("The required time is", time);
       lprint("The used memory", memory);
       lprint("The number of reductions to zero", count1);
       lprint("The number of removed pair is", count2);
       lprint("The number of Groebner basis elements", nops(GG));
       #lprint("Elimination basis", EB);
       #lprint("IsGrobner?", IsGrobner(FF,GG,tord));
       #lprint("IsStaggered?", staggered(EB,tord));
       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(Basis(LeadingMonomial(F, T),T),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:
##############################
staggered:=proc(F,tord)
#option trace:
local count,n,i,j,Lcm,m1,m2,T,t,flag,l;
n:=nops(F):
count:=0:
for i from 1 to n do
    for j from i+1 to n do
        count:=count+1:
        Lcm:=lcm(LM(F[i][1],tord),LM(F[j][1],tord)):
        m1:=(Lcm)/(LM(F[i][1],tord)):
        m2:=(Lcm)/(LM(F[j][1],tord)):
        if not (m1 in LeadingMonomial(<op(F[i][2])>, tord)) and  not (m2 in LeadingMonomial(<op(F[j][2])>, tord)) then
           print("B is not an independent set and not Staggered basis",count);
           RETURN();
        fi:
    od:
od:
#print("B is independent set",count);
for i from 1 to n do
    T:=F[i][2]:
    while op(T)<> NULL do
    t:=T[1]:
    T:=T minus {t}:
    flag:=false:
    l:=1:
    while flag=false and l<= n do
          if divide(t*LM(F[i][1],tord), LM(F[l][1],tord), 'q') then
             if (q in LeadingMonomial(<op(F[l][2])>, tord)) then
                flag:=true:
             fi:
          fi:
          l:=l+1:
    od:
    if flag=false then
       print("B is not a generating set and not Staggered");
       RETURN(0):
    fi:
od:
od:
#print("B is a generating set");
RETURN(true):
end proc:
###################
###################This is the Maple implementation of Hashemi-Javanbakht to compute staggered linear bases.
###################
with(Groebner):
with(PolynomialIdeals):
lm := proc (f, tord) options operator, arrow; LeadingMonomial(f, tord) end proc:
lc := proc (f, tord) options operator, arrow; LeadingCoefficient(f, tord) end proc:
lt := proc (f, tord) options operator, arrow; LeadingTerm(f, tord)[1]*LeadingTerm(f,tord)[2] end proc:

LM := proc(u)
local i;
for i to nops(u) do
   if u[i] <> 0 then
   RETURN([u[i], i])
   end if:
end do:
RETURN(0)
end proc:

#`#########################`
pot := proc (a, b, tord)
i := a[2];
j := b[2];
t1 := a[1];
t2 := b[1];
if i < j then
   RETURN(true):
elif j < i then
   RETURN(false):
elif TestOrder(t1, t2, tord) then
   RETURN(true):
else
   RETURN(false):
end if:
end proc:
#`####################`
Pot := proc (a, b, tord)
#option trace:
 global Fsub:
 t1 := a[1];
 i := a[2];
 t2 := b[1];
 j := b[2];
 if Fsub[j][3][2] < Fsub[i][3][2] then
      RETURN(false):
 elif Fsub[i][3][2] < Fsub[j][3][2] then
      RETURN(true):
 elif TestOrder(t1*Fsub[i][3][1], t2*Fsub[j][3][1], tord) then
      RETURN(true):
 else
      RETURN(false):
 end if:
 end proc:
 #`########## Cover function ############`
 Cover := proc (p1, p2, tord)
 local u1, u2, v1, v2;
 u1 := p1[1];
 u2 := p2[1];
 if not evalb(u1[2] = u2[2]) then
      RETURN(false):
 elif divide(u1[1], u2[1],'q') then
      t := q*p2[2];
      if TestOrder(t, p1[2],tord) and t <> p1[2] then
               RETURN(true):
      end if:
 end if;
 RETURN(false):
 end proc:
#`########## Regtopred function #########`
 regtopred := proc (p1,tord)
 #option trace:
 local Flag;
 global Fsub:
 Flag := true;
 G:=Fsub:
 n := op(2, op(2, G));
 #Q:=[seq(0,i=1..n)]:
Q:=Array(1..n,i->0):
r:=[0,0]:
 p := [p1[1], p1[2]];
 while p[2]<> 0 do
      i := 1;
      flag := false;
      u1, v1:= p[1], p[2];
      while flag = false and i <= n do
      u2, v2 := G[i][3], G[i][1];
      if divide(lm(v1, tord), G[i][5], 'q') then
          t := q;
          c := lc(v1, tord)/G[i][6];
          if pot([t*u2[1], u2[2]], u1, tord) then
               if [t*u2[1], u2[2]]<>u1 or c <> 1 then
                    p := [u1, v1-expand(c*t*v2)];
                    Q[i]:=Q[i]+c*t:
                    flag := true;
                    if [t*u2[1], u2[2]] = u1 then
                       p[2] := p[2]/(1-c):
                    end if:
               end if:
          end if:
      end if;
      i := i+1:
 end do;
 if flag = false then
         r:=[u1,simplify(r[2]+lt(p[2],tord))]:
         p[2]:=simplify(p[2]-lt(p[2],tord)):
 end if:
 end do;
 RETURN(r,Q):
 end proc:
 #########Checkpair##############
 checkpair:=proc(m,A)
 #option trace:
 flag:=true:  i:=1:
 while flag and i<= nops(A) do
       if divide(m,A[i],'q')=true then
          RETURN(false):
       fi:
       i:=i+1:
 od:
 RETURN(true);
 end proc:
 
######### Staggered linear basis algorithm based on the incrimental process######
 Stag:=proc(F1,tord)
 #option trace:
 global Fsub:
 tim1,b1:=kernelopts(cputime,bytesused):
 NumOfTop:=0:
 NumOfZero:=0:
 NumOfCov:=0:
 NumOfStaggeredRem:=0:
 csp:=0:
 F:=[seq([F1[i],{},[1,i],i,lm(F1[i],tord),lc(F1[i],tord)],i=1..nops(F1))]:
 Fsub:=[F[1]]:
 Fsub:=Array(1..1,Fsub):
 F:=F[2..-1]:
 H:=[]:
 JP:=Array([]):
 while op(F)<>NULL do
     n:=op(2, op(2, Fsub)):
     f:=F[1]:
     F:=F[2..-1]:
     f[2]:={seq(g[5],g in Fsub)}:
     f[4]:=n+1:
     LMG:=[seq(Fsub[i][5],i=1..n)]:
     #print("f",f):
     #print("Fsub",Fsub):
     for p in Fsub do
         Lcm:=lcm(p[5],f[5]):
         m1:=(Lcm)/(p[5]):
         m2:=(Lcm)/(f[5]):
         m1u1:=[m1*p[3][1],p[3][2]]:
         m2u2:=[m2*f[3][1],f[3][2]]:
         if m1u1<>m2u2  then
            Flag:= pot(m1u1,m2u2,tord):
            Num:=op(2,op(2,JP)):
            
            if Flag and checkpair(m2,f[2]) then
	       JP:=Array(1..Num+1,JP):
               JP[Num+1]:=[m2,n+1]:
            elif Flag=false and checkpair(m1,p[2]) then
	       JP:=Array(1..Num+1,JP):
               JP[Num+1]:=[m1,p[4]]:
            #elif (Flag and checkpair(m2,f[2])=false) or (Flag=false and checkpair(m1,p[2])=false) then
               #NumOfCov:=NumOfCov+1:
            else
                #prnit("Removed by staggered");
                NumOfStaggeredRem:=NumOfStaggeredRem+1:
            fi:
            elif checkpair(m2,f[2])  then
               csp:=csp+1:
               p[2]:=p[2] union {m1}:
         fi:
     od:
     Fsub:=Array(1..n+1,Fsub):
     Fsub[n+1]:=f:
     n:=n+1:
     ##In this loop we try to obtain a Groebner basis for the ideal <Fsub>
     while ArrayNumElems(JP)<>0 do
        JP:=sort(JP,(a,b)->Pot(a,b,tord)):
        t,i:=JP[1][1],JP[1][2]:
        #print("####################################################################");
        #print("t,i", t,i);
        p:=Fsub[i]:
        #print("p",p);
	JP:=Array(1..ArrayNumElems(JP)-1,i->JP[i+1]):
	jp:=[[t*p[3][1],p[3][2]],expand(t*p[1]),t*p[5]]:
        flag:=false:
        if checkpair(t,p[2])=false  then
          flag:=true:
          NumOfStaggeredRem:=NumOfStaggeredRem+1:
          #print("Remove by staggered");
          ###### elif checkpair`(t,LMG)=false then
          ###### flag:=true:
          #####NumOfCov:=NumOfCov+1:
          #print("Remove by forbidden LM");
        fi:
        if flag=false then
           r,Q:=regtopred(jp,tord):
           #print("r",r):
           #print("Q",Q):
           NumOfTop:=NumOfTop+1:
           if r[2]=0 then
               NumOfZero:=NumOfZero+1:
               Fsub[i][2]:=Fsub[i][2] union {t}:
           else
               qoutiontt:={}:
               for su in Fsub[i][2] do
                   l:=numer(simplify(su/t)):
                   qoutiontt:={op(qoutiontt),l}:
               end do:
               r:=[r[2],qoutiontt,r[1],n+1,lm(r[2],tord),lc(r[2],tord)]:
               Fsub[i][2]:=Fsub[i][2] union {t}:
               #print("r",r):
               JPnew:=Array([]):
               for p in Fsub do
                  Lcm:=lcm(p[5],r[5]):
                  m1:=Lcm/p[5]:
                  m2:=Lcm/r[5]:
                  t1u1:=[m1*p[3][1],p[3][2]]:
                  t2u2:=[m2*r[3][1],r[3][2]]:
                  if t1u1<>t2u2  then
                      Flag:= pot(t1u1,t2u2,tord):
                      if Flag and checkpair(m2,r[2]) then
		          endd:=op(2,op(2,JPnew)):
                          JPnew:=Array(1..endd+1,JPnew):
                          JPnew[endd+1]:=[m2,n+1]:
                      elif Flag=false and checkpair(m1,p[2]) then
		          endd:=op(2,op(2,JPnew)):
			  JPnew:=Array(1..endd+1,JPnew):
                          JPnew[endd+1]:=[m1,p[4]]:
                     #elif(Flag and checkpair(m2,r[2])=false) or (Flag=false and checkpair(m1,p[4])=false) then                                    #`
                     #NumOfCov:=NumOfCov+1:
                      else
                      NumOfStaggeredRem:=NumOfStaggeredRem+1:
                     fi:
                  elif checkpair(m2,r[2])  then
                     csp:=csp+1:
                     p[2]:=p[2] union {m1}:
                  fi:
               od:
               n:=n+1:
               Fsub:=Array(1..n,Fsub):
               Fsub[n]:=r:
               #print("Fsub",Fsub):
               #print("JPnew",JPnew,JP):
               endJp:=ArrayNumElems(JP);
               nnn:=ArrayNumElems(JPnew):
               JP:=Array(1..endJp+nnn,JP):
               for k from endJp+1 to endJp+nnn do
                  JP[k]:=JPnew[k-endJp];
               od:
               #JP:=[op(JP),op(JPnew)]:
               JP:=sort(JP,(a,b)->Pot(a,b,tord)):
               #print("JP",JP):
         fi:
       fi:
     od:
  #print("JP is NUll"):
end do:
tim2,b2:=kernelopts(cputime,bytesused):
L:=[seq(f[1],f in Fsub)]:
EB:=[seq(Fsub[i][1..2],i=1..ArrayNumElems(Fsub))];
       lprint("The required time is", tim2-tim1);
       lprint("The used memory", evalf((b2-b1)/10^6));
       lprint("The number of reductions to zero", NumOfZero);
       lprint("The number of removed pairs by cover criterion", NumOfCov);
       lprint("The number of removed pairs by csp criterion", csp+NumOfStaggeredRem);
       lprint("The number of Groebner basis elements", nops(L));
       #lprint("Is Grobner?", IsGrobner(F1,L,tord));
       #lprint("Is Staggrred?", staggered(EB,tord));
       RETURN();
end proc:
##############################
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:
##############################
staggered:=proc(F,tord)
#option trace:
local count,n,i,j,Lcm,m1,m2,T,t,flag,l;
n:=nops(F):
count:=0:
for i from 1 to n do
    for j from i+1 to n do
    	count:=count+1:
	Lcm:=lcm(lm(F[i][1],tord),lm(F[j][1],tord)):
	m1:=(Lcm)/(lm(F[i][1],tord)):
	m2:=(Lcm)/(lm(F[j][1],tord)):
	if not (m1 in LeadingMonomial(<op(F[i][2])>, tord)) and  not (m2 in LeadingMonomial(<op(F[j][2])>, tord)) then
	   print("B is not an independent set and not Staggered basis",count);
	   RETURN();
	fi:
    od:
od:
#print("B is independent set",count);
for i from 1 to n do
    T:=F[i][2]:
    while op(T)<> NULL do
    t:=T[1]:
    T:=T minus {t}:
    flag:=false:
    l:=1:
    while flag=false and l<= n do
    	  if divide(t*lm(F[i][1],tord), lm(F[l][1],tord), 'q') then
	     if (q in LeadingMonomial(<op(F[l][2])>, tord)) then
	     	flag:=true:
	     fi:
	  fi:
	  l:=l+1:
    od:
    if flag=false then
       print("B is not a generating set and not Staggered");
       RETURN(0):
    fi:
od:
od:
#print("B is a generating set");
RETURN(true):
end proc:
###################
################### Examples
###################
#######Cyclic5
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];
Stag(F, tdeg(a, b, c, d, e)):

#######Sturmfels and Eisenbud
F :=[bb*vv+ss*uu, bb*ww+tt*uu, ss*ww+tt*vv, bb*yy+ss*xx, bb*zzz+tt*xx, ss*zzz+tt*yy, uu*yy+vv*xx, uu*zzz+ww*xx, vv*zzz+ww*yy];
Stag(F, tdeg(bb, xx, yy, zzz, ss, tt, uu, vv, ww));

#######################Eco8 
F:=[(x1 + x1*x2 + x2*x3 + x3*x4 + x4*x5 + x5*x6 + x6*x7)*x8 - 1,
         (x2 + x1*x3 + x2*x4 + x3*x5 + x4*x6 + x5*x7)*x8 - 2,
                 (x3 + x1*x4 + x2*x5 + x3*x6 + x4*x7)*x8 - 3,
                         (x4 + x1*x5 + x2*x6 + x3*x7)*x8 - 4,
                                 (x5 + x1*x6 + x2*x7)*x8 - 5,
                                         (x6 + x1*x7)*x8 - 6,
                                                   x7*x8 - 7,
                        x1 + x2 + x3 + x4 + x5 + x6 + x7 + 1];
			Stag(F, tdeg(x1, x2, x3, x4, x5,x6,x7,x8));
##############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];
 Stag(F, tdeg(x1, x2, x3, x4, x5,x6,x7));

#######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]:
Stag(F, tdeg(a, b, c, d, e,f)):

###################Cohn2
F:=[x^3*y^2+4*x^2*y^2*z-x^2*y*z^2+288*x^2*y^2+207*x^2*y*z+1152*x*y^2*z+
   156*x*y*z^2+x*z^3-3456*x^2*y+20736*x*y^2+19008*x*y*z+82944*y^2*z
   +432*x*z^2- 497664*x*y+62208*x*z+2985984*x,
 y^3*t^3+4*y^3*t^2-y^2*z*t^2+4*y^2*t^3-48*y^2*t^2-5*y*z*t^2
   +108*y*z*t+ z^2*t+144*z*t-1728*z,
 -x^2*z^2*t+4*x*z^2*t^2+z^3*t^2+x^3*z+156*x^2*z*t+207*x*z^2*t+
   1152*x*z*t^2+288*z^2*t^2+432*x^2*z+19008*x*z*t-3456*z^2*t+82944*x*t^2+
   20736*z*t^2+62208*x*z-497664*z*t+2985984*z,
 y^3*t^3-x*y^2*t^2+4*y^3*t^2+4*y^2*t^3-5*x*y^2*t-48*y^2*t^2+x^2*y+
   108*x*y*t+144*x*y-1728*x];
  Stag(F, tdeg(x,y,z,t)): 

#############Huneke
F := [s^15, t0^15, u^15, u^5-s^3*t0*x+s^2*t0^2*x+s^2*t0^2*y-s*t0^3*y]:
Stag(F, tdeg(x, y, s, t0, u, z)):

##############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]:
Stag(F, tdeg(x,y,t)):

##############Chemkin homogenized
F := [ z2+z3+z4+z5, 33461*c3-94642*h, c2+y2+y3+y4+y5, 9*c2-8*h, 36926*c1-208885*h, -3*c3*y2+3*x3*h+3*x4*h+8*h^2, -c1*y2+9*y2^2+z2*h, z5^2+y5^2-c2*h, 3*z4*z5+3*y4*y5+x4*h-h^2, z4^2+y4^2+x4^2-h^2, 3*z3*z4+3*y3*y4+3*x3*x4-h^2, z3^2+y3^2+x3^2-h^2, -3*c3*y2*x3+3*z2*z3*h+3*y2*y3*h+3*x3*h^2-h^3]:
Stag(F,tdeg(c1, c2, c3, z2, z3, z4, z5, y2, y3, y4, y5, x3, x4, h)):

############HAAS3
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]:
Stag(F,tdeg(x,y,a,b)):
#############################