############################################################################################################################################
################################################################################################################################################
################################################################################################################################################
################################The following sub-algorithms use in the PFGLM algorithm and PGBMain (KSW) algorithm#############################
################################################################################################################################################
################################################################################################################################################
################################################################################################################################################
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 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:
###################################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;
#N:=LL;
#W := M;
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;
#L:=map(u->LeadingMonomial(u,T),F);
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 #print(11111,[N,op(F[i+1..-1])]);
        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;
LinearAlgebra; 
#print(E, f, R);
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); #print(E,EG,Rr);
ns, rv := NormalSet(EG, Rr); 
Linear*Algebra; 
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; 
#option trace;
E := [op(EE)]; 
Rr:=tdeg(op(R));
LE:=LeadingMonomial(E,Rr);
V:=Groebner:-MaximalIndependentSet(LE,{op(R)}); 
d := nops(V); 
abar := [seq(i, i = 1 .. d)]; 
EV := seq(subs(V[i] = abar[i], E), i = 1 .. nops(V)); 
spE := Groebner:-Basis([EV][nops([EV])], Rr); 
#if `minus`({op(R)}, V) = IdealInfo[Variables](`<,>`(op(spE))) then 
   if IsZeroDimensional(`<,>`(op(spE))) then 
      if Zdim(spE, f, R) then 
         RETURN(true); 
      end if; 
   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, Inf) then 

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 #print(icheck); 
#   RETURN(false); 
#end if; #print(E, [f], R);
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(new corrected version)#################################
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; #print(N, M, R,salaaaaamconsist1111111111);
flag := consist1(N, M, R);  #print(byeeeeeconsist1111111111);
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; #print(whilewhilewhilewhilewhile);
   end do; #print(endendendendendendendend);
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)]; #print(seq(Variables(F[i], p), i = 1 .. nops(F)),ffff);
Ps := `union`(seq(Variables(F[i], p), i = 1 .. nops(F))); 
RETURN(Ps); 
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:
##################################################################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;  
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:
#######################################################Buchberger algorithm#################
GBasis:=proc(L,tord)
local cont,T,p,Polys,n,G,B,s,r,t1,t2,b1,b2:
global numzero;
cont:=0:
T:=tord:
Polys:=L:
n:=nops(L):
G:=[seq(i,i=1..n)]:
B:=[seq(seq([i,j],j=i+1..n),i=1..n)]:
while nops(B)<>0 do
    p:=B[1]:
    B:=B[2..-1]:
    s:=SPolynomial(Polys[p[1]],Polys[p[2]],T);
    r:=NormalForm(s,[seq(Polys[i], i in G)],T);
    if r<>0 then
        Polys:=[op(Polys),r]:
        B:=[op(B),seq([i,nops(Polys)],i in G)]:
        G:=[op(G),nops(Polys)]:
    else
        cont:=cont+1:
        numzero:=numzero+1;
    fi:
od:
G:=[seq(Polys[i], i in G)]:
RETURN(G);
end:


#############################################################PGBmain with Canspec and consist checks(without branch [1])######################KAPUR
pgbmainc := proc (EE, N, F, R, T) 
local 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:=GBasis([op(F),op(CC[2])],prod(T,R));
   G:=InterReduce(BGB,prod(T,R));
   
end if; 
#if G = [1] then 
#if evalb(1 in G) then
#   RETURN({[CC[2], factor(CC[3]), [1]]}); 
#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)):