with(Groebner):
with(Algebraic):
with(PolynomialIdeals):
##########################################################
##################### Montes 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; 
#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; 
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)})]; 
        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:
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:
######Condpgb Algoritm########################
condpgb := proc (B, N, W, R, T) 
local test, NN, WW, BB, i, j, t, s, J, x, c, Y, S, ii; 
global setB, Tord, DDD, iii; 
#option trace;
test := true; 
s := nops(B); 
J := []; 
setB := B; 
Tord := T; 
for i to s do 
    for j from i+1 to s do 
        J := [op(J), {i, j}]: 
    end do: 
end do; 
BB := B; 
NN := N; 
WW := W; 
t := s; 
J := sort(J, f2); 
while J <> [] and test do 
    x := J[1]; 
    J := J[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, R)); 
    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; 
                J := [op(J), seq({t, ii}, ii = 1 .. t-1)];
            end if: 
        elif S <> 0 then 
            test := false; 
            BB := [op(BB), S]: 
        end if: 
    else 
        iii := iii+1: 
    end if: 
end do; 
if test then 
    BB:= InterReduce(BB, prod(T, R)): 
end if; 
RETURN([test, BB, NN, WW]) 
end proc: 

##########################################
branch := proc (v, BBB, N, W, R, T) 
local 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); 
B :=BBB;
Bb := B; 
for i from 1 to nops(B) while cd = [] do 
    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]; 
TT[vv] := [v, Bb, NN, WW]; 
if cd = [] then 
    X := condpgb(Bb, NN, WW, R, T); 
    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) 
    end if: 
else 
    newvertex(1, [vv], cd, Bb, NN, WW, pivot, R, T, gg); 
    newvertex(0, [vv], cd, Bb, NN, WW, pivot, R, T, gg): 
end if; 

RETURN([LIST], iii, DDD): 
end proc: 
####New vertex Algoritm#########################
newvertex := proc (n, u, cd, B, N, W, pivot, R, T, gg) 
local i, r,v, WW, NN, BB, cond, CN, TT, vv; 
#option trace;
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)
            
           # if BB[i] <> gg and divide(LeadingMonomial(BB[i], T), LeadingMonomial(gg, T)) = true then###############
           # BB := subs(BB[i] = expand(simplify(SPolynomial(BB[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):#######
    else
         branch([vv], BB, CN[2], CN[3], R, T):##########
    end if:
else
    RETURN([vv, BB, {1}, WW])
end if: 
end proc:
######DISPGB Algoritm########################
DISPG := proc (F, R, T) 
local B,t1,t2,b1,b2; 
global LIST, iii, DDD; 
t1,b1:=kernelopts(cputime,bytesused);
LIST := NULL; 
iii := 0; 
DDD := 0; 
B := InterReduce(F, prod(T, R)); 
A:=branch([], B, [], [], R, T);
t2,b2:=kernelopts(cputime,bytesused);
printf("%-1s %1s %1s             : %2g\n",The,CPU,time,t2-t1):
printf("%-1s %1s %1s          : %2g\n",The,used,memory,b2-b1):
printf("%-1s %1s %1s             : %2g\n",Num,of,Rds,iii):
 printf("%-1s %1s %1s         : %2a\n",The,max,degree,DDD):
RETURN(A): 
end proc:                                                                 
##########################################################
#############Improved Montes Algorithm ###################
##########################################################
###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; 
#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; 
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)})]; 
        ff := simplify(ff-simplify(expand(lt(ff, T)))): 
        #NN := [op({op(NN), LeadingCoefficient(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:=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:
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:
######Update Algoritm########################
update := proc (h, LL, BB) 
local L1, L, L2, L3, L4, L5, i, lth, flag, j, p, E,nn; 
global setB, Tord; 
#option trace;
lth := LeadingTerm(setB[h], Tord)[2]; 
L := LL; 
L1 := NULL; 
L2 := NULL; 
for i to nops(L) do 
    if gcd(LeadingTerm(setB[i], Tord)[2], lth) = 1 then 
        L1 := L1, i: 
    else 
        L2 := L2, i: 
    end if: 
end do; 
L1 := [L1]; 
nn:=nops(L1):
L2 := [L2]; 
for i in L2 do 
    flag := false;
    member(i,L2,'zz');
    L2:=subsop(zz=NULL,L2); 
    for j in L1 while not flag do 
        if divide(lcm(lth, LeadingTerm(setB[i], Tord)[2]), LeadingTerm(setB[j], Tord)[2]) then 
            flag := true: 
        end if: 
    end do; 
    for j in L2 while not flag do 
        if divide(lcm(lth, LeadingTerm(setB[i], Tord)[2]), LeadingTerm(setB[j], Tord)[2]) then 
            flag := true: 
        end if: 
    end do; 
    if not flag then 
        L1 := [op(L1), i]: 
    end if: 
end do; 
L5 := NULL; 
for p in BB do 
    if not divide(lcm(LeadingTerm(setB[p[1]], Tord)[2], LeadingTerm(setB[p[2]], Tord)[2]), lth) or lcm(LeadingTerm(setB[p[1]], Tord)[2], LeadingTerm(setB[p[2]], Tord)[2]) = lcm(LeadingTerm(setB[p[1]], Tord)[2], lth) or lcm(LeadingTerm(setB[p[1]], Tord)[2], LeadingTerm(setB[p[2]], Tord)[2]) = lcm(LeadingTerm(setB[p[2]], Tord)[2], lth) then 
        L5 := L5, p 
    end if: 
end do; 
if BB=[] and nops(LL)=1 then
    RETURN([{h,op(LL)}]);
fi;
E := [seq({i, h}, `in`(i, L1[nn+1..nops(L1)])), L5]; 
RETURN(sort(E, f2)): 
end proc:
#####CondPGB algorithm ####################
condpgb := proc (B, N, W, R, T,CPP)
local test, NN, WW, BB, i, j, t, s, J, x, c, Y, S, ii;
global setB, Tord, DDD, iii,inde;
#option trace;
felag:=true;
CP:=CPP:
inde:=true;
test := true;
s := nops(B);
if s=0 then RETURN([true, [], N, W,1,[]]);fi;
setB := B;
Tord := T;
BB := B;
NN := N;
WW := W;
t := s;
if CP={} then  
    for tt from 1 to nops(B) do
	CP:=update(tt, [seq(i, i = 1 .. tt-1)], CP):
   od:
else
    CP:=update(t, [seq(i, i = 1 .. t-1)], CP):
fi:    
J := sort(CP, f2);
while nops(J) <> 0 and test do
    x := J[1];
    J := J[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, 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;
                J := update(t, [seq(i, i = 1 .. t-1)], J): 
            end if
        elif S <> 0 then
            test := false;
            BB := [op(BB), S]:
        end if:
    else
        iii := iii+1:
    end if:
end do;
if test then 
    BB := InterReduce(BB, prod(T, R)):
end if;
if felag=true then 
   RETURN([test, B, NN, WW,S,J]):
else 
    RETURN([test, BB, NN, WW,S,J]):
fi;  
end proc:
#####Selection strategy for polynomials############
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:
######Branch algorithm ###############################
branch := proc (v, BBB, N, W, R, T,CPP) 
local 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,inde; 
#option trace:
cd := []; 
termord1 := R; 
termord2 := T; 
termord3 := prod(termord2, termord1); 
B:=BBB;
Bb := B; 
if inde then 
    Bb:=nrm([seq(NormalForm(h,N,R), h in Bb)]);
    CP:=CPP;
    f := B[-1];
    Y := new(N, W, f, R, T);
    if Y[2]=0 then 
	    fun:=proc (p, j) if j in p then RETURN(true) else RETURN(false) end if end proc;
	    CP:=remove(fun,CP,nops(B));
	    funi:=proc(p,j)  if p<j then RETURN(p); else RETURN(p-1); fi;end:  
	    CP:={seq({funi(p[1],nops(B)), funi(p[2],nops(B))}, p in CP)}:
    fi:
	cd := [op(Y[1])]; 
	ff := simplify(expand(Y[2])); 
	NN := Y[3]; 
	WW := Y[4]; 
	Bb[-1] := Y[2]; 
	if cd <> [] then 
	    pivot := nops(B): 
	end if: 
else
    CP:=CPP;
    for i from nops(B) to 1 by -1 while cd = [] do 
	f := B[i];
	Y := new(N, W, f, R, T);
	if Y[2]=0 then 
	    fun:=proc (p, j) if j in p then RETURN(true) else RETURN(false) end if end proc;
	    CP:=remove(fun,CP,i);
	    funi:=proc(p,j)  if p<j then RETURN(p); else RETURN(p-1); fi;end:  
	    CP:={seq({funi(p[1],i), funi(p[2],i)}, p in CP)}:
	fi:
	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;
fi: 
if nops(B)=0 then
    NN := N; 
    WW := W;
fi:
gg := Bb[pivot]; 
LL := NULL; 
for i to nops(Bb) do 
    if Bb[i] <> 0 then 
        LL := LL, Bb[i]: 
	fun:=proc (p, j) if j in p then RETURN(true) else RETURN(false) end if end proc;
	    CP:=remove(fun,CP,i);
	    funi:=proc(p,j)  if p<j then RETURN(p); else RETURN(p-1); fi;end:  
	    CP:={seq({funi(p[1],i), funi(p[2],i)}, p in CP)}:	
    end if: 
end do; 
Bb := [LL];
if v = [] then 
    vv := NULL: 
else 
    vv := op(v): 
end if; 
TT[vv] := [v, Bb, NN, WW]; 
if cd = [] then 
    X := condpgb(Bb, NN, WW, R, T,CP); 
    test := X[1]; 
    BB := X[2]; 
    NN := X[3]; 
    WW := X[4]; 
    CP:=X[6];
    if test then 
        TT[vv] := [[vv], BB, NN, WW]; 
        LIST := LIST, TT[vv]: 
    else 
        branch([vv], BB, NN, WW, R, T,CP) 
    end if: 
else 
    newvertex(1, [vv], cd, Bb, NN, WW, pivot, R, T, gg,CP); 
    newvertex(0, [vv], cd, Bb, NN, WW, pivot, R, T, gg,CP): 
end if; 
RETURN([LIST]): 
end proc: 
#######NewVertex algorithm ##################
newvertex := proc (n, u, cd, B, N, W, pivot, R, T, gg,CPP) 
local i, r,v, WW, NN, BB, cond, CN, TT, vv; 
#option trace;
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); 
CP:=CPP:
if CN[1] then
LLL:=NULL:
    TT[vv] := [cond, BB, NN, nrm(WW)]; 
    if vv = NULL then
        v := 0;
          branch(v, BB, CN[2], CN[3], R, T,CP):#######
    else
         branch([vv], BB, CN[2], CN[3], R, T,CP):##########
    end if:
else
    RETURN([vv, BB, {1}, WW])
end if: 
end proc:

######Impeoved DisPGB algorithm###########
DISPG := proc (F, R, T) 
local B,t,q; 
global LIST, iii, DDD,NR,NF,termord1,termord2,termord3,inde;
t1,b1:=kernelopts(cputime,bytesused);
NR:=0:
NF:=0: 
LIST := NULL; 
iii := 0; 
DDD := 0;
inde:=false;
termord1 := R; 
termord2 := T; 
termord3 := prod(termord2, termord1);  
B := sort(InterReduce(F, prod(T, R)), selpoly1);
J:={};
X:=branch([], B, [], [], R, T, J):
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 %1s    : %3a\n",The, max, deg, is,DDD):
 printf("%-1s %1s %1s %1s    : %3a\n",The, num, red, is,iii):
RETURN(X): 
end proc:
#####################################
#####################################
#######Big Reduction Number Algorithm
#####################################
#####################################
#restart:
with(PolynomialIdeals):
with(Groebner):
#####################################
Monom:=proc(Vars,d)
local i,j,L;
if d=0 then
    RETURN([1]);
else
    L:=Monom(Vars,d-1);
    RETURN([op({seq(seq(i*j,i=Vars),j=L)})]);
fi:
end:
#####################################
MulList:=proc(F,Vars,d)
#option trace;
local i,j,k,L,l,f;
l:=min(seq(degree(f),f=F));
if d<l then
   RETURN([]);
else
   L:=NULL:
   for f in F do
      if degree(f)<=d then
         L:=L,seq(expand(f*i),i=Monom(Vars,d-degree(f))); 
      fi:
   od:
   RETURN([L]);
fi:
end:
###############
COEFF:=proc(f,m)
#option trace;
local i,j,L;
global tord,Vars,v,dd,n;
if f=0 then
   RETURN(0);
else
   i:=LeadingTerm(f,tord):
   j:=i[1]*i[2]:
   if m=i[2] then
       RETURN(i[1]);
   else
       RETURN(COEFF(simplify(f-j),m));
   fi:
fi;
end:
################
MatrixForm:=proc(F)
#option trace;
local i,j,L,l,f,d,M,V;
global tord,Vars,v,dd,n;
d:=degree(F[1],Vars);
V:=Monom(Vars,d):
V:=sort(V, proc (a, b) options operator, arrow; TestOrder(b, a, tord) end proc);
#with(LinearAlgebra):
member(Vars[n-dd]^d,V,'k');;
M:=Matrix(nops(F),k,(i,j)->COEFF(F[i],V[j]));interface(rtablesize=infinity);
M:=GaussianElimination(M,'method'='FractionFree'):
RETURN(convert(M.convert(V[1..k],Vector),list));
end:
#################
LM:=proc(f)
global tord;
if f<>0 then
    RETURN(LeadingMonomial(f,tord));
fi:
RETURN(0);
end:
#################
Borel := proc (AA, Vars)
#option trace;
local i,a,e,u,Ind; 
    for a in AA do 
	Ind := indets(a); 
	i := 1; 
	for e in Ind do 
	    member(e, Vars, 'q'); 
	    i := max(i, q) 
	od; 
	for u in Vars[1 .. i] do 
	    if not IdealMembership(a*u/Vars[i], <op(AA)>) then 
		RETURN(false, a, u, Vars[i]) 
	    fi: 
	od: 
    od: 
RETURN(true) 
end:

Borel2:= proc (AA, Vars) 
local B,Var,i;
    B := AA; 
    Var := Vars; 
    for i from 0 to nops(Vars)-1 do 
	if Quotient(<op(B)>, <op(Var)>) <> Quotient(<op(B)>, <Vars[nops(Vars)-i]>) then 
	    RETURN(false, B, Vars[nops(Vars)-i]) 
	else 
	    B := subs(Vars[nops(Vars)-i] = 0, B); 
	    Var := Var[1 .. -2] 
	fi: 
    od: 
RETURN(true) 
end:


Nested:= proc (AA, Vars) 
local B,Var,i,Deg,polys,ee,dd,ListVar,mm,LL,bb,ii;
#option trace;
    B := Basis(AA,tdeg(op(Vars))); 
    Var:= Vars; 
    Deg:=max(seq(degree(polys), polys in B));
    ee:=0;
    for i from 1 to nops(Vars) do
	if <seq(Vars[j]^Deg,j=1..i)> subset <op(B)> then
	    ee:=i;
	fi;
    od:
    dd:=nops(Vars)-ee;
    for bb in B do
	ListVar:=indets(bb);
	 mm:=NULL:
	 for ii from 1 to nops(Vars) do
	    if Vars[ii] in ListVar then
		mm:=ii:
	    fi:
	od:
	for ii from 1 to mm-1 while mm>nops(Vars)-dd do
	    if not IdealMembership(bb*Vars[ii]^Deg/(Vars[mm]^degree(bb,Vars[mm])), <op(B)>) then
		RETURN(false,bb,bb*Vars[ii]^Deg/(Vars[mm]^degree(bb,Vars[mm])));
	    fi;
	od;
    od:
RETURN(true) 
end:

dstable2:= proc (AA, Vars) 
local B,i,Deg,polys,ee,dd,mm,LL,bb,ii;
#option trace;
    B := Basis(AA,tdeg(op(Vars))); 
    Deg:=max(seq(degree(polys), polys in B));
    ee:=0;
    for i from 1 to nops(Vars) do
	if <seq(Vars[j]^Deg,j=1..i)> subset <op(B)> then
	    ee:=i;
	fi;
    od:
    dd:=nops(Vars)-ee;
    for bb in B do
	 mm:=NULL:
	 for ii from 1 to nops(Vars) do
	    if Vars[ii] in indets(bb) then
		mm:=ii:
	    fi:
	od:
	for ii from 1 to mm-1 while mm>=nops(Vars)-dd do
	    if not IdealMembership(bb*Vars[ii]/Vars[mm], <op(B)>) then
		RETURN(false,bb,Vars[ii],Vars[mm]);
	    fi;
	od;
    od:
RETURN(true);
end:

Weakdstable:= proc (AA, Vars) 
local B,i,Deg,polys,ee,dd,mm,LL,bb,ii;
#option trace;
    B := Basis(AA,tdeg(op(Vars))); 
    Deg:=max(seq(degree(polys), polys in B));
    ee:=0;
    for i from 1 to nops(Vars) do
	if <seq(Vars[j]^Deg,j=1..i)> subset <op(B)> then
	    ee:=i;
	fi;
    od:
    dd:=nops(Vars)-ee;
    for bb in B do
	 mm:=NULL:
	 for ii from 1 to nops(Vars) do
	    if Vars[ii] in indets(bb) then
		mm:=ii:
	    fi:
	od:
	for ii from 1 to nops(Vars)-dd while mm>=nops(Vars)-dd do
	    if not IdealMembership(bb*Vars[ii]/Vars[mm], <op(B)>) then
		RETURN(false,bb,Vars[ii],Vars[mm]);
	    fi;
	od;
    od:
RETURN(true);
end:


VarTest:=proc(AA, Vars,S)
#option trace;
local B,bb,ii;
global tord,dd,n:
    B:=Basis(AA,tdeg(op(Vars)));
    B:=subs({seq(Vars[i]=0,i=S+1..nops(Vars))},B);
    B:=sort(B, proc (a, b) options operator, arrow; TestOrder(a, b, tord) end proc);
    for bb in B do 
	   for ii from n-dd to 1 by -1 while degree(bb,Vars[S])>0 do
        #for ii from nops(Vars) to 1 by -1 while degree(bb,Vars[S])>0 do
	       if not IdealMembership(bb*Vars[ii]/Vars[S], <op(B)>) then
		      RETURN(false,Vars[S],Vars[ii]);
	       fi;
	   od;
    od:
RETURN(true,1);
end:

class:=proc(F,Vars)
local i, FF, L,f;
FF:=F:
L:=[seq(0,i=1..nops(Vars))];
for i from nops(Vars) to 1 by -1 do
    for f in FF do
        if Vars[i] in indets(f) then
            L[i]:=L[i]+1;
        fi:
    od;
    FF:=subs(Vars[i]=0,FF);
    FF:=remove(has,FF,0);
od:
RETURN([seq(L[nops(L)-i+1],i=1..nops(L))]);
end:

Reg:=proc(F)
local J,s;
global tord,Vars,v,dd,n;
       J:=<op(F)>:
       if HilbertDimension(J,{op(Vars[1..n-dd])})=0 then
       	  RETURN(degree(HilbertSeries(F,s),s));
       else
	       RETURN(0);
       fi:

end:




LinearChange:=proc(Id,VARS)
#option trace;
local B,mm,LL,bb,ii,NewId,TemId,SS,hh,vars,chen,J,A,Ideal,WW,Wt,A1,A2,A3,firsttime,firstbytes,secondtime,secondbytes,ListVar,b,A4,a,aa,tt,N,flag,i,U,AA,T,reg,d,FLAG,u;
global tord,Vars,v,dd,n;
firsttime,firstbytes:=kernelopts(cputime,bytesused);
Vars:=VARS:
n:=nops(VARS):
ListVar:=NULL:
tord:=tdeg(op(Vars)):
NewId:=Basis(Id,tord):
Ideal:=NewId:
B:=LM(NewId);
dd:=HilbertDimension(<op(Ideal)>);
SS:=VarTest(B,Vars,nops(Vars)):
J:=NewId:
chen:=NULL:
while SS[1]=false do 
	   hh:=SS[2]+rand(-1..1)()*SS[3];
	   J:=NewId:
	   NewId:=Basis(subs(SS[2]=hh,NewId),tord);
	   B:=LM(NewId);
	   WW:=SS[3]:
	   Wt:=SS[2]:
	   SS:=VarTest(B,Vars,nops(Vars)):
	   if SS[1]=true then
	       chen:=chen,Wt=hh; 
	       Ideal:=subs(Wt=hh,Ideal);
	       NewId:=Ideal:
	   elif SS[1]=false and SS[3]<>WW then 
	       chen:=chen,Wt=hh;
	       Ideal:=subs(Wt=hh,Ideal);
	       NewId:=Ideal:
	   else
	       NewId:=J:
	   fi: 
od:
vars:=Vars:
for ii from 1 to dd do 
	   ListVar:=vars[-1],ListVar:
	   NewId:=subs(vars[-1]=0,NewId);
	   vars:=vars[1..-2];
	   tord:=tdeg(op(vars)):
	   NewId:=Basis(NewId,tord):
	   B:=LM(NewId);
	   for b in B while nops(vars)-dd+1>0 do 
	       if nops(indets(b))=1 and indets(b) subset {seq(VARS[i],i=nops(vars)-dd+1..nops(vars))} then 
		      tord:=tdeg(op(vars)):
		      NewId:=Basis(NewId,tord):
		      B:=LM(NewId); 
	       fi;
	   od;
	   SS:=VarTest(B,vars,nops(vars)):
	   J:=NewId:
	   while SS[1]=false do
	       hh:=SS[2]+rand(-1..1)()*SS[3];
	       J:=NewId:
	       NewId:=Basis(subs(SS[2]=hh,NewId),tord);
	       B:=LM(NewId);
	       WW:=SS[3]:
	       Wt:=SS[2]:
	       SS:=VarTest(B,vars,nops(vars)):
	       if SS[1]=true then
		      chen:=chen,Wt=hh;
		      Ideal:=subs(Wt=hh,Ideal);
	       elif SS[1]=false and SS[3]<>WW then
		      chen:=chen,Wt=hh;
		      Ideal:=subs(Wt=hh,Ideal);
	       else
		      NewId:=J:
	       fi: 
       od:	
od:

Vars:=[op(vars),ListVar]:
tord:=tdeg(op(Vars)):
AA:=Basis(Ideal,tord):
A:=LM(AA):print(A);
reg:=max(seq(degree(a),a=A)):
T:=Array(1..n,1..n,(i,j)->t[i,j]):
flag:=false:
for d from min(seq(degree(A[i]),i=1..nops(A))) to max(seq(degree(A[i],Vars[n-dd]),i=1..nops(A)))-1 while flag=false do 
    lprint("We check now the degree",d);
    U:=subs(seq(Vars[n-dd+i]=Vars[n-dd+i]-add(T[i,j]*Vars[j],j=1..n-dd),i=1..dd),MulList(AA,Vars,d));
    U:=expand(U):
    U:=subs(seq(Vars[i]=0,i=n-dd+1..n),U):
    U:=LM(InterReduce(U,tord)):
    FLAG:=false:
    for u in U while FLAG=false do
      if indets(u)={Vars[n-dd]} then
         FLAG:=true:lprint("We found a smaller reduction number:",degree(u)-1);
      fi:
    od:
    if FLAG=true then
        flag:=true:
    fi:
od:

T:=Array(1..n,1..n,(i,j)->t[i,j]):
U:=subs(seq(Vars[n-dd+i]=Vars[n-dd+i]-add(T[i,j]*Vars[j],j=1..n-dd),i=1..dd),AA);
U:=subs(seq(Vars[i]=0,i=n-dd+1..n),U):
U:=expand(U):print(U);
U:=DISPG(U,tdeg(seq(seq(T[i,j],i=1..n),j=1..n)),tord):print(seq(LM(a[2]),a=U));
a:=max(seq(Reg(LM(a[2])),a=U)):




secondtime,secondbytes:=kernelopts(cputime,bytesused);
A1:=Nested(A,Vars);
A2:=dstable2(A,Vars):
A4:=Weakdstable(A,Vars):
A3:=Borel(A,Vars):
printf("\n",""):
printf("\n",""):
printf("%1s %1s %1s %1s %1s:\n",Some,information,about,the,computation):
printf("\n",""):
printf("%-1s %1s %1s %1s:         %3a %3a\n",The, cpu, time, is,secondtime-firsttime,(sec)):
printf("%-1s %1s %1s:         %3a %3a\n",The,used,memory,secondbytes-firstbytes,(bytes)):
printf("%1s %1s %1s:%16a\n",List,of,variables,Vars):
printf("%1s :%15a\n",Dimension,dd):
printf("%1s %1s:%9a\n",Reduction,number,max(seq(degree(A[i],Vars[nops(Vars)-dd]),i=1..nops(A)))-1):
printf("%1s %1s %1s :%5a\n",Big, Reduction,number,a):
printf("%1s :%14s\n",WeakDstablity,A4):
printf("%1s %1s:%12s\n",Delta,regularity,A1):
#printf("%1s :%18s\n",dstable,dstable(A,Vars)):
printf("%1s :%18s\n",Dstablity,A2):
printf("%1s %1s:%17s\n",Borel,fixed,A3):
printf("%1s %1s %1s %1s:%1a\n",Used,change,of,variables,[chen]):
printf("\n",""):
printf("\n",""):
RETURN(); 
end:














print("####################Green");
F := [x[3]*x[1],x[2]^2+ x[2]*x[1], x[1]^2]:
LinearChange(F, [x[1],x[2],x[3]]):