#This file presents the Maple code of F5 and F5M algorithms to detect regular sequences. 
#To apply this code, 
#1) Save this file with .mpl suffix (such as F5M.mpl), in a path (name it MyPath).
#2) Open Maple and write [> read "MyPath/F5M.mpl"; and press enter. 
#3) To apply F5, use F5RegSeq and to apply F5M, use F5MRegSeq commands. 
#The argument of each one is a list of polynomials.
#The output determines the input is regular or not, together with
#some complementary informations such as time, memory, deg, etc. 
#Don't hesitate to contact us, if there is any trouble.

with(PolynomialIdeals):with(Groebner):
######################
LM:=proc(f)
global T:
if f=0 then
    RETURN(0);
fi:
RETURN(LeadingMonomial(f,T));
end:
######################
SigOrder:=proc(s1,s2)
global Rules,Sig,Polys,Grob,T:
if s1[2]<>s2[2] then
    RETURN(evalb(s1[2]>s1[2]));
else
    RETURN(TestOrder(s1[1],s2[1],T));
fi:
end:
######################
SigSort:=proc(i,j)
global Sig;
RETURN(SigOrder(Sig[i],Sig[j]));
end:
######################
DegOrder:=proc(cp1,cp2)
RETURN(evalb(degree(cp1[1])< degree(cp2[1])));
end:
######################
CritPair:=proc(i,j,k)
global Rules,Sig,Polys,Grob,T,NF5:
local ua,ub,sa,sb,a,b,ui,uj,si,sj,t:
t:=lcm(LM(Polys[i]),LM(Polys[j])):
#ui,uj:=t/(LeadingCoefficient(Polys[i],T)*LM(Polys[i])),t/(LeadingCoefficient(Polys[j],T)*LM(Polys[j])):
ui,uj:=t/LM(Polys[i]),t/LM(Polys[j]):
si,sj:=Sig[i],Sig[j]:
a,b:=i,j:
ua,ub:=ui,uj:
sa,sb:=si,sj:
if SigOrder([ui*si[1],si[2]],[uj*sj[1],sj[2]]) then
    a,b:=j,i:
    ua,ub:=uj,ui:
    sa,sb:=sj,si:
fi:

if NormalForm(ua*sa[1],[seq(Polys[l], l in Grob[k+1])],T)<>ua*sa[1] then
    NF5:=NF5+1:
    RETURN(NULL);
elif sb[2]=k and NormalForm(ub*sb[1],[seq(Polys[l], l in Grob[k+1])],T)<>ub*sb[1] then
    NF5:=NF5+1:
    RETURN(NULL);
fi:
RETURN([t,ua,a,ub,b]);
end:
######################
AddRule:=proc(s)
global Rules,Sig,Polys,Grob,T:
Rules[s[2]]:=[op(Rules[s[2]]),s]:
end:
######################
IsRewritten:=proc(s)
#option trace;
global Rules,Sig,Polys,Grob,T,NIR:
local L,n,i:
#RETURN(false);
L:=Rules[s[2]]:
n:=nops(L):
for i from n by -1 to 1 while L[i][3]>s[3] do
    if divide(s[1],L[i][1]) then
        NIR:=NIR+1:
        RETURN(true);
    fi:
od:
RETURN(false);
end:
######################
Spol:=proc(F)
global Rules,Sig,Polys,Grob,T,NRZ:
local CP,N,i,S1,S2,r:
#option trace;
N:=NULL:
for i from 1 to nops(F) do
    S1:=Sig[F[i][3]]:
    S2:=Sig[F[i][5]]:
    if (not IsRewritten([F[i][2]*S1[1] , S1[2], S1[3]])) and (not IsRewritten([F[i][4]*S2[1] , S2[2], S2[3]])) then
        r:=simplify(1/LeadingCoefficient(Polys[F[i][3]],T)*F[i][2]*Polys[F[i][3]]-1/LeadingCoefficient(Polys[F[i][5]],T)*F[i][4]*Polys[F[i][5]]):
        if r<>0 then
            Polys(ArrayNumElems(Polys)+1):=r:
            Sig(ArrayNumElems(Sig)+1):=[F[i][2]*S1[1],S1[2],ArrayNumElems(Polys)];
            AddRule(Sig[ArrayNumElems(Sig)]);
            N:=N,ArrayNumElems(Polys):
        else 
            NRZ:=NRZ+1;
            RETURN(NRS);
            
        fi:
    fi:
od:
RETURN([N]);
end:
######################
IsRdc:=proc(r,G0,k,h)
global Rules,Sig,Polys,Grob,T:
local G,m,j,u,S:
#option trace;
G:=sort(G0,SigSort);
m:=nops(G):
for j from 1 to m do
    if divide(LM(r), LM(Polys[G[j]])) then
        u:=LM(r)/LM(Polys[G[j]]):
        S:=Sig[G[j]]:
        if NormalForm(u*S[1],[seq(Polys[l], l in Grob[k+1])],T)=u*S[1] then 
            if  u*S[1]<>Sig[h][1] and not IsRewritten([u*S[1],S[2],S[3]])then
                RETURN(G[j]);
            fi:
        fi:
    fi:
od:
RETURN(0);
end:
######################
TopRdc:=proc(r,G,k,h)
#option trace;
global Rules,Sig,Polys,Grob,T,NRZ:
local Tmp,r1,u,S,S1:
if r=0 then
    RETURN(NRS,1);
    print("The system is not a regular sequence.");
    NRZ:=NRZ+1:
    RETURN([],[]);
fi:
r1:=IsRdc(r,G,k,h):
if r1=0 then 
    RETURN([h],[]):
else
    u:=LeadingCoefficient(r,T)*LM(r)/(LeadingCoefficient(Polys[r1],T)*LM(Polys[r1])):
    S,S1:=Sig[h],Sig[r1]:
    Tmp:=simplify(r-u*Polys[r1]);
    if Tmp<>0 then
        if SigOrder([u*S1[1],S1[2],S1[3]],S) then
            Polys[h]:=Tmp;
            RETURN([],[h]);
        else
            Polys(ArrayNumElems(Polys)+1):=Tmp:
            Sig(ArrayNumElems(Sig)+1):=[u*S1[1],S1[2],ArrayNumElems(Polys)];
            AddRule(Sig[-1]):
            RETURN([],[ArrayNumElems(Polys),h]): 
        fi;
    else
        NRZ:=NRZ+1;
        RETURN(NRS,1);
        
    fi:
fi:
end:
######################
Rdc:=proc(F,k)
#option trace;
global Rules,Sig,Polys,Grob,T:
local ToDo,ToDo1,h,h1,Done,h0:
Done:=NULL:
ToDo:=sort(F,SigSort):
while nops(ToDo)<>0 do
    h:=ToDo[1]:
    ToDo:=ToDo[2..-1]:
    
    h0:=NormalForm(Polys[h],[seq(Polys[l], l in Grob[k+1])],T):
    Polys[h]:=h0;
    h1,ToDo1:=TopRdc(h0,[op(Grob[k]),Done],k,h):
    if h1=NRS then
        RETURN(NRS);
    fi;
    Done:=Done,op(h1):
    ToDo:=[op(ToDo),op(ToDo1)]:
    ToDo:=sort(ToDo,SigSort):
od:
RETURN([Done]);
end:
######################
F5:=proc(k)
#option trace;
global Rules,Sig,Polys,Grob,T,Deg:
local i,CP,d,CPd,F,Rd,r:
Grob[k]:=[op(Grob[k+1]),k]:
CP:=[seq(CritPair(k,i,k), i in Grob[k+1])]:
CP:=sort(CP,DegOrder):
while nops(CP)<>0 do
    d:=degree(CP[1][1]):Deg:=max(d,Deg):
    CPd:=NULL:
    for i from 1 to nops(CP) while degree(CP[i][1])=d do
        CPd:=CPd,CP[i]:
    od:
    CPd:=[CPd]:
    CP:=CP[i..-1]:
    F:=Spol(CPd):
    if F=NRS then 
        RETURN(NRS);
    fi;
    Rd:=Rdc(F,k):
    if Rd=NRS then 
        RETURN(NRS);
    fi;
    for r in Rd do
            CP:=[op(CP),seq(CritPair(r,i,k), i in Grob[k])]:
            Grob[k]:=[op(Grob[k]),r]:
    od:
    CP:=sort(CP,DegOrder):
od:
RETURN(Grob[k]);
end:
######################
F5RegSeq:=proc(Ideal)
#option trace;
local F,m,k,i,b1,b2,t1,t2,H:
global Rules,Sig,Polys,Grob,T,NRZ,NF5,NIR,Deg:
b1,t1:=kernelopts(bytesused,cputime):
NF5:=0:NRZ:=0:NIR:=0:
T:=tdeg(op(indets(Ideal))):
m:=nops(Ideal):
Polys:=Array(1..m,i->Ideal[i]):
Rules:=[seq([],i=1..m)]:
Sig:=Array(1..m,i->[1,i,i]):
Grob:=[seq([],i=1..m-1),[m]]:
for k from m-1 by -1 to 1 do
    print("In progress for");
    print([seq(f[i],i=k..m)]);
    Deg:=0:
    Grob[k]:=F5(k):
    if Grob[k]=NRS then
        print("Non-Regular-Sequence");
        b2,t2:=kernelopts(bytesused,cputime):
        printf("%-1s %1s %1s %1s   : %3a %3a\n",The, cpu, time, is,t2-t1,(sec)):
        printf("%-1s %1s %1s   : %3a %3a\n",The,used,memory,b2-b1,(bytes)):
        printf("%-1s %1s %1s %1s    : %3a\n",The, max, deg, is,Deg):
        printf("%-1s %1s %1s %1s     : %3a\n",The,Num, of, F5,NF5):
        printf("%-1s %1s %1s %1s    : %3a\n",Doing, via, Hilbert, Series,IsRSHilbert(Ideal)):
        print("Non-Regular-Sequence");
        RETURN();
    fi;   
od:
H:=[seq(Polys[i],i in Grob[1])]:
print("Regular-Sequence");
b2,t2:=kernelopts(bytesused,cputime):
printf("%-1s %1s %1s %1s   : %3a %3a\n",The, cpu, time, is,t2-t1,(sec)):
printf("%-1s %1s %1s   : %3a %3a\n",The,used,memory,b2-b1,(bytes)):
printf("%-1s %1s %1s %1s    : %3a\n",The, max, deg, is,Deg):
printf("%-1s %1s %1s %1s     : %3a\n",The,Num, of, F5,NF5):
printf("%-1s %1s %1s %1s    : %3a\n",Doing, via, Hilbert, Series,IsRSHilbert(Ideal)):
print("Regular-Sequence");
end:
#########################################
##### F5 Equipped by Macaulay bound #####
#########################################

F5M:=proc(k)
#option trace;
global Rules,Sig,Polys,Grob,T,Deg, MacBound:
local i,CP,d,CPd,F,Rd,r:
Grob[k]:=[op(Grob[k+1]),k]:
CP:=[seq(CritPair(k,i,k), i in Grob[k+1])]:
CP:=sort(CP,DegOrder):
while nops(CP)<>0 do
    d:=degree(CP[1][1]):
    Deg:=max(d,Deg):
    if d>MacBound then 
      if k=1 then  
        RETURN(RS);
      fi;
    fi;
    CPd:=NULL:
    for i from 1 to nops(CP) while degree(CP[i][1])=d do
        CPd:=CPd,CP[i]:
    od:
    CPd:=[CPd]:
    CP:=CP[i..-1]:
    F:=Spol(CPd):
    if F=NRS then
        RETURN(NRS);
    fi;
    Rd:=Rdc(F,k):
    if Rd=NRS then 
        RETURN(NRS);
    fi;
    for r in Rd do
            CP:=[op(CP),seq(CritPair(r,i,k), i in Grob[k])]:
            Grob[k]:=[op(Grob[k]),r]:
    od:
    CP:=sort(CP,DegOrder):
od:
RETURN(Grob[k]);
end:
######################
F5MRegSeq:=proc(Ideal)
#option trace;
local F,m,k,i,b1,b2,t1,t2,H:
global Rules,Sig,Polys,Grob,T,NRZ,NF5,NIR,Deg,Degs, MacBound:
b1,t1:=kernelopts(bytesused,cputime):
NF5:=0:NRZ:=0:NIR:=0:
Degs:=[seq(degree(p),p=Ideal)];
T:=tdeg(op(indets(Ideal))):
m:=nops(Ideal):
MacBound:=add(Degs[i],i=1..m)-m+1;
Polys:=Array(1..m,i->Ideal[i]):
Rules:=[seq([],i=1..m)]:
Sig:=Array(1..m,i->[1,i,i]):
Grob:=[seq([],i=1..m-1),[m]]:
for k from m-1 by -1 to 1 do
    print("In progress for");
    print([seq(f[i],i=k..m)]);
    Deg:=0:
    Grob[k]:=F5M(k):
    if Grob[k]=NRS then
         print("Non-Regular-Sequence");
         b2,t2:=kernelopts(bytesused,cputime):
         printf("%-1s %1s %1s %1s   : %3a %3a\n",The, cpu, time, is,t2-t1,(sec)):
         printf("%-1s %1s %1s   : %3a %3a\n",The,used,memory,b2-b1,(bytes)):
         printf("%-1s %1s %1s %1s    : %3a\n",The, max, deg, is,Deg):
         printf("%-1s %1s %1s %1s     : %3a\n",The,Num, of, F5,NF5):
         printf("%-1s %1s %1s %1s    : %3a\n",Doing, via, Hilbert, Series,IsRSHilbert(Ideal)):
         print("Non-Regular-Sequence");    
         RETURN();      
    fi;
    
    if Grob[k]=RS then
        print("Regular-Sequence");
        b2,t2:=kernelopts(bytesused,cputime):
        printf("%-1s %1s %1s %1s   : %3a %3a\n",The, cpu, time, is,t2-t1,(sec)):
        printf("%-1s %1s %1s   : %3a %3a\n",The,used,memory,b2-b1,(bytes)):
        printf("%-1s %1s %1s %1s    : %3a\n",The, max, deg, is,Deg):
        printf("%-1s %1s %1s %1s     : %3a\n",The,Num, of, F5,NF5):
        printf("%-1s %1s %1s %1s    : %3a\n",Doing, via, Hilbert, Series,IsRSHilbert(Ideal)):
        print("Regular-Sequence");     
        RETURN();
    fi;
    
od:
print("Regular-Sequence");
b2,t2:=kernelopts(bytesused,cputime):
printf("%-1s %1s %1s %1s   : %3a %3a\n",The, cpu, time, is,t2-t1,(sec)):
printf("%-1s %1s %1s   : %3a %3a\n",The,used,memory,b2-b1,(bytes)):
printf("%-1s %1s %1s %1s    : %3a\n",The, max, deg, is,Deg):
printf("%-1s %1s %1s %1s     : %3a\n",The,Num, of, F5,NF5):
printf("%-1s %1s %1s %1s    : %3a\n",Doing, via, Hilbert, Series,IsRSHilbert(Ideal)):
print("Regular-Sequence");
end:
#######################IsGrobner#############
IsGrobner:=proc(A,H,T)
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:
#######################
IsRSHilbert := proc (F) 
local f, g; 
f := HilbertSeries(F, var); 
g := mul(1-var^degree(F[i]), i = 1 .. nops(F))/(1-var)^nops(indets(F)); 
#print(series(f-g, var, 10)); 
RETURN(evalb(expand(simplify(f-g)) = 0)); 
end proc:
####################### Examples
Sequence1:=[-x_1^5*x_2^5*x_5^2+x_2*x_4^4*x_5^7-x_3*x_5^11, -x_1^3*x_4^3+x_1^2*x_2^3*x_3-x_1*x_3^5+x_2^2*x_3^3*x_4+x_2^2*x_3*x_4^3, -x_1^3+x_1*x_2^2+x_1*x_3^2]:
Sequence2:=[x_1^5*x_4^5-x_1^4*x_4^3*x_5^3-x_1^2*x_2*x_3^2*x_4^4*x_5+x_1^2*x_2*x_5^7+x_1^2*x_4^5*x_5^3-x_2^2*x_4^6*x_5^2 ,x_1*x_2^2*x_4+x_1*x_2*x_4^2+x_1*x_3^3+x_2*x_3^3-x_4^4 ,x_1^3-x_1^2*x_2+x_1^2*x_3+x_1*x_2^2-x_1*x_3^2-x_2^3+x_2^2*x_3]:
Sequence3:=[-x_1^4*x_3^3*x_4^3-x_1^3*x_4^5*x_5^2+x_1*x_2*x_3^4*x_4^2*x_5^2-x_1*x_2*x_3^3*x_4^5+x_1*x_3^8*x_4-x_2^7*x_3^3-x_2^3*x_3^2*x_4^2*x_5^3-x_2*x_3^2*x_4^3*x_5^4-x_3^3*x_4^7 ,x_1^2*x_2*x_4+x_2^2*x_3^2+x_2*x_3^3-x_2*x_3^2*x_4+x_2*x_4^3+x_3*x_4^3+x_4^4, -x_1^2*x_3-x_1*x_2*x_3+2*x_1*x_3^2-x_2^3-x_2^2*x_3+2*x_2*x_3^2]:
Sequence4:=[-x_1^9*x_4*x_5^2-x_1*x_2*x_4^9*x_5+x_2^4*x_4^6*x_5^2, -x_1^4*x_2^2-x_1*x_2*x_3^4-x_2^3*x_3^2*x_4, -x_1^3-x_1*x_2*x_3-x_3^3]:
Sequence5:=[x_1^6*x_2*x_3^2*x_4*x_5^2+x_1^4*x_3^7*x_4+x_3^4*x_4*x_5^7, -x_1^4*x_2*x_4+x_1^4*x_4^2-x_1^2*x_2^2*x_3*x_4+x_1*x_2^4*x_4-x_2^4*x_4^2+x_2^3*x_4^3+x_2*x_3^5 ,x_1^3-x_1*x_2*x_3-x_1*x_3^2-x_2^2*x_3]:
Sequence6:=[x_1^8*x_4*x_5-x_1^5*x_2^4*x_5-x_1^4*x_2^3*x_3^2*x_4+x_1^2*x_2^6*x_5^2-x_1^2*x_2^4*x_3*x_4^2*x_5-x_1^2*x_2^2*x_3^4*x_5^2 ,x_1^3*x_4+x_1^2*x_2^2+x_1^2*x_2*x_3-x_1*x_2^3+x_1*x_3^3+x_1*x_3*x_4^2-x_1*x_4^3+x_2^2*x_4^2+x_3*x_4^3, -x_1^3+2*x_1^2*x_2+x_1*x_2^2-x_2^3-x_2^2*x_3]:
Sequence7:=[x_1^6*x_3*x_5^5+x_1^5*x_2*x_4^2*x_5^4-x_1^3*x_3*x_4^2*x_5^6+x_1*x_2^7*x_5^4-x_2^2*x_3^7*x_4^2*x_5, -x_1^3*x_3^2*x_4-x_1^2*x_2*x_4^3+x_2^4*x_3*x_4+x_2^3*x_3^2*x_4+x_3^6-x_4^6, -x_1^3+x_1^2*x_2-2*x_1*x_2^2-x_1*x_3^2]:
Sequence8:=[-x_1^5*x_2*x_4^3*x_5-x_1^4*x_2^2*x_3^2*x_4*x_5+x_1^2*x_2^4*x_3^3*x_4+x_1^2*x_2^2*x_3*x_4^4*x_5+x_1*x_3^4*x_4^5-x_1*x_3*x_5^8+x_2^6*x_5^4+x_2^2*x_4^3*x_5^5, -x_1^2*x_2^2+x_1^2*x_2*x_4+x_1*x_2*x_3^2-x_1*x_2*x_4^2-x_1*x_3*x_4^2-x_1*x_4^3-x_2^3*x_3-x_2*x_3^2*x_4+x_3^4, -x_1^3+x_2*x_3^2]:
Sequence9:=[x_1^5*x_2^4*x_4+x_1^3*x_2*x_3*x_4*x_5^4-x_1^2*x_2^6*x_3*x_4+x_1^2*x_2^5*x_4^3-x_1^2*x_2^2*x_5^6+x_2*x_3^4*x_4^2*x_5^3, -x_1^3*x_3+x_1*x_2^2*x_4+x_2^3*x_3+x_2*x_3^2*x_4+x_4^4, 2*x_1^3-x_1^2*x_2+x_1^2*x_3+x_1*x_2^2-x_1*x_2*x_3+x_2^3]:
Sequence10:=[-x_1^6*x_2*x_5^5-x_1^3*x_2^6*x_3*x_4*x_5-x_2^6*x_3*x_4^2*x_5^3-x_4^2*x_5^10, -x_1^3*x_3^2*x_4+x_1^2*x_3*x_4^3+x_1*x_2^4*x_3+x_1*x_2^3*x_4^2-x_1*x_4^5 ,x_1^3-x_1^2*x_2-x_1*x_2*x_3-x_2*x_3^2]: