with(CodeTools):with(LinearAlgebra):with(Groebner): with(Algebraic): with(PolynomialIdeals):

#############################################################Convert a list of Polynomials in to the linear cases polynomial###################
Homogenizz := proc (FF, T::ShortMonomialOrder, prefix::symbol)
local SS,mp,n, k, Mons, Cofs, Cons, p, M, m, AL, S, V,F;
global VV;
#option trace;
F:=expand(FF);
n := nops(F);
Mons := table(); Cofs := table(); Cons := table(sparse);
V := [op(T)];
for k to n do
   Cofs[k] := [coeffs(F[k], V, 'M')];
   M := [M];
   if member(1, M, 'p') then
      M := subsop(p = 0, M);
      Cons[k] := Cofs[k][p];
   end if;
   Mons[k] := M[];
end do;
M := sort([(`minus`({entries(Mons, 'nolist')}, {0}))[]], proc (a, b) options operator, arrow; Groebner:-TestOrder(b, a, T) end proc);
m := nops(M);
S := `~`[`=`](M, [A || (1 .. m)]);
AL := A || (m+1);
VV:=[seq(A || i, i = 1 .. m+1)];
SS:=[op(S),1=AL];
mp := map(rhs = lhs, SS);
RETURN([seq(`+`(`~`[`*`](Cofs[k], subs(S, [Mons[k]]))[])+Cons[k]*AL, k = 1 .. n)],mp):
#[seq(`+`(`~`[`*`](Cofs[k], subs(S, [Mons[k]]))[])+Cons[k]*AL, k = 1 .. n)];
end proc:
####################################################################################
lt := proc (f, T)
RETURN(LeadingTerm(f, T)[1]*LeadingTerm(f, T)[2]) :
end proc:
##################################################################################
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:
######################################################################################
sortlist := proc (L, T)
local LL, LLL;
LL := [op({seq(simplify(expand(L[i]/LeadingCoefficient(L[i], T))), i = 1 .. nops(L))})];
LLL := sort(LL, proc (b, a) options operator, arrow; TestOrder(a, b, T) end proc);
RETURN(LLL) :
end proc:
########################################################################
coefff := proc (F, flm, T)
local LLL, LLLL, l, L, LL, s, LV, cs, i, j;
LLL := NULL;
for i to nops(F) do
   l[i] := NULL;
   L[i] := FL(F[i], T);
   LL[i] := [[seq(LeadingCoefficient(L[i][n], T), n = 1 .. nops(L[i]))], [seq(LeadingMonomial(L[i][n], T), n = 1 .. nops(L[i]))]];
   for j to nops(flm) do
      s := flm[j];
      cs := member(s, LL[i][2], 'q');
      if cs then
         l[i] := l[i], LL[i][1][q] :
      else
         l[i] := l[i], 0 :
      end if
   end do;
   LLL := LLL, [l[i]] :
end do;
LLLL := [LLL];
LV := convert(LLLL, Matrix);
RETURN(LV):
end proc:
###################################################################
produceM := proc (L, G, T)
local H, don, MonH, Menha, SortMenha, xbeta, n, M;
H := L;
don := [seq(LeadingMonomial(H[i], T), i = 1 .. nops(H))];
don := {op(don)};
MonH := {seq(op(FL(H[i], T)), i = 1 .. nops(H))};
MonH := sortlist([op(MonH)], T);
while don <> {op(MonH)} do
   Menha := `minus`({op(MonH)}, don);
   SortMenha := sortlist([op(Menha)], T);
   xbeta := SortMenha[1];
   don := `union`(don, {xbeta});
   for n to nops(G) do
      if divide(xbeta, LeadingMonomial(G[n], T), 'q') then
         H := [op(H), simplify(expand(q*G[n]))];
         MonH := {seq(op(FL(H[i], T)), i = 1 .. nops(H))};
         MonH := sortlist([op(MonH)], T);
      end if;
   end do;
end do;
M := coefff(H, MonH, T);
RETURN(H, M, MonH) ;
end proc:
#################################################################
sortdeg := proc (P, Q)
RETURN(evalb(P[2] <= Q[2]));
end proc:
#####################################################################
SelectBp := proc (Bpairs)
local Bp, i, B, n;
B := Bpairs;
n := B[1][2];
Bp := B[1];
for i from 2 to nops(B) while B[i][2] = n do
   Bp := Bp, B[i];
end do;
RETURN([Bp]);
end proc:
############################################################################
Listpoly := proc (M, Mon)
local ML, Listpol;
ML := convert(M, list, nested);
Listpol := [seq(sum([seq(ML[j][i]*Mon[i], i = 1 .. nops(Mon))][ii], ii = 1 .. nops(Mon)), j = 1 .. nops(ML))];
RETURN(subs(0 = NULL, Listpol));
end proc:
######################################################################
newpolys := proc (NList, MList, T)
local LM, newpol, j;
LM := `<,>`(seq(LeadingMonomial(MList[i], T), i = 1 .. nops(MList)));
newpol := NULL;
for j to nops(NList) do
   if IdealMembership(LeadingMonomial(NList[j], T), LM) = false then
      newpol := newpol, NList[j];
   end if;
end do;
RETURN(newpol);
end proc:
########################################################################
F4Basis := proc (F, T)
local n, NL, Npositive, G, t, B, M, L, L1, L2, N, Bp;
G := F;
t := nops(F);
B := [seq(seq([[i, j], degree(lcm(LeadingMonomial(F[i], T), LeadingMonomial(F[j], T)))], j = i+1 .. t), i = 1 .. t-1)];
while B <> [] do
   B := sort(B, sortdeg);
   Bp := SelectBp(B);
   B := B[nops(Bp)+1 .. -1];
   L1 := {seq(simplify(expand(lcm(LeadingMonomial(G[Bp[i][1][1]], T), LeadingMonomial(G[Bp[i][1][2]], T))*G[Bp[i][1][1]]/lt(G[Bp[i][1][1]], T))), i = 1 .. nops(Bp))};
   L2 := {seq(simplify(expand(lcm(LeadingMonomial(G[Bp[i][1][1]], T), LeadingMonomial(G[Bp[i][1][2]], T))*G[Bp[i][1][2]]/lt(G[Bp[i][1][2]], T))), i = 1 .. nops(Bp))};
   L := `union`(L1, L2);
   M := produceM(L, G, T);
   N := ReducedRowEchelonForm(M[2]);
   NL := Listpoly(N, M[3]);
   Npositive := [newpolys(NL, M[1], T)];
   for n to nops(Npositive) do
      t := t+1;
      G := [op(G), Npositive[n]];
      B := [op(B), seq([[i, t], degree(lcm(LeadingMonomial(G[i], T), LeadingMonomial(G[t], T)))], i = 1 .. t-1)];
   end do;
end do;
RETURN(G);
end proc:
#############################################################################
##############################################################################
#######################################################################################
#######################################################################################
############################################compute the normalform and quotient of f on F##########################
normalform := proc (f, F, T)
local nf1, DN, Lc, NUM, flc, Q;
if F = [] then
   nf1 := Groebner:-NormalForm(f, F, T);
   Q := []
else
   DN := [seq(denom(F[i]), i = 1 .. nops(F))];
   Lc := lcm(op(DN));
   NUM := simplify(expand(Lc*F));
   flc := simplify(expand(Lc*f));
   nf1 := Groebner:-NormalForm(flc, NUM, T, 'Q');
   nf1 := simplify(expand(nf1/Lc));
end if;
RETURN([nf1, Q]);
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:
########################################################input: poly f and null condition set N and monomial order on variables T###################################
########################################################output: substitute elements of N in f######################################################
subss := proc (f, N, T)
local M, S, g;
#option trace;
#print(input,f, N, T);
M := Poly(f, [op(T)]);
S := subs(seq(N[i] = 0, i = 1 .. nops(N)), M);
g := add(s, `in`(s, S)); #print(output,g);
RETURN(g);
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:
##############################################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:

##################################################################Return true if a is constant#########################################
isconstant := proc (a)
RETURN(type(a, constant));
end proc:
##################################Return true if there is a constant coefficient in a list of terms#########################
isthereconstant := proc (A)
local i;
for i to nops(A) do
   if isconstant(A[i]) then
      RETURN(true);
   end if;
end do;
RETURN(false);
end proc:
#######################################################################LDS algorithm################################################################
LinGB13 := proc (N, W, F, f, R, T)
global nfp, LIST, contor;
#option trace;
LIST := NULL; contor := 1;
RETURN(IsPlinMH(N, W, F, f, R, T));
end proc:

IsPlinMH := proc (NNN, WW, FF, ff, R, T)
local fp, CNS, flag2, f, TERMS, COEFS, flag, LISTT, i, M,F,F1,F2,Sys,h;
global LIST, contor;
#option trace;
#print(NNN, WW, FF, ff, R, T);
f := Groebner:-NormalForm(numer(ff), simplify(expand(denom(ff)*NNN)), R);
f := simplify(expand(f/denom(ff))); #print(f, FF, T);

F:=FF;
F1:=simplify(expand([seq(Groebner:-NormalForm(numer(F[i]),NNN,R),i=1..nops(F))]));
F1:=subs(0=XX,F1);
F2:=[seq(denom(F[i]),i=1..nops(F))];
F := [seq(simplify(expand(F1[i]/F2[i])), i = 1 .. nops(F))];

fp := simplify(expand(normalform(f, F, T)));
if fp[1] = 0 then
   LIST := LIST, [[true, fp[2],fp[1]], numer(NNN), numer(WW)];
else
   M := Poly(fp[1], [op(T)]);
   TERMS := sort(M, proc (b, a) options operator, arrow; TestOrder(a, b, T) end proc);
   COEFS := LeadingCoefficient(TERMS, T);
   for i to nops(COEFS) while not type(COEFS[i], constant) do
      LIST := LIST, [[false, fp[2],fp[1]], numer([op(NNN), seq(COEFS[j], j = 1 .. i-1)]), numer([op(WW), COEFS[i]])];
   end do;
   if i <= nops(COEFS) then
      if type(COEFS[i], constant) then
         LIST := LIST, [[false, fp[2],fp[1]], numer([op(NNN), seq(COEFS[j], j = 1 .. i-1)]), numer(WW)];
      end if;
   end if;
   flag := isthereconstant(COEFS);
   if flag = false then
      LIST := LIST, [[true, fp[2],fp[1]], numer([op(NNN), op(COEFS)]), WW];
   end if;
end if;
LIST := [LIST];
LISTT := NULL;
for i to nops(LIST) do
   #CNS := ccc(LIST[i][2], LIST[i][3], R);
  # print(LIST[i][2], LIST[i][3], R,111111111111111);
   CNS := canspec(LIST[i][2], LIST[i][3], R);
   flag2 := CNS[1];
   if flag2 = true then
      LISTT := LISTT, [LIST[i][1], CNS[2], CNS[3]];
   end if;
end do;

LIST:=[LISTT];#print(LIST);
Sys:=NULL;
for i to nops(LIST) do
   h := subss(LIST[i][1][3], LIST[i][2], T);
   Sys := Sys, [[LIST[i][1][1], LIST[i][1][2], h],LIST[i][2], LIST[i][3]];
end do;
RETURN(Sys);
end proc:
####################################################convert a poly f in to its collect term w.r.t. its variables Var###############
Poly := proc (f, Var)
local g, fl;
g := collect(f, Var);
fl := FL(g, plex(op(Var)));
RETURN(fl);
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:
#####################################################(new version)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:
############################################GGES algorithm, its output is {(null,notnull,list of general polynomials)}####################################
GGES := proc (FF,TT)
local HM,F, P, A, i, a, g, PL, Syss, T, G,R,var;
global VV,HM2;
#option trace;
Syss := NULL;
#F := SYS(M);
HM:=Homogenizz(FF,TT);
HM2:=HM[2];
F:=HM[1];#print(F);
T := plex(op(VV));
var := vars(F, T);
R := plex(op(var)); #print(R,T);
A := [[[], [], [], F[1], F]];
while A <> [] do
   a := A[1];
   A := subs(A[1] = NULL, A);
   if a[5] = [] then
      #Syss := Syss, [a[1], a[2], a[3]]; #print(subs(seq(t, t in [HM][2]), a[3]),444444444);
      Syss := Syss, [a[1], a[2], subs(seq(t, t in [HM][2]), a[3])];
   else
      #G := subs(a[5][1] = NULL, a[5]); ##sabz example#
      G := subsop(1 = NULL, a[5]);
      if G <> [] then
         g := G[1];
      end if; #print(a[1], a[2], a[3], a[4], R, T,1111111111);
      PL := [LinGB13(a[1], a[2], a[3], a[4], R, T)];
      for i to nops(PL) do
         P := PL[i];
         if P[1][1] = true then
            A := [op(A), [P[2], P[3], subs(0 = NULL, a[3]), g, G]];
         else
            A := [op(A), [P[2], P[3], subs(0 = NULL, [op(a[3]), P[1][3]]), g, G]];
         end if;
      end do;
   end if;
end do;
RETURN(Syss);
end proc:
##############################################################
NormalList := proc (F, N, R)
local GG, FF;
GG := NULL;
FF := [seq(normalform(op(nrm([f])), N, R)[1], f in F)];
RETURN(FF);
end proc:
######################################################################
HasSubexpression := proc (s::algebraic, e) local S; has(algsubs(s = S, e), S) end proc;
########################################################################
NormalList1 := proc (GS, R) local GG, FF, funcc, gs, func2; GG := NULL;
for gs in GS do
   funcc := proc (S) options operator, arrow; map(has, S, gs[1]) end proc;
   if gs[1] = [] then
      GG := GG, [gs[1], gs[2], [seq(op(nrm([f])), f in gs[3])]];
   elif true in funcc(gs[3]) then
      GG := GG, [gs[1], gs[2], [seq(normalform(op(nrm([f])), gs[1], R)[1], f in gs[3])]];
   elif ormap(HasSubexpression, gs[1], gs[3]) then
      GG := GG, [gs[1], gs[2], [seq(normalform(op(nrm([f])), gs[1], R)[1], f in gs[3])]];
   else
      GG := GG, [gs[1], gs[2], [seq(op(nrm([f])), f in gs[3])]];
  end if;
end do;
RETURN([GG]);
end proc:

############################################GGES algorithm, its output is {(null,notnull,list of general polynomials)}####################################
GGESNW := proc (N,W,FF,TT)
local AN,ANN,HM,F, P, A, i, a, g, PL, Syss, T, G,R,var,Mons;
global VV,HM2;
#option trace;
Syss := NULL;
#F := SYS(M);
HM:=Homogenizz(FF,TT);#print(HM[1],N,W,FF);
HM2:=HM[2];
F:=HM[1];#print(F,HM);
T := plex(op(VV));
var := vars(F, T); #print(vvvvvvvv,{op(var),op(vars(N, T)),op(vars(W, T))});
var:={op(var),op(vars(N, T)),op(vars(W, T))};
R := plex(op(var)); #print(R,T);
A := [[N, W, [], F[1], F]];
while A <> [] do
   a := A[1];
   A := subs(A[1] = NULL, A);
   if a[5] = [] then
      #Syss := Syss, [a[1], a[2], a[3]]; #print(subs(seq(t, t in [HM][2]), a[3]),444444444);
      Syss := Syss, [a[1], a[2], subs(seq(t, t in [HM][2]), a[3])];
   else
      #G := subs(a[5][1] = NULL, a[5]); ##sabz example#
      G := subsop(1 = NULL, a[5]);
      if G <> [] then
         g := G[1];
      end if; #print(a[1], a[2], a[3], a[4], R, T,1111111111);
      PL := [LinGB13(a[1], a[2], a[3], a[4], R, T)];
      for i to nops(PL) do
         P := PL[i];
         if P[1][1] = true then
            A := [op(A), [P[2], P[3], subs(0 = NULL, a[3]), g, G]];
         else
            A := [op(A), [P[2], P[3], subs(0 = NULL, [op(a[3]), P[1][3]]), g, G]];
         end if;
      end do;
   end if;
end do;
Mons:=map(rhs, HM2):
AN:=[Syss];#print(beforeeeeee,AN);
ANN:=NormalList1(AN,R);
#ANN:=[seq([AN[i][1], AN[i][2], NormalList(AN[i][3], AN[i][1], R)], i = 1 .. nops(AN))];#print(afterrrrrr,ANN);
#RETURN([Syss],Mons);
RETURN(ANN,Mons);
end proc:
########################################gives us the number of zero in a list###################################
zer := proc (L)
local i, j;
i := 0;
for j to nops(L) do
   if evalb(L[j] = 0) then
      i := i+1;
   end if;
end do;
RETURN(i);
end proc:
########################################its input is a list of lists and sort this list such as upper triangular form######################
szer := proc (LL1, LL2)
local p, q, flag1, flag2, n1, L1, L2, n2;
L1 := LL1;
L2 := LL2;
if member(0, L1) then
   flag1 := member(0, L1, 'q');
else
   q := 1000;
end if;
if member(0, L2) then
   flag2 := member(0, L2, 'p');
else
   p := 1000;
end if;
if q = p then
   n1 := zer(L1);
   n2 := zer(L2);
   if n1 <> n2 then
      RETURN(evalb(n1 < n2));
   else
      L1 := subsop(1 = NULL, L1);
      L2 := subsop(1 = NULL, L2);
      if L1 <> [] then
         szer(L1, L2);
      else
         RETURN(true);
      end if;
   end if;
else
RETURN(evalb(p < q));
end if;
end proc:

################################Forming a coefficient matrix from the list of polynomial F w.r.t. a monomial ordering T ###############
coefmat := proc (F, T) local flm, LLL, LLLL, l, L, LL, s, LV, cs, i, j;
flm := FLm(F, T);
LLL := NULL;
for i to nops(F) do
   l[i] := NULL;
   L[i] := FL(F[i], T);
   LL[i] := [[seq(LeadingCoefficient(L[i][n], T), n = 1 .. nops(L[i]))], [seq(LeadingMonomial(L[i][n], T), n = 1 .. nops(L[i]))]];
   for j to nops(flm) do
      s := flm[j];
      cs := member(s, LL[i][2], 'q');
      if cs then
         l[i] := l[i], LL[i][1][q];
      else
         l[i] := l[i], 0;
      end if;
   end do;
   LLL := LLL, [l[i]];
end do;
LLLL := [LLL];
LV := convert(LLLL, Matrix);
RETURN(LV);
end proc:

####################################################Input: A sqree matrix $A$, Output: A linear system with coefficient matrix $A$########################
SYS := proc (A)

local r, m, S, VS, ADD;
global VV;
#option trace;
r := RowDimension(A);
m := ColumnDimension(A);
VV := [seq(x[i], i = 1 .. m)];
S := [seq(convert(A[i], list), i = 1 .. r)];
VS := [seq(`~`[`*`](S[i], VV), i = 1 .. r)];
ADD := [seq(add(VS[i][j], j = 1 .. nops(VS[i])), i = 1 .. r)];
RETURN(ADD);
end proc:
##################################################GES algorithm, its output is {(null,notnull,list of polynomials)}####################################
GES := proc (M)
local F, P, A, i, a, g, PL, Syss, T, G,R,var;
global VV;
#option trace;
Syss := NULL;
F := SYS(M);
T := plex(op(VV));
var := vars(F, T);
R := plex(op(var)); #print(R,T);
A := [[[], [], [], F[1], F]];
while A <> [] do
   a := A[1];
   A := subs(A[1] = NULL, A);
   if a[5] = [] then
      Syss := Syss, [a[1], a[2], a[3]];
   else
      #G := subs(a[5][1] = NULL, a[5]); ##sabz example#
      G := subsop(1 = NULL, a[5]);
      if G <> [] then
         g := G[1];
      end if;
      PL := [LinGB13(a[1], a[2], a[3], a[4], R, T)];
      for i to nops(PL) do
         P := PL[i];
         if P[1][1] = true then
            A := [op(A), [P[2], P[3], subs(0 = NULL, a[3]), g, G]];
         else
            A := [op(A), [P[2], P[3], subs(0 = NULL, [op(a[3]), P[1][3]]), g, G]];
         end if;
      end do;
   end if;
end do;
RETURN(Syss);
end proc:
M := GenerateMatrix([a*x-b*y,c*x+d*z], [x, y, z])[1]:
GES(M): AaA := [%]: nops(AaA):
###################The following code is also GES algorithm but its input is a list of parametric linear polynomial and its output is matrix form#############
#############################Parametric Gaussian elimination (Matrix form)(our algorithm based on above "GES" algo)##########################
#################################################The output is {(null,non-null,matrix representation)}##############################
PGES1 := proc (F, R, T)
local TT, PG, PGBB, DPGB, div, t, L, LL, LLL, l, dpgb, j, FF, cs, LV, LLLL, flm, s, PGB, cmat, M, rd, cd, H, HH, CMat, Cmat, i;
#option trace;
M := coefmat(F, T); #print(M);
rd := RowDimension(M);
cd := ColumnDimension(M);
PGBB := [GES(M)];
DPGB := NULL;
TT := plex(op(VV)); #print(TT,VV,111111111111);
flm := VV;
for j to nops(PGBB) do #print(PGBB,1111111,R);
   div := [seq(NormalForm(numer(PGBB[j][3][i]), PGBB[j][1], R)/denom(PGBB[j][3][i]), i = 1 .. nops(PGBB[j][3]))];
   FF := div;
   LLL := NULL;
   for i to nops(FF) do
      l[i] := NULL;
      L[i] := FL(FF[i], TT);
      LL[i] := [[seq(LeadingCoefficient(L[i][n], TT), n = 1 .. nops(L[i]))], [seq(LeadingMonomial(L[i][n], TT), n = 1 .. nops(L[i]))]];
      for t to nops(flm) do
         s := flm[t];
         cs := member(s, LL[i][2], 'q');
         if cs then
            l[i] := l[i], LL[i][1][q];
         else
            l[i] := l[i], 0;
         end if;
      end do;
      LLL := LLL, [l[i]];
   end do;
   LLLL := [LLL];
   LV := convert(LLLL, Matrix);
   cmat := LV;
   Cmat := Matrix(rd, cd, cmat, fill = 0);
   HH := NULL;
   for i to RowDimension(Cmat) do
      HH := HH, convert(Cmat[i], list);
   end do;
   H := sort([HH], szer); #print(H);
   CMat := convert(H, Matrix);
   dpgb := [PGBB[j][1], PGBB[j][2], CMat];
   DPGB := DPGB, dpgb;
end do;
RETURN(DPGB);
end proc:
PGES1([a*x+b*y,z+x+c*y,(a-b)*x+d*z],tdeg(a,b,c,d),tdeg(x,y,z)):

#################################################Sort the monimoals appeared in polynomials of F w.r.t. a monomial order T ##########################
FLm := proc (F, T)
local L, LL, LLL;
L := [op({seq(op(FL(F[i], T)), i = 1 .. nops(F))})];
LL := [op({seq(simplify(expand(L[i]/LeadingCoefficient(L[i], T))), i = 1 .. nops(L))})];
LLL := sort(LL, proc (b, a) options operator, arrow; TestOrder(a, b, T);
end proc);
RETURN(LLL);
end proc:
################################Forming a coefficient matrix from the list of polynomial F w.r.t. a monomial ordering T ###############
coefmat := proc (F, T) local flm, LLL, LLLL, l, L, LL, s, LV, cs, i, j;
flm := FLm(F, T);
LLL := NULL;
for i to nops(F) do
   l[i] := NULL;
   L[i] := FL(F[i], T);
   LL[i] := [[seq(LeadingCoefficient(L[i][n], T), n = 1 .. nops(L[i]))], [seq(LeadingMonomial(L[i][n], T), n = 1 .. nops(L[i]))]];
   for j to nops(flm) do
      s := flm[j];
      cs := member(s, LL[i][2], 'q');
      if cs then
         l[i] := l[i], LL[i][1][q];
      else
         l[i] := l[i], 0;
      end if;
   end do;
   LLL := LLL, [l[i]];
end do;
LLLL := [LLL];
LV := convert(LLLL, Matrix);
RETURN(LV);
end proc:
###################The following code is also GES algorithm but its input is a list of parametric linear polynomial and its output is matrix form#############
#############################Parametric Gaussian elimination (Matrix form)(our algorithm based on above "GES" algo)##########################
#################################################The output is {(null,non-null,matrix representation)}##############################
MGES := proc (F, R, T)
local TT, PG, PGBB, DPGB, div, t, L, LL, LLL, l, dpgb, j, FF, cs, LV, LLLL, flm, s, PGB, cmat, M, rd, cd, H, HH, CMat, Cmat, i;
option trace;
M := coefmat(F, T); #print(M);
rd := RowDimension(M);
cd := ColumnDimension(M);
PGBB := [GGES(F,T)];
DPGB := NULL;
#TT := plex(op(VV));flm := VV;
TT := plex(op(T));
flm:=[op(T)];
for j to nops(PGBB) do #print(PGBB,1111111,R);
   div := [seq(NormalForm(numer(PGBB[j][3][i]), PGBB[j][1], R)/denom(PGBB[j][3][i]), i = 1 .. nops(PGBB[j][3]))];
   FF := div;
   LLL := NULL;
   for i to nops(FF) do
      l[i] := NULL; #print(FF[i], TT);
      L[i] := FL(FF[i], TT);
      LL[i] := [[seq(LeadingCoefficient(L[i][n], TT), n = 1 .. nops(L[i]))], [seq(LeadingMonomial(L[i][n], TT), n = 1 .. nops(L[i]))]];
      for t to nops(flm) do
         s := flm[t];
         cs := member(s, LL[i][2], 'q');
         if cs then
            l[i] := l[i], LL[i][1][q];
         else
            l[i] := l[i], 0;
         end if;
      end do;
      LLL := LLL, [l[i]];
   end do;
   LLLL := [LLL];
   LV := convert(LLLL, Matrix);
   cmat := LV;
   Cmat := Matrix(rd, cd, cmat, fill = 0);
   HH := NULL;
   for i to RowDimension(Cmat) do
      HH := HH, convert(Cmat[i], list);
   end do;
   H := sort([HH], szer);
   CMat := convert(H, Matrix);
   dpgb := [PGBB[j][1], PGBB[j][2], CMat];
   DPGB := DPGB, dpgb;
end do;
RETURN(DPGB);
end proc:
############################################3
#############################################
#####################################################compute the normalform of elements of F w.r.t. null set N and order R#########
NormalList := proc (F, N, R)
local GG, FF;
GG := NULL;
FF := [seq(normalform(op(nrm([f])), N, R)[1], f in F)];
RETURN(FF):
end proc:
#####################################################L: is some lists of polynomials & B: is a set of monomials########################
#####################################################computes a matrix representation of list of L w.r.t. B############################
GM := proc (L, B)
local coeffff, T;
T := indets(B);
coeffff := proc (P, T, t)
local C, H, i, k;
C := [coeffs(P, T, 'h')]; H := [h];
k := 0;
for i to nops(H) do
   if H[i] = t then
      k := C[i]:
   end if:
end do;
k:
end proc;
map(proc (p) options operator, arrow; Matrix(nops(p), nops(B), proc (i, j) options operator, arrow;
coeffff(p[i], T, B[j]):
end proc):
end proc, L):
end proc:
#######################################
#################################################Computes Grobner system of F w.r.t. RR and TT by using the following algorithm(PF4Basis)######
PF4 := proc (F, RR, TT)
local N, W, t, B;
global SYS, R, T, OUTSYS,setB ,Tord;
R := RR;
T := TT;
Tord:=TT;
N := [];
W := LeadingCoefficient(F, TT);
t := nops(F);
B := [seq(seq([[i, j], degree(lcm(LeadingMonomial(F[i], T), LeadingMonomial(F[j], T)))], j = i+1 .. t), i = 1 .. t-1)];
SYS := [N, W, F, t, B];
OUTSYS := NULL;
RETURN(PF4Basis(op(SYS))):
end proc:
#######################################decide that L is contain polynomial with degree atmost 1 or not##########################
islinear := proc (L, R)
local vars;
vars := `minus`(indets(L), {op(R)});
RETURN(andmap(type, L, Or(freeof(vars), linear(vars))));
end proc:
###################################Convert a matrix to the corresponding linear system with original variables##########################################
SYSS := proc (A, T)
#option trace;
local r, m, S, VS, ADD, VVV;
global VV;
VV := [op(T)];
r := RowDimension(A);
m := ColumnDimension(A);
if nops(VV) < m then
   VV := [op(VV), 1];
end if;
S := [seq(convert(A[i], list), i = 1 .. r)];
VS := [seq(`~`[`*`](S[i], VV), i = 1 .. r)];
ADD := [seq(add(VS[i][j], j = 1 .. nops(VS[i])), i = 1 .. r)];
RETURN(ADD);
end proc:
##################################### This is PGES1 with polynomial system presentation instead of matrix representation###################################
LineGes := proc (F, R, T)
local AB, ssyyss;
#global VV;
#option trace;
#VV := [op(T)];
AB := [PGES1(F, R, T)]; #print(111111111);
ssyyss := seq([AB[i][1], AB[i][2], SYSS(AB[i][3], T)], i = 1 .. nops(AB));
RETURN(ssyyss);
end proc:
#################################################its input is a Grobner system and term ordering on variables#################################
#################################################its output is two list of input GS s.t first list contain branches with grobner basis [1]#####
#################################################and second list is contain those branches s.t. their Grobner basis is not [1]#################
hasconstant := proc (GS, T)
local oneGB, GB, i, V, func;
oneGB := NULL;
GB := NULL;
V := {op(T)};
func := proc (S) options operator, arrow; not andmap(has, S, V) end proc:
for i to nops(GS) do
   if func(GS[i][3]) = true then
      oneGB := oneGB, [GS[i][1], GS[i][2], [1]];
   else
      GB := GB, GS[i];
   end if;
end do;
RETURN([oneGB], [GB]);
end proc:
##############this###################################Computes Grobner system of F w.r.t. RR and TT by using the following algorithm(PF4Basis)######
###########################################
#############################################################################
##############################################################################
DisjMon := proc (m1, m2, vars::(list(name))) options operator, arrow; not depends(gcd(m1, m2), vars) end proc:
firstupdate := proc (L, R, T)
local B, bb, FF, t, F, tt, tn, BB, V, numb, m, n;
numb := NULL;
F := L;
t := nops(F);
B := [seq(seq([[i, j], degree(lcm(LeadingMonomial(F[i], T), LeadingMonomial(F[j], T)))], j = i+1 .. t), i = 1 .. t-1)];
V := [op(T)];
for bb in B do
   if DisjMon(LeadingMonomial(F[bb[1][1]], T), LeadingMonomial(F[bb[1][2]], T), V) then
      B := subs(bb = NULL, B)
   end if
end do;
RETURN([F, B]);
end proc:

#########################oldversion#####
firstupdate1 := proc (L, R, T)
local B, bb, FF, t, F, tt, tn, BB, V, numb,m,n;
numb := NULL;
F := L;
tt := nops(F);
BB := [seq(seq([[i, j], degree(lcm(LeadingMonomial(F[i], T), LeadingMonomial(F[j], T)))], j = i+1 .. tt), i = 1 .. tt-1)];
for bb in BB do
    m, n := [op(Variables(LeadingMonomial(F[bb[1][1]], T), [op(R)]))], [op(Variables(LeadingMonomial(F[bb[1][2]], T), [op(R)]))];
   if divide(LeadingMonomial(F[bb[1][1]], T), LeadingMonomial(F[bb[1][2]], T)) and m=n then
      numb := numb, bb[1][1]
   elif divide(LeadingMonomial(F[bb[1][2]], T), LeadingMonomial(F[bb[1][1]], T)) and m=n then
      numb := numb, bb[1][2]
   end if;
end do;
F := subs(seq(F[num] = NULL, num in [numb]), F);
t := nops(F);
B := [seq(seq([[i, j], degree(lcm(LeadingMonomial(F[i], T), LeadingMonomial(F[j], T)))], j = i+1 .. t), i = 1 .. t-1)];
V := [op(T)];
for bb in B do
   if DisjMon(LeadingMonomial(F[bb[1][1]], T), LeadingMonomial(F[bb[1][2]], T), V) then
      B := subs(bb = NULL, B)
   end if;
end do;
RETURN([F, B]);
end proc:
#################################################Computes Grobner system of F w.r.t. RR and TT by using the following algorithm(PF4Basis)######
NEWPF4 := proc (F, RR, TT)
local AA, n, N, W, t, B, BB, twolists,FF;
global SYS, R, T, OUTSYS,setB ,Tord;
#option trace;
#print(input,newpf4);
R := RR;
T := TT;
Tord:=TT;
if islinear(F,R) then #print(F,R,T);
   RETURN(LineGes(F,R,T));
fi;
OUTSYS := NULL;#print([],[],F,TT);
BB:=GGESNW([],[],F,TT)[1]; #print(1111111111111,BB,nops(BB));
twolists:=hasconstant(BB,T);
AA:=twolists[2];print(notone,nops(AA));
OUTSYS:=OUTSYS,op(twolists[1]);#print(outsys,nops([OUTSYS]));
for n from 1 to nops(AA) do print(nops([OUTSYS]),nnnnnnnnnnnnnnnnnn,n);
     FF:=firstupdate(AA[n][3],R,T);
     t:=nops(FF[1]);
     B:=FF[2];
     SYS:=[AA[n][1],AA[n][2],FF[1],t,B];

     #t:=nops(AA[n][3]);
     #B:=[seq(seq([[i, j], degree(lcm(LeadingMonomial(AA[n][3][i], T), LeadingMonomial(AA[n][3][j], T)))], j = i+1 .. t), i = 1 .. t-1)];
     #SYS:=[AA[n][1],AA[n][2],AA[n][3],t,B];
     PF4Basis(op(SYS))
od;#print(output,newpf4);

#RETURN(OUTSYS);
end proc:
###########################################
###########################################
################################################
PF4Basis := proc (N, W, F, nopsF, Bpairs)
local n, NL, Npositive, G, t, B, M, L, L1, L2, Bp, Ges, i, NPi, j, allNPi, sys, Bsys,BJ,tt,SYSnull;
global OUTSYS, SYS, R, T, setB ,Tord;
#option trace;
#print(inpuuuuuuuuuuuuuuuuuut,pf4basis,N, W, F, nopsF, Bpairs);
B := Bpairs;
#flagnew:=false;
while [SYS] <> [] do
   sys := [SYS][1];
   SYS := [SYS][2 .. -1];
   SYS := op(SYS);
   G := sys[3];
   if sys[5] = [] then
      OUTSYS := OUTSYS, [sys[1], sys[2], sys[3]]
   else
      B := sort(sys[5], sortdeg);
      Bp := SelectBp(B); B := B[nops(Bp)+1 .. -1];
      Bsys := B; L1 := {seq(simplify(expand(lcm(LeadingMonomial(G[Bp[i][1][1]], T), LeadingMonomial(G[Bp[i][1][2]], T))*G[Bp[i][1][1]]/lt(G[Bp[i][1][1]], T))), i = 1 .. nops(Bp))};
      L2 := {seq(simplify(expand(lcm(LeadingMonomial(G[Bp[i][1][1]], T), LeadingMonomial(G[Bp[i][1][2]], T))*G[Bp[i][1][2]]/lt(G[Bp[i][1][2]], T))), i = 1 .. nops(Bp))};
      L := `union`(L1, L2);
      M := produceM(L, G, T); 
      Ges := GGESNW(sys[1], sys[2], M[1], T);  
      allNPi := NULL;
      for i to nops(Ges[1]) do
         NPi := [newpolys(Ges[1][i][3], M[1], T)];
         allNPi := allNPi, NPi;
      end do; #print(allnpi,[allNPi],3333333333333);
      for j to nops([allNPi]) do
         t := sys[4];
         G := NormalList(sys[3], Ges[1][j][1], R);
         B := Bsys;
         for n to nops([allNPi][j]) do
            t := t+1;
            G := [op(G), [allNPi][j][n]];
            BJ := [seq({op(B[i][1])}, i = 1 .. nops(B))];
            setB:=G;#print(t, [seq(i, i = 1 .. t-1)], BJ,setB);
            B := update(t, [seq(i, i = 1 .. t-1)], BJ):
         end do;
         SYS := SYS, [Ges[1][j][1], Ges[1][j][2], G, t, B]; 
         SYSnull:=NULL;
         for tt from 1 to nops([SYS]) do
            sys := [SYS][tt];
            if sys[5] = [] then
               OUTSYS := OUTSYS, [sys[1], sys[2], sys[3]];
               SYSnull:=SYSnull,sys;
            fi;
         od;
         SYS := op({SYS}minus{SYSnull});#print([SYS],555555555555555555555555);

      end do;
      seq(PF4Basis(op(syys)), syys in [SYS]);
   end if;
end do; 
 if [SYS]=[] then
      RETURN(OUTSYS);
 fi;
RETURN(OUTSYS):
end proc:
######Update Algoritm: return critical pairs only########################
update := proc (h, LL, BB)
local L1, L, L2, L3, L4, L5, i, lth, flag, j, p, E,nn,E1,E2;
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; #print(BB,bbbbb);
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;  #print(555555555555555);
if BB=[] and nops(LL)=1 then
    RETURN([{h,op(LL)}]);
fi;
E1 := [seq([{h, i}, degree(lcm(LeadingMonomial(setB[h], Tord), LeadingMonomial(setB[i], Tord)))], i in L1[nn+1 .. nops(L1)])];#print(L5,setB,Tord);
E2 := [seq([nm, degree(lcm(LeadingMonomial(setB[nm[1]], Tord), LeadingMonomial(setB[nm[2]], Tord)))], nm in [L5])];
E := [op(E1), op(E2)];
RETURN(E):
end proc:
######Update Algoritm: returns critical pairs B and set G########################
update1 := proc (h, LL, BB)
local L1, L, L2, L3, L4, L5, i, lth, flag, j, p, E,nn,E1,E2;
global setB, Tord;
#option trace;
#print(inputUPDATE,h, LL, BB);#print(setB, Tord,2222222);
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; #print(BB,bbbbb);
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;
E1 := [seq([{h, i}, degree(lcm(LeadingMonomial(setB[h], Tord), LeadingMonomial(setB[i], Tord)))], i in L1[nn+1 .. nops(L1)])];#print(L5,setB,Tord);
E2 := [seq([nm, degree(lcm(LeadingMonomial(setB[nm[1]], Tord), LeadingMonomial(setB[nm[2]], Tord)))], nm in [L5])];
E := [op(E1), op(E2)];
RETURN(E):
end proc: