#############################################################################################################################################
#############################################################################################################################################
#############################################################################################################################################
########################################## This file contains the implementation of three algorithms: #######################################
########################################## 1. GVWDisPGB algorithm (Combination of GVW (Gao...) and DisPGB (Montes) algorithm) ############### 
########################################## 2. GVW algorithm of Gao et al. 2016 ##############################################################
########################################## 3. PGBMain (KSW) algorithm of Kapur et al. 2010 ##################################################
#############################################################################################################################################
#############################################################################################################################################
#############################################################################################################################################
############################################################################################################################
############################################################################################################################
                                          ##########################
                                          #        GVWDisPGB       #
                                          #        algorithm       #
                                          ##         |  |         ##
                                          ###        |  |        ###
                                          ####      _|  |_      ####
                                          #####     \    /     #####
                                          ######     \  /     ######
                                          #######     \/     #######
                                          ########          ########
                                          ##########################

with(Groebner):
with(Algebraic):
with(PolynomialIdeals):
############################Select Algorithm###################### 
sff := proc (x) 
evalb(type(x, 'constant')): 
end proc:
################
lt := proc (f, T) 
RETURN(LeadingTerm(f, T)[1]*LeadingTerm(f, T)[2]): 
end proc: 
#############Normalize Algorithm################
nrm := proc (F) 
local KK, h, i, NM; 
#option trace; 
if F = [] then 
RETURN(F); 
end if; 
KK := F; 
if type(denom(KK[1]), monomial) then 
h := LeadingMonomial(denom(KK[1]), plex(op(indets(denom(KK[1]))))); 
else 
h := denom(KK[1]); 
end if; 
for i from 2 to nops(KK) do 
if type(denom(KK[i]), monomial) then 
h := lcm(h, LeadingMonomial(denom(KK[i]), plex(op(indets(denom(KK[i])))))); 
else 
h := lcm(h, denom(KK[i])); 
end if; 
end do; 
NM := simplify(expand(h*KK)); 
RETURN(NM); 
end proc:
###############convert a list of polynomial to the set "FACVAR"##########
fac := proc (L, T) 
local A, AA, AAA, N, P, p, B, C, i; 
#option trace;
A := L; 
AA := [seq(factor(i), `in`(i, A))]; 
AAA := NULL; 
B := []; 
for i to nops(AA) do 
    p := AA[i]; 
    if irreduc(p) = true then 
        AAA := AAA, p: 
    else 
        AAA := AAA, op(convert(p, list)): 
    end if 
end do; 
P := NULL; 
for i to nops([AAA]) do 
    if `not`(type(AAA[i], 'constant')) then 
        P := P, [AAA][i]: 
    end if 
end do; 
N := NULL; 
for i to nops([P]) do 
    if irreduc([P][i]) = false and type([P][i], 'constant') <> true then 
        N := N, op([P][i])[1]: 
    else 
        N := N, [P][i]: 
    end if 
end do; 
RETURN({N}): 
end proc:

################## New Condition #################
new := proc (N, W, f, R, T) 
local ff, test, NN, WW, cd, KK, ww, i, cc, ccc, k, M, www,fff; 
#option trace:
ff := simplify(expand(f)); 
test := true; 
NN := N; 
KK := []; 
while test and NN <> [1] and ff <> 0 do 
    while type(LeadingCoefficient(ff, T), 'constant') = true and ff <> 0 do 
        KK := [op(KK), lt(ff, T)]; 
        ff := simplify(expand(ff-simplify(expand(lt(ff, T))))) 
    end do; 
    if RadicalMembership(LeadingCoefficient(simplify(ff), T), simplify(`<,>`(NN))) = true then 
        NN := [op({op(NN), LeadingCoefficient(ff, T)})]; 
        ff := simplify(ff-simplify(expand(lt(ff, T)))): 
    else 
        test := false:
    end if :
end do; 
NN := Basis(NN, R); 
WW := [op({seq(NormalForm(W[i], NN, R), i = 1 .. nops(W))})]; 
ff := op(nrm([NormalForm(ff, NN, R)])); 
if WW <> [] then 
    #`WW&Assign;nrm`(WW);
     WW:=nrm(WW):
    ww := NULL; 
    for i to nops(WW) do 
        if not type(WW[i], 'constant') = true then 
            ww := ww, WW[i] 
        end if 
    end do; 
    WW := [ww] 
end if; 
fff := ff; 
while type(LeadingCoefficient(fff, T), 'constant') = true and fff <> 0 do 
    fff := simplify(expand(fff-simplify(expand(lt(fff, T))))) 
end do; 
cc := fac([LeadingCoefficient(fff, T)]); 
ccc := {seq(op(nrm([k/LeadingCoefficient(k, R)])), `in`(k, cc))};
cd := `minus`({seq(expand(kk), `in`(kk, ccc))}, {-op(WW), op(WW)}); 
k := add(KK[j], j = 1 .. nops(KK)); 
M := [op(fac(WW, R))]; 
WW := [op({seq(ii/LeadingCoefficient(ii, R), `in`(ii, M))})]; 
www := [seq(nrm([ee]), `in`(ee, WW))]; 
if www <> [] then 
    www := [op(www)]; 
www := [seq(op(jj), `in`(jj, www))] 
end if; 
if nops(cd)<>0 then
     RETURN([{op([nrm([op(cd minus select(sff,cd))][1])])}, ff+k, NN, www]); 
else
     RETURN([{}, ff+k, NN, www]);
fi:
end proc:
##################
sf := proc (x) 
evalb(x <> 0): 
end proc: 

##############################
##                          ##
##############################
LM:=proc(f)
global tord;
if f<>0 then
    RETURN(LeadingMonomial(op(nrm([f])),tord));
fi:
RETURN(0);
end:
LC:=proc(f)
global tord;
RETURN(LeadingCoefficient(op(nrm([f])),tord));
end:
#############################
refine:=proc(LL,Jsig)
local u,f:
u:=Jsig:
f:=proc(Q)
global u:
if not divide(LM(Q[1]),u) then
    RETURN(true);
fi:
RETURN(false);
end:
RETURN(select(f,LL)):
end:
#############################
JPair2:=proc(P,Q,LMH)
#option trace;
global tord,Polys,Deg,NumOfFirstBuch,m,DDD;
local t,t1,t2,f,L,T:
####
f:=proc(r)
if r<>0 then 
    RETURN(true);
fi:
RETURN(false);
end:
####
T:=tord;
L:=select(f,LMH):
t:=lcm(P[3],Q[3]):
DDD:=max(DDD,degree(t));
t1,t2:=simplify(t/P[3]) , simplify(t/Q[3]):
if TestOrder(t1*LM(P[1]) , t2*LM(Q[1]),T) and NormalForm(t2*LM(Q[1]),L,T)<>0  then
    RETURN(expand([LeadingCoefficient(Polys[P[2]],T)*t2*Q[1],Q[2],t2*Q[3],t2*Q[4],P[2]]));
elif TestOrder(t2*LM(Q[1]) , t1*LM(P[1]),T) and NormalForm(t1*LM(P[1]),L,T)<>0 then 
    RETURN(expand([LeadingCoefficient(Polys[Q[2]],T)*t1*P[1],P[2],t1*P[3],t1*P[4],Q[2]]));
fi:
RETURN(NULL);
end:
#############################
tidy:=proc(P,Q)
global tord;
local T:
T:=tord;
if TestOrder(LM(P[1]),LM(Q[1]),T) then
    RETURN(true);
fi:
RETURN(false);
end:
#############################
RegRed:=proc(P)
#option trace:
local u1,v1,i,flag,Q,u2,v2,t,c,cmpt,T:
global tord,r,Grob,Polys,NumOfSup,R,GG,tp:
T:=tord;
u1,v1:=P[1],expand(P[4]*Polys[P[2]]):
i:=1:
cmpt:=0;
while i<=nops(R) and v1<>0 do
    flag:=false:
    Q:=R[i]:
    u2,v2:=Q[1],expand(Q[4]*Polys[Q[2]]):
    t:=LM(v1)/LM(v2):
    if divide(LM(v1),LM(v2)) and TestOrder(LM(expand(t*u2)) , LM(u1),T) and (P[2]<>Q[2] or cmpt=1) then# and Q[1]<>t*P[1] then
        cmpt:=1;
        c:=LeadingCoefficient(op(nrm([v1])),T)/LeadingCoefficient(op(nrm([v2])),T):
        if LM(u1)=t*LM(u2) then 
           RETURN(1); 
        fi;
        u1:=op(nrm([simplify(u1-c*t*u2)])): 
        v1:=SPolynomial(v1, v2, T);
        if LM(expand(t*u2)) = LM(u1) and c<>1 then
            u1:=op(nrm([u1/(1-c)])):
            v1:=op(nrm([v1/(1-c)])):
        fi:
        u1:=op(nrm([NormalForm(u1,Grob,T)])):
        v1:=op(nrm([NormalForm(v1,Grob,T)])):
        i:=1:
        flag:=true:
    else
        i:=i+1:
    fi:
od:
RETURN([op(nrm([u1])),op(nrm([v1]))]);
end:
############################
gvw:=proc(fff,GGG,N,W,RR,T,CPP,RRR)
#option trace:
local n,U,V,LMH,v,u,JP,NF,UV,Pm,g,flag,nr,m,P,k,Y,c,S,test,CP,f,fun,i,funi,G;
global tp,Grob,L,Polys,R,NumOfZero,NumOfF5,H,tord,null,notnull,fflag,iii,DDD,GG,rdz:
f:=fff[2];
CP:=CPP:
GG:=NULL:
tp:=RR;
R:=NULL:
fun:=proc (p, i) if i in {p[2],p[-1]} then RETURN(true) else RETURN(false) end if end proc;
for i from 1 to nops(GGG) do
    g:=NormalForm(GGG[i],N,RR):
    if g=0 then 
    CP:=remove(fun,CP,i);
    funi:=proc(p,j) if p<j then RETURN(p); else RETURN(p-1); fi;end:  
    CP:={seq([p[1],funi(p[2],i),p[3],p[4],funi(p[5],i)], p in CP)}:
    else
       GG:=GG,g;
       R:=R,[RRR[i][1],RRR[i][2],LM(g),RRR[i][4]];
    fi;
od;
G:=[GG];if N=[1] then RETURN(true,G,N,W,0,{},{}); ;fi;
R:=[R];
tord:=T;
null, notnull:=N,W:
Polys:=G:
m:=nops(Polys);
LMH:=[seq(LM(Grob[j]), j=1..nops(Grob))]:
g:=expand(NormalForm(f,N,RR)):
v:=op(nrm([NormalForm(g,Grob,T)])):
if v=0 and CP={} then 
    RETURN(true,G,N,W,0,{},{}):
elif v=0 then
    JP:=CP:
else
    
    u:=fff[1];
    Y := new(null, notnull, op(nrm([v])), RR, T);
                    c[dec] := Y[1];
                    S:=Y[2]: 
                    null := Y[3]; 
                    notnull := Y[4]; 
                    if S=0 then RETURN(true,G,N,W,0,{},{}): fi:
                    if c[dec] = {} then 
                        if S <> 0 then 
                             v:=simplify(op(nrm([S])));
                             nr:=nops(R):
                             R:=[op(R),[u,nr+1,LM(v),1]]:
                             Polys:=[op(Polys),v]:
                             JP:=[op({op(CP),seq(JPair2(R[-1],R[j],LMH),j=1..nr)})]:
                             JP:=sort(JP,tidy);
                         end if: 
                    elif S <> 0 then 
                        RETURN(false,Polys,null,notnull,[u,S],CP,R);
                    end if: 
fi:
while nops(JP)<>0 do 
        P:=JP[1]:
        JP:=JP[2..-1]:
        NF:=RegRed(P):
        if NF<>1 then
            u:=NF[1]:
            v:=NF[2]:
            if v=0 then
                #H:=[op(H),u]:
                LMH:=[op(LMH),LM(u)]:
                JP:=refine(JP,LM(u)):
                rdz:=rdz+1;
            else
                    v:=NormalForm(v,N,RR);
                    Y := new(null, notnull, op(nrm([v])), RR, T);
                    c[dec] := Y[1];
                    S:=Y[2]: 
                    null := Y[3]; 
                    notnull := Y[4]; 
                    if c[dec] = {} then 
                        if S <> 0 then 
                             v:=simplify(op(nrm([S])));
                             nr:=nops(R):
                             R:=[op(R),[u,nr+1,LM(v),1]]:
                             Polys:=[op(Polys),v]:
                             JP:=[op({op(JP),seq(JPair2(R[-1],R[j],LMH),j=1..nr)})]:
                             JP:=sort(JP,tidy);
                         end if: 
                    elif S <> 0 then 
                        RETURN(false,Polys,null,notnull,[u,S],JP,R);
                    end if: 
            fi:
        fi:
od:
RETURN(true,Polys,null,notnull,0,{},{});
end:
#######################################
#######################################
canspec := proc (LL, M, R) 
local N, W, WW, WWW, NN, NNN, test, i,h, t, facN, facW, x, NNNN, L, w, flag, ww, www, MM,NnN,n; 
#option trace;
N := Basis(LL,R); 
W := M; 
WW := seq(NormalForm(W[i], N, R), i = 1 .. nops(W));
WWW := seq(factor([WW][i]), i = 1 .. nops([WW])); 
test := true; 
t:= true; 
NN := N; 
h := product(W[i], i = 1 .. nops(W)); 
if RadicalMembership(h, `<,>`(NN)) = true then
    test := false;
    NN := {1};
    RETURN([test, NN,[WW]]):
end if; 
while t do
    t := false;
    facN := [seq([op(i)], `in`(i, [seq(fac([NN[i]]), i = 1 .. nops(NN))]))];
    facW := [op(fac([WWW], R))];
    NNN := [];
    L := NULL;
    for i from 1 to nops(facN) do
        for x in facN[i] do
           flag := true;
           for w in facW while flag do
             if divide(w, x) then
                flag := false;
                L := L, x:
             end if:
           end do:
           
        end do;
        NnN[i]:= [op(`minus`({op(facN[i])}, {L}))];
        n[i]:=simplify(expand(product(NnN[i][j], j = 1 .. nops(NnN[i]))));
     end do;   
     NNNN :=[seq(n[i],i=1..nops(facN))];
    if {op(NNNN)} <> {op(NN)} then 
        t := true;
        NN := Basis(NNNN, R);
        WW := seq(NormalForm([WWW][i], NN, R), i = 1 .. nops([WWW]));
        WWW := seq(factor([WW][i]), i = 1 .. nops([WW]))
    end if:
end do; 
ww := NULL; 
for i to nops([WWW]) do
    if type(WWW[i], 'constant') = true then
    else
        ww := ww, [WWW][i]:
    end if
end do; 
WWW := ww; 
MM := [op(fac([WWW], R))]; 
WWW :=[op({seq(i/LeadingTerm(i, R)[1], `in`(i, MM))})]; 
www :=[seq(nrm([ee]), `in`(ee, WWW))]; 
if www <> [] then
    www := [op(www)];
    www := [seq(op(jj), `in`(jj, www))]:
end if; 
RETURN([test, NN, www]): 
end proc:

##########################################
branch := proc (v, ffff,BBB, NNN, WWW, R, T,CPP,RR) 
local cd,l,ll,i,ff,Y,X,BB,pivot,test,TTT,Bb,g,TT,vv,LL,gg,B,CN,N,W,NN,WW,f,CP,RRRR;
global LIST, Grob,termord1, termord2, termord3, iii, DDD,ind,rdz; 
#option trace:
f:=ffff[2];
CN := canspec(NNN, WWW, R); 
B:=BBB:
N:=CN[2];
W:=CN[3];
cd := []; 
Bb:=B: 
Y := new(N, W, f, R, T);
cd := [op(Y[1])];
ff := simplify(expand(Y[2])); 
NN := Y[3]; 
WW := Y[4];  
if v = [] then 
    vv := NULL: 
else 
    vv := op(v): 
end if;  
TT[vv] := [v, Bb, NN, WW]; 
if cd = [] then 
    X:= gvw([ffff[1],ff],Bb, NN, WW, R, T,CPP,RR); 
    test := X[1]; 
    BB := X[2]; 
    NN := X[3]; 
    WW := X[4]; 
    ff:=X[5];
    CP:=X[6];
    RRRR:=X[7]:
    if test  then 
        #TT[vv] := [[vv], nrm(InterReduce(BB,T)), NN, WW]; 
        TT[vv] := [[vv], nrm(BB), NN, WW]; 
    if CP={} then
        LIST := LIST, TT[vv]:
        
    fi;
    else 
        branch([vv],ff, BB, NN, WW, R, T,CP,RRRR); 
    end if: 
else 
    branch([0,vv], [ffff[1],ff],Bb, [op(N),op(cd)], W, R, T,CPP,RR); 
    branch([1,vv], [ffff[1],ff],Bb, N, [op(W),op(cd)], R, T,CPP,RR);
end if; 
end proc:
################## Incremental Montes Algorithm (IncDisPGB algorithm)##############
IncMontes:=proc(f,Grobsys,R,T)
local L,i,G,N,W,ff,CPP;
global LIST, DDD,rdz,Grob;
#option trace;
L:=NULL:
for i from 1 to nops(Grobsys) do 
    CPP:={};
    LIST:=NULL:
    G:=Grobsys[i];
    Grob:=G[2]:
    N:=G[3]:
    W:=G[4]:
    ff:=NormalForm(f,N,R);
    if ff=0 then
        LIST:=LIST,[G[1],Grob,G[3],G[4]]: 
    else 
     branch([], [1,ff],Grob, N, W, R, T, {},[seq([0,i,LM(Grob[i]),1],i=1..nops(Grob))]);
    fi:
    L:=L,LIST;
od:
LIST:=L;
RETURN(LIST);
end:

#######################
selpoly1 := proc (f, g) 
local ZPf, ZPg;
global termord1, termord2; 
ZPf := LeadingCoefficient(f, termord2); 
ZPg := LeadingCoefficient(g, termord2); 
RETURN(TestOrder(ZPg, ZPf, termord1)): 
end proc:
################################### GVWDisPGB Algorithm ##############
GVWDisPGB:=proc(FF,R,T)
local f,t1,t2,b1,b2,F;
global LIST,DDD,ind,rdz,Grob,termord1, termord2;
local IndPolys,POLYS;
#option trace;
t1,b1:=kernelopts(cputime,bytesused);
termord1, termord2:=R,T:
IndPolys:=[seq(i,i=1..nops(FF))];
F:=InterReduce(FF,prod(T,R));
POLYS:=FF;
rdz:=0;
DDD:=0;
LIST:=NULL;
f:=F[1];
Grob:=[];
ind:=nops(F)-1;
branch([], [1,f],[], [], [], R, T,{},[]);
for f in F[2..-1] do
    ind:=ind-1;print(ind);
    IncMontes(f,[LIST],R,T); 
od;
t2,b2:=kernelopts(cputime,bytesused);
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        : %3g\n",Num,of,Rds,rdz):
 printf("%-1s %1s %1s        : %3a\n",Num,of,sys,nops([LIST])):
 printf("%-1s %1s %1s %1s    : %3a\n",The, max, deg, is,DDD):
RETURN(LIST);
end:
                                          




with(Groebner):with(PolynomialTools):with(PolynomialIdeals):
############################################################################################################################
                                          ##########################
                                          #          GVW           #
                                          ##         |  |         ##
                                          ###        |  |        ###
                                          ####      _|  |_      ####
                                          #####     \    /     #####
                                          ######     \  /     ######
                                          #######     \/     #######
                                          ########          ########
                                          ##########################
########################################## AutoReduced Algorithm #################
AutoReduction2:=proc(f,F,TermOrder)
#option trace;
local LMF,Vars,p,r,flag,g,u,Mult,InvDiv,i,FF;
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:

##########################
##########################
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,FF,i,j;
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:
##############################
Invgvw:=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,Q,Q2,q,g,p,h,lmp,flag,G:
#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();
RETURN(G);
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:

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]:
Invgvw(F,tdeg(x,y,z),J):

############################################################################################################################
############################################################################################################################
                                          ##########################
                                          #         PGBMain        #
                                          #        algorithm       #
                                          ##         |  |         ##
                                          ###        |  |        ###
                                          ####      _|  |_      ####
                                          #####     \    /     #####
                                          ######     \  /     ######
                                          #######     \/     #######
                                          ########          ########
                                          ##########################
#################################################### Notice #################################################################################
#### The output Grobner system may contain several branches so that the corresponding Grobner basis is {1}. By combining all these branches##
#### in only one branch one can reduce the consistency check and this can improve significantly the performance of the algorithm. So, we ####
#### removed the fifth line of the original Kapur algorithm (where the branches with the Grobner basis {1} are detected), and instead, ######
#### we can add the branch "(Otherwise, {1})" into the output.#################################################################################
#############################################################################################################################################
with(LinearAlgebra):with(CodeTools):with(Groebner): with(Algebraic): with(PolynomialIdeals):
#####################PRIMEPART Algorithm###################
ppart := proc (K, T) 
local h, i, ti, o, oo, ooo;
if {op(K)}<>{0}and{op(K)} <> {} then
    h := LeadingTerm(K[1], T)[1];
    for i from 2to nops(K) do
        ti := LeadingTerm(K[i], T)[1];
        h := gcd(h, ti)
    end do;
    RETURN(simplify(expand(K/h)))
else
    RETURN(K)
end if end proc:

 #########################LT- Algorithm######################
lt := proc (f, T)
    RETURN(LeadingTerm(f, T)[1]*LeadingTerm(f,T)[2]):
end proc:
#########################Normalize Algorithm################
nrm := proc (F) 
#option trace;
local KK, h, i, NM:
if F = [] then
    RETURN(F)
end if;
KK := F;
h := denom(KK[1]);
for i from 2 to nops(KK) do
    h := simplify(lcm(denom(KK[i]), h));
end do;
NM := simplify(expand(h*KK));
RETURN(NM):
end proc:
#####################################convert a polynomial into the list of its terms############
FL := proc (f, T) 
local L, p; 
L := []; 
p := f; 
while p <> 0 do
    L := [op(L), lt(p, T)];
    p := simplify(p-lt(p, T))
end do; 
RETURN(L) 
end proc:
###################################convert a list of polynomial to the set "FACVAR"##########
fac := proc (L, T)
local A, AA, AAA, N, P, p, B, C, i;
A := L;

AA:= [seq(factor(i), `in`(i, A))];
AAA := NULL;
B := [];

for i to nops(AA) do
    p := AA[i];
    if irreduc(p) = true then
        AAA := AAA, p:
    else
        AAA := AAA, op(convert(p, list)):
    end if
end do;
P := NULL;
for i to nops([AAA]) do
    if `not`(type(AAA[i], 'constant')) then
        P := P, [AAA][i]:
    end if
end do;
N := NULL;
for i to nops([P]) do
    if irreduc([P][i]) = false and type([P][i], 'constant') <> true then
        N := N, op([P][i])[1]:
    else
        N := N, [P][i]:
    end if
end do;
RETURN({N}):
end proc:
##############################################CANSPEC ALGORITM############################
canspec := proc (LL, M, R) 
local N, W, WW, WWW, NN, NNN, test, i,h, t, facN, facW, x, NNNN, L, w, flag, ww, www, MM,n,NnN; 
#option trace;
N := Basis(LL,R); 
W := M; 
WW := seq(Groebner:-NormalForm(W[i], N, R), i = 1 .. nops(W));
WWW := seq(factor([WW][i]), i = 1 .. nops([WW])); 
test := true; 
t:= true; 
NN := N; 
h := product(W[i], i = 1 .. nops(W)); 
if RadicalMembership(h, `<,>`(NN)) = true then
    test := false;
    NN := {1};
    RETURN(test, NN):
end if; 
while t do
    t := false;    
    facN := [seq([op(i)], `in`(i, [seq(fac([NN[i]]), i = 1 .. nops(NN))]))];
   ### facN := [op(fac(NN, R))];
    facW := [op(fac([WWW], R))];
    NNN := [];
    L := NULL;
    for i from 1 to nops(facN) do
        for x in facN[i] do
           flag := true;
           for w in facW while flag do
             if divide(w, x) then
                flag := false;
                L := L, x:
             end if:
           end do:
           
        end do;
        NnN[i]:= [op(`minus`({op(facN[i])}, {L}))];
        n[i]:=simplify(expand(product(NnN[i][j], j = 1 .. nops(NnN[i]))));
     end do;   
     NNNN :=[seq(n[i],i=1..nops(facN))];  
     if {op(NNNN)} <> {op(NN)} then 
        t := true;
        NN := Basis(NNNN, R);
        WW := seq(NormalForm([WWW][i], NN, R), i = 1 .. nops([WWW]));
        WWW := seq(factor([WW][i]), i = 1 .. nops([WW]))
     end if:
end do; 
ww := NULL; 
for i to nops([WWW]) do
    if type(WWW[i], 'constant') = true then
    else
        ww := ww, [WWW][i]:
    end if
end do; 
WWW := ww; 
MM := [op(fac([WWW], R))]; 
WWW :=[op({seq(i/LeadingTerm(i, R)[1], `in`(i, MM))})]; 
www :=[seq(nrm([ee]), `in`(ee, WWW))]; 
if www <> [] then
    www := [op(www)];
    www := [seq(op(jj), `in`(jj, www))]:
end if; 
if www=[1] then
www:=[];
fi;
RETURN([test, NN, www]): 
end proc:
#########################Consistent ALGORITM############################
consist := proc (LL, M, R) 
local N, W, WW, WWW, NN, NNN, test, i,h,t, facN, facW, x, NNNN, L, w, flag, ww, www, MM;
#option trace;
h := product(M[i], i = 1 .. nops(M));
if RadicalMembership(h, `<,>`(LL)) = true then
    #NN := {1};
    RETURN(false, LL):
end if;
RETURN([true, LL, M]):
end proc:  
#################################################N*M######compute the binary product of two list L and M (new version) that is used#######################          
zarb:=proc(A,B)
RETURN(simplify(expand([seq(seq(u*v, u = A), v = B)])));
end:
########################################################## Minimal Dickson Basis ######################################B
MDBasis:=proc(F,T)
local g,N,L,i,flag;
#option trace;
N:=NULL;
for i from 1 to nops(F) do
    flag:=false;
    for g in [N,op(F[i+1..-1])] while not flag do 
        if divide(LeadingMonomial(F[i],T),LeadingMonomial(g,T)) then 
            flag:=true;
        fi;
    od;
    if not flag then
        N:=N,F[i];
    fi;
od;
RETURN([N]);
#RETURN(Basis([N],T));
end:
############################################### Zero Dim-Check ################################################OK
Zdim := proc (E, f, R) 
local ns, rv, poly, d, M, EE, Inf, Rr, EG,ms,sms; 
#option trace;
EE := `<,>`(op(E)); 
Inf := IdealInfo[Variables](EE); 
ms:=[op({op(R)} minus Inf)];
sms:=subs(seq(ms[i]=NULL,i=1..nops(ms)),R);
Rr:=sms; 
#EG := Groebner:-Basis(E, R); 
EG := E;
ns, rv := NormalSet(EG, Rr); 
M := MultiplicationMatrix(f, ns, rv, EG, Rr); 
poly := CharacteristicPolynomial(M, x); 
d := degree(poly, x); 
if poly <> x^d then 
   RETURN(true): 
else 
   RETURN(false): 
end if: 
end proc:
###################################################################### CCheck sub-algorithm ###############################
Ccheck := proc (EE, f, R) 
local LM, V, d, abar, EV, spE, E,Rr,LE,ff; 
#option trace;
E := [op(EE)]; 
#Rr:=tdeg(op(R));
LE:=LeadingMonomial(E,R);
V:=Groebner:-MaximalIndependentSet(LE,{op(R)}); 
d := nops(V); 
abar := [seq(i, i = 1 .. d)]; 
ff:=subs(seq(V[i] = abar[i], i = 1 .. nops(V)), f);
EV:=subs(seq(V[i] = abar[i], i = 1 .. nops(V)), E);
spE := Groebner:-Basis([EV][nops([EV])], R); 

   if IsZeroDimensional(`<,>`(op(spE))) then 
      if Zdim(spE, [ff][-1], R) then
         RETURN(true); 
      end if; 
   end if; 
 
RETURN(false); 
end proc:
########################################################## ICheck sub-algorithm ######################################################
Icheck := proc (E, f, R) 
local loop, p, i, s, M, m,Rr; 
#option trace;
Rr:=tdeg(op(R));
loop := 1; 
p := f; 
for i to loop do 
   M := FL(p, R); 
   M:=LeadingMonomial(M,R);
   s := 0; 
   for m in M do 
      s := simplify(expand(s+NormalForm(p*m, E, Rr))); 
   end do; 
   if s = 0 then 
      RETURN(true); 
   end if; 
   p := s; 
end do; 
RETURN(false); 
end proc:
################################################## Consistent  Kapur tests(Zdim, Ccheck, Icheck) for a polynomial################
consist2 := proc (E, f, R) 
local EE, Inf, Rr, flag1, flag2; 
#option trace;
if nops(E) = 0 then 
   RETURN(true) 
end if; 
if IdealMembership(f, `<,>`(op(E))) then 
   RETURN(false) 
end if; 
EE := `<,>`(op(E)); 
Inf := IdealInfo[Variables](EE); 
Rr:=op(0,R)(op(Inf));

if IsZeroDimensional(EE) then 
   RETURN(Zdim(E, f, R)); 
end if; 
flag1 := Ccheck(E, f, R); 
if flag1 then  
   RETURN(true); 
end if; 
#flag2 := Icheck(E, f, R); 
#if flag2 then 
#   RETURN(false); 
#end if; 
RETURN(consist(E, [f], R)[1]); 
end proc:
########################################consist1 new(check consistency of E and product of element of N)############################
consist1 := proc (EE, N, R) 
local f;
#option trace;
f:=mul(i,i in N);
RETURN(consist2(EE,f,R)); 
end proc:
#########################################################Combination of Canspec and Consist checks#################################
ccc := proc (LL, M, R) 
local NnN, n, t, facN, facW, x, NNNN, L, w, ww, www, MM, N, W, WW, WWW, NN, NNN, test, flag, i; 
#option trace;
if nops(LL) = 0 then 
   RETURN([true, LL, M]); 
end if; 
N:=LL;
if M = [] or M={} or M=[1] then 
   if Groebner:-Basis(LL,R)<>[1] then
      RETURN([true, LL, M]); 
   fi;
end if;
 
W := M; 
WW := seq(NormalForm(W[i], N, R), i = 1 .. nops(W)); 
WWW := seq(factor([WW][i]), i = 1 .. nops([WW])); 
test := true; 
t := true; 
NN := N; 
i := 1; 
flag := consist1(N, M, R);  
if flag = false then 
   RETURN([flag, N, [WWW]]) 
else 
   while t do 
      t := false; 
      facN := [seq([op(i)], `in`(i, [seq(fac([NN[i]]), i = 1 .. nops(NN))]))]; 
      facW := [op(fac([WWW], R))]; 
      NNN := []; 
      L := NULL; 
      for i to nops(facN) do 
         for x in facN[i] do 
            flag := true; 
            for w in facW while flag do 
               if divide(w, x) then 
                  flag := false; 
                  L := L, x; 
               end if; 
            end do; 
         end do; 
         NnN[i] := [op(`minus`({op(facN[i])}, {L}))]; 
         n[i] := simplify(expand(product(NnN[i][j], j = 1 .. nops(NnN[i])))); 
      end do; 
      NNNN := [seq(n[i], i = 1 .. nops(facN))]; 
      #if {op(NNNN)} <> {op(NN)} then 
      if flag=false then
         t := true; 
         NN := Groebner:-Basis(NNNN, R); 
         WW := seq(NormalForm([WWW][i], NN, R), i = 1 .. nops([WWW])); 
         WWW := seq(factor([WW][i]), i = 1 .. nops([WW])); 
      end if; 
   end do; 
end if; 
ww := NULL; 
for i to nops([WWW]) do 
   if type(WWW[i], 'constant') = true then  
   else 
      ww := ww, [WWW][i]; 
   end if; 
end do; 
WWW := ww; 
MM := [op(fac([WWW], R))]; 
WWW := [op({seq(i/LeadingTerm(i, R)[1], `in`(i, MM))})]; 
www := [seq(nrm([ee]), `in`(ee, WWW))]; 
if www <> [] then 
   www := [op(www)]; 
   www := [seq(op(jj), `in`(jj, www))]; 
end if; 
RETURN([flag, NN, www]); 
end proc:
#####################################################input: parametric polynomial f, output: variables appeared in f###############################
Variables := proc (e, p::(({list, rtable, set})(name))) 
`minus`(indets(e, And(name, Not(constant))), convert(p, set)):
end proc:
#####################################################input: parametric poly ideal F, output: parameters appeared in F ###############################
pars := proc (F, T) 
local p, Ps; 
#option trace;
p := [op(T)]; 
Ps := `union`(seq(Variables(F[i], p), i = 1 .. nops(F))); 
RETURN(Ps); 
end proc:
#####################################################input: parametric poly ideal F, output: variables appeared in F###############################
vars := proc (F, R) 
local p, Vs; 
p := [op(R)]; 
Vs := `union`(seq(Variables(F[i], p), i = 1 .. nops(F))); 
RETURN(Vs); 
end proc:
##################################################################Return true if a is constant#########################################
isconstant := proc (a) 
RETURN(type(a, constant)); 
end proc:
##############################################Return true if there is some variable in parametric polynomial f and vice versa##############################
issubsetv:=proc(f)  
local var;
global ordp,ordv;  
var:={op(ordv)};  
if indets(f) subset var then      
   RETURN(true);  
else       
   RETURN(false);  
fi;  
end:
####################################################Return true if there is any variable in parametric polynomial f and vice versa##########################
issubsetp:=proc(f)  
local par;
global ordp,ordv;  
#option trace;
par:={op(ordp)};  
if indets(f) subset par then      
   RETURN(true);  
else       
   RETURN(false);  
fi;  
end:
###################################################use in the end of PGBMain algo. for product of h_i's###################################################
cunion := proc (H, E, N) 
local EE, NN, k, List, i, h; 
#option trace;
k := nops(H); 
List := NULL; 
for i to k do 
   EE := E; 
   NN := N; 
   if i = 1 then 
      h := 1; 
   else 
      h := product(H[j], j = 1 .. i-1); 
   end if; 
   EE := [op(EE), H[i]]; 
   NN := zarb(NN, [h]); 
   List := List, [EE, NN]; 
end do; 
RETURN(List); 
end proc:


#############################################################PGBmain algorithm(without branch [1])######################KAPUR
pgbmainc := proc (EE, N, F, R, T) 
local BGB,G, Gg, i, j, Nj, Hj, hj, GG, Gm, Gr, Z, Zh, H, GB, h, Grr, Grrr, C, CC, CCC, CCCC,g,GgG, U,Zj,E,k,CU; 
global ordp,ordv,PGB,LIST;
#option trace;
ordp:=R;
ordv:=T;
E:=Groebner:-Basis(EE,R);
CC := ccc(E, N, R); 
if CC[1] = false then 
   RETURN(); 
else 
#G:=Groebner:-Basis([op(F),op(CC[2])],prod(T,R)); 
BGB:=Invgvw([op(F),op(CC[2])],prod(T,R));
G:=InterReduce(BGB,prod(T,R));
end if;  
if pars(F,T)={} then 
   RETURN([E,N,F]);
fi;
Gr:=select(issubsetp,G);
 
GG:=remove(issubsetp,G);
CCC := ccc(Gr, [op(N)], R); 
if CCC[1] = false then 
   RETURN(PGB); 
else 
   Gm := MDBasis(GG, T);  
   H := [op({seq(LeadingCoefficient(Gm[i], T), i = 1 .. nops(Gm))})];
   H := remove(isconstant, H); 
   h:=lcm(op(H));
   Zh := zarb(N, [h]); 
   CCCC := ccc(Gr, [op(Zh)], R); 
   if CCCC[1] then 
      PGB := PGB, [CCCC[2], factor(CCCC[3]), Gm];
      CU:=[cunion(H,CCCC[2],N)];
   else
      CU:=[cunion(H,CCC[2],CCC[3])];
   end if; 
   LIST:=[seq(pgbmainc(CU[i][1],CU[i][2],GG,R,T),i=1..nops(CU))];
end if;
RETURN(PGB); 
end proc:
####################################PGBMAIN based on above pgbmainc algorithm(we call this for computation of Grobner system)#########################
PGBMAIN := proc(E,N,F,R,T)
local AA;
global PGB,LIST,numzero;
numzero:=0;
PGB:=NULL;
AA:=pgbmainc(E,N,F,R,T);
#RETURN(LIST);
end:
##############################################################Example#########################
#PGBMAIN([], [1], [a*x^2*y+b*x+y^3, a*x^2*y+b*x*y, y^2+b*x^2*y+x*y], plex(a, b, c, d), plex(x, y, z)):
############################################################################################################################