##################################################################################################################FINAL#####
############################################################################################################################
############################################################################################################################
                                          ##########################
                                          #          DisPGB        #
                                          #        algorithm       #
                                          ##         |  |         ##
                                          ###        |  |        ###
                                          ####      _|  |_      ####
                                          #####     \    /     #####
                                          ######     \  /     ######
                                          #######     \/     #######
                                          ########          ########
                                          ##########################
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) 
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 to the list of it 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:
                      ################PRIMITIVATION one list 'G' of polynomials################
plp := proc (G, T) 
local FFFF, FFF, FF, FFi, Gi, i, j, gj, k; 
FFFF:= []; 
FFF := []; 
for i to nops(G) do
    Gi := FL(G[i], T);
    FFi := ppart(Gi, T);
    FFF := [op(FFF), FFi]
end do; 
for j to nops(FFF) do
    gj := 0;
    for k to nops(FFF[j]) do
        gj := gj+FFF[j][k]
    end do;
    FFFF := [op(FFFF), gj]
end do; 
RETURN(FFFF) 
end proc:
                      ##############################GGE ALGORITM###############################
#gge := proc (F, T) 
#RETURN(InterReduce(F, T)): 
#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,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):
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; 
if www=[1] then
www:=[];
fi;
RETURN([test, NN, www]): 
end proc:
                      #######################NEW COND ALGORITM##########################
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)})]; 
       # NN := Basis(NN, R);###
        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 := sum(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:
#RETURN([{op(nrm([op(cd)]))}, ff+k, NN, www]); 
end proc:
                      ############################Select Algoritm###################### 
sff := proc (x) 
evalb(type(x, 'constant')): 
end proc:
                      ############################Select Algoritm######################
sf := proc (x)
evalb(x <> 0):
end proc:
                      ###########################Normal Sterategy######################
f2 := proc (p, q)
local t;
global setB, Tord;
t := Tord;
RETURN(TestOrder(lcm(lt(setB[p[1]], t), lt(setB[p[2]], t)),lcm(lt(setB[q[1]], t), lt(setB[q[2]], t)), t))
end proc:
                      ###########################normal strategy########################
f2 := proc (p, q)
local t;
global setB, Tord;
t := Tord;
RETURN(TestOrder(lcm(lt(setB[p[1]], t), lt(setB[p[2]], t)), lcm(lt(setB[q[1]], t), lt(setB[q[2]], t)), t))
end proc:

#################################################################  mini algorithm ##################

mini := proc (F, T)
local H,i,f,HH,test,j;
  H := F;
  for i to nops(H) do 
    f := H[1]; 
    #H := subs(H[1] = NULL, H); 
     H:=H[2..-1];
    HH := [seq(LeadingMonomial(H[i], T), i = 1 .. nops(H))]; 
    test := false; 
        for j to nops(H) do 
            if divide(LeadingMonomial(f, T), HH[j]) = true then 
                test := true 
            end if 
        end do; 
        if test = false then 
            H := [op(H), f] 
        end if 
  end do;
  RETURN(H) 
  end proc:

                      ###########################CONDPGB(with normal strategy)########################
condpgb := proc (B, N, W, R, T, JJ)
local test, NN, WW, BB, i, j, t, s, J, x, c, Y, S, ii,SOP,felag;
global setB, Tord, DDD, iii;
#option trace;

SOP:=JJ;

felag:=true;
test := true;
s := nops(B);
J := [];
setB := B;
Tord := T;
if SOP=[] then
for i to s do
    for j from i+1 to s do
         SOP:=[op(SOP), {i, j}];
    end do
end do;
fi;
BB := B;
NN := N;
WW := W;
t := s;
#J := sort(J, f2);
SOP:= sort(SOP, f2);
#while J <> [] and test do
while SOP <> [] and test do
    #x := J[1];
     x := SOP[1];
    #J := J[2 .. -1];
     SOP := SOP[2 .. -1];
    i := x[1];
    j := x[2];
    S := expand(simplify(NormalForm(op(nrm([SPolynomial(BB[i], BB[j], T)])), select(sf, BB), T)));
    S := simplify(NormalForm(op(nrm([S])), NN, prod(T,R)));
    if S<>0 then 
       felag:=false;
    fi;
    DDD := max(DDD, degree(lcm(BB[i], BB[j]), {op(T)}));
    if S <> 0 then 
        Y := new(NN, WW, S, R, T);
        c[dec] := Y[1];
        S := Y[2];
        NN := Y[3];
        WW := Y[4];
        if c[dec] = {} then
            if S <> 0 then
                t := t+1;
                BB := [op(BB), S];
                setB := BB;
                
                SOP := [op(SOP), seq({t, ii}, ii = 1 .. t-1)];
                SOP := sort(SOP, f2):
            end if
        #elif S <> 0 then
        else
            test := false;
            BB := [op(BB), S]:
            #SOP:=[op(SOP),S]
        end if:
    else
        iii := iii+1:
    end if:
end do;
if test then 
    BB := InterReduce(mini(BB,T), prod(T, R)): 
end if;
RETURN([test, BB, NN, WW, SOP]):
end proc:

                      #######################Branch Algoritm#######################
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:
selpoly2 := proc (f, g)
local ZPf, ZPg;
global termord1, termord2;
ZPf := LeadingCoefficient(f, termord1);
ZPg := LeadingCoefficient(g, termord1);
RETURN(TestOrder(ZPg, ZPf, termord2)):
end proc:
selpoly3 := proc (f, g)
local ZPf, ZPg;
global termord1, termord2;
ZPf := LeadingCoefficient(f, termord2);
ZPg := LeadingCoefficient(g, termord2);
RETURN(evalb(degree(ZPg, termord1) < degree(ZPf, termord1))):
end proc:
selpoly4 := proc (f, g)
local ZPf, ZPg;
global termord1, termord2;
ZPf := LeadingCoefficient(f, termord2);
ZPg := LeadingCoefficient(g, termord2);
RETURN(evalb(degree(ZPf, termord1) < degree(ZPg, termord1))):
end proc:
selpoly5 := proc (f, g)
local ZPf, ZPg;
global termord1, termord2;
ZPf := LeadingCoefficient(f, termord1);
ZPg := LeadingCoefficient(g, termord1);
RETURN(evalb(degree(ZPg, termord2) < degree(ZPf, termord2))):
end proc:
selpoly6 := proc (f, g)
local ZPf, ZPg;
global termord1, termord2;
ZPf := LeadingCoefficient(f, termord1);
ZPg := LeadingCoefficient(g, termord1);
RETURN(evalb(degree(ZPf, termord2) < degree(ZPg, termord2))):
end proc:
selpoly7 := proc (f, g)
local ZPf, ZPg;
global termord1, termord2, termorder3;
ZPg := lt(g, termord3);
RETURN(evalb(degree(ZPg, termord3) < degree(ZPf, termord3))):
end proc:
selpoly8 := proc (f, g)
local ZPf, ZPg;
global termord1, termord2, termorder3;
ZPf := lt(f, termord3);
ZPg := lt(g, termord3);
RETURN(evalb(degree(ZPf, termord3) < degree(ZPg, termord3))):
end proc:
selpoly9 := proc (f, g)
local ZPf, ZPg;
global termord1, termord2;
ZPf := LeadingCoefficient(f, termord1);
ZPg := LeadingCoefficient(g, termord1);
RETURN(divide(ZPg, ZPf)):
end proc:
selpoly10 := proc (f, g)
local ZPf, ZPg;
global termord1, termord2;
ZPf := LeadingCoefficient(f, termord1);
ZPg := LeadingCoefficient(g, termord1);
RETURN(divide(ZPf, ZPg)):
end proc:
selpoly11 := proc (f, g)
local ZPf, ZPg;
global termord1, termord2;
ZPf := lt(f, termord1);
ZPg := lt(g, termord1);
RETURN(divide(ZPg, ZPf)):
end proc:
selpoly12 := proc (f, g)
local ZPf, ZPg;
global termord1, termord2;
ZPf := lt(f, termord1);
ZPg := lt(g, termord1);
RETURN(divide(ZPf, ZPg)):
end proc:
##############################                   
branch := proc (v, BBB, N, W, R, T,JJ) 
local SOP,cd, l, ll, i, f, ff, Y, X, BB, NN, WW, pivot, test, TTT, Bb, g, TT, vv, LL,gg,B; 
global LIST, termord1, termord2, termord3, iii, DDD; 
#option trace:
cd := []; 
termord1 := R; 
termord2 := T; 
termord3 := prod(termord2, termord1); 
SOP:=JJ;

B:=BBB;
B := sort(BBB, selpoly1);   
for i from 1 to nops(B) while cd = [] do
Bb := B;

    f := B[i]; 
    Y := new(N, W, f, R, T); 
    cd := [op(Y[1])]; 
    ff := simplify(expand(Y[2])); 
    NN := Y[3]; 
    WW := Y[4]; 
    Bb[i] := Y[2]; 
    if cd <> [] then 
        pivot := i: 
    end if: 
end do; 
if nops(B)=0 then
    NN := N; 
    WW := W;
fi:
LL := NULL;

for i to nops(Bb) do 
    if Bb[i] <> 0 then 
        LL := LL, Bb[i]:
    end if: 
end do; 
Bb := [LL];
if v = [] then 
   vv := NULL: 
else 
    vv := op(v): 
end if; 
gg := B[pivot]; 
if B=[] then
Bb:=[];
fi;
TT[vv] := [v, Bb, NN, WW]; 
if cd = [] then 
    
    X := condpgb(Bb, NN, WW, R, T,SOP); 
    test := X[1]; 
    BB := X[2]; 
    NN := X[3]; 
    WW := X[4]; 
    if test then 
        TT[vv] := [[vv], BB, NN, WW]; 
        LIST := LIST, TT[vv]: 
    else 
        branch([vv], BB, NN, WW, R, T, SOP) 
    end if: 
else 
    newvertex(1, [vv], cd, Bb, NN, WW, pivot, R, T, gg,SOP); 
    newvertex(0, [vv], cd, Bb, NN, WW, pivot, R, T, gg,SOP):
end if; 
RETURN([LIST], iii, DDD): 
end proc:
                     ##########################New vertex Algoritm#########################
newvertex := proc (n, u, cd, B, N, W, pivot, R, T, gg, JJ) 
local SOP,i, r,v, WW, NN, BB, cond, CN, TT, vv; 
#option trace;
SOP:=JJ;
if u = 0 then
    vv := [NULL]
else
    vv := op(u), n
end if; 
BB := simplify(B); 
r := product([op(cd)][i], i = 1 ..nops(cd)); 
if n = 0 then
    cond := r;
    WW := W;
    NN := Basis([op(cd), op(N)], R);
else
    cond := r;
    WW := {op(cd), op(W)};
    NN := N;
    BB := B; 
    if new(NN, WW, lt(gg, T), R, T)[1]={} then #####
    for i to nops(B) do
         if B[i] <> gg and divide(LeadingMonomial(B[i], T), LeadingMonomial(gg, T)) = true then
           BB := subs(B[i] = expand(simplify(SPolynomial(B[i], gg, T))), BB);
         
        end if
    end do
    end if;#####
end if; 
CN := canspec(NN, WW, R); 
if CN[1] then
   
    BB:=simplify(InterReduce([seq(NormalForm(simplify(BB[i]),CN[2],R),i=1..nops(BB))],prod(T,R))):
    TT[vv] := [cond, BB, NN, nrm(WW)];
    if vv = NULL then
        v := 0;
          branch(v, BB, CN[2], CN[3], R, T, SOP):#######
    else
         branch([vv], BB, CN[2], CN[3], R, T, SOP):##########
    end if:
else
    RETURN([vv, BB, {1}, WW])
end if: 
end proc:
                      ##############################DISPGB Algoritm########################                     
DISPG := proc (F, R, T) 
local B; 
global LIST, iii, DDD; 
LIST := NULL; 
iii := 0; 
DDD := 0; 
B := InterReduce(F, prod(T, R)):   
RETURN(branch([], B, [], [], R, T,[])): 
end proc: 
#################################################example###########################
#DISPG([a*x+y,b*y+c],plex(a,b,c),plex(x,y));