#*******************************************************************************
#*******************************************************************************
#***************Kapur et al. Algorithm (PGBMAIN 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:
#################################################its input is a Grobner system(the output of pgbmain), This algo remove some additional branch#################
filter := proc (AA) 
#option trace;
local A, i; 
A := NULL; 
for i to nops(AA) do 
   if AA[i] <> [] and AA[i] <> {} and AA[i] <> {[{}, [], [1]]} and AA[i] <> {[[1], [], [1]]} then 
   
      if nops(AA)=1 then 
        A:=A,AA[1];
      else
        A := A, op(AA[i]);
      fi;
       
   end if; 
end do; 
RETURN(A); 
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:

                      ##############################GGE ALGORITM###############################
                      ###############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,WWWW;
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])); 
WWWW:=subs(0=NULL,[WWW]					);
test := true; 
t:= true; 
NN := N; 
h := product(WWWW[i], i = 1 .. nops(WWWW)); 
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:

                                 #########################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;
h := product(M[i], i = 1 .. nops(M));
if RadicalMembership(h, `<,>`(LL)) = true then
    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##################################
MDB := proc (L, T) 
local K, p, n, N, H, ll, M, i; 
N := NULL; 
H := NULL; 
K := L; 
while K <> [] do 
   ll := NULL; 
   p := K[1]; 
   for i to nops(K) do 
      if divide(LeadingMonomial(K[i], T), LeadingMonomial(p, T)) then 
         ll := ll, K[i]; 
      elif divide(LeadingMonomial(p, T), LeadingMonomial(K[i], T)) then 
         H := H, p; 
      end if; 
   end do; 
   K := [op(`minus`({op(K)}, {ll}))]; 
   N := N, p; 
end do; 
M := [op(`minus`(`union`({op(K)}, {N}), {H}))]; 
RETURN(M); 
end proc:
##########################################################Benyamin Minimal Dickson Basis######################################
MDBasis:=proc(F,T)
local g,N,L,i,flag;
N:=NULL;
for i from 1 to nops(F) do
    flag:=false;
    for g in [N,op(F[i+1..-1])] while not flag do
        if divide(LeadingMonomial(F[i],T),LeadingMonomial(g,T)) then 
            flag:=true;
        fi;
    od;
    if not flag then
        N:=N,F[i];
    fi;
od;
RETURN([N]);
end:
###################################################Compare the max of degree of factors of LC of f,g for sort a list of poly################
comp := proc (f, g) 
local maxm, maxn,MAXM,MAXN; 
global tord1, tord2; 
MAXM := fac([LeadingCoefficient(f, tord2)], tord1); 
maxm := max(seq(degree(MAXM[i]), i = 1 .. nops(MAXM))); 
MAXN := fac([LeadingCoefficient(g, tord2)], tord1); 
maxn := max(seq(degree(MAXN[i]), i = 1 .. nops(MAXN))); 
RETURN(evalb(maxm < maxn)); 
end proc:
#################################################################Improved Minimal Dickson Basis ###############################
MDB1 := proc (L, R, T) 
local K, p, n, N, H, ll, M, MM, m, maxm, maxn, mm, M1, MAXM, MAXN, i, Mm, flag1, flag, j; 
global tord1, tord2; 
tord1 := R; 
tord2 := T; 
N := NULL; 
H := NULL; 
K := L; 
while K <> [] do 
   ll := NULL; 
   p := K[1]; 
   for i to nops(K) do 
      if divide(LeadingMonomial(K[i], T), LeadingMonomial(p, T)) then 
         ll := ll, K[i];
      elif divide(LeadingMonomial(p, T), LeadingMonomial(K[i], T)) then 
         H := H, p; 
      end if; 
   end do; 
   K := [op(`minus`({op(K)}, {ll}))]; 
   N := N, p; 
end do; 
M := [op(`minus`(`union`({op(K)}, {N}), {H}))]; 
Mm := [op(`minus`({op(L)}, {op(M)}))]; 
MM := sort(Mm, comp); 
flag1 := false; 
for i to nops(M) do 
   flag := false; 
   while flag = false do 
      MAXM := fac([LeadingCoefficient(M[i], T)], R); 
      maxm := max(seq(degree(MAXM[i]), i = 1 .. nops(MAXM))); 
      for j to nops(MM) do 
         if flag = false then 
            if evalf(LeadingMonomial(M[i], T) = LeadingMonomial(MM[j], T)) then 
               MAXN := fac([LeadingCoefficient(MM[j], T)], R); 
               maxn := max(seq(degree(MAXN[i]), i = 1 .. nops(MAXN)));
               if maxn < maxm then 
                  M1 := subs(M[i] = MM[j], M); 
                  flag := true; 
                  flag1 := true; 
               end if; 
            end if; 
         end if; 
      end do; 
      flag := true; 
   end do; 
end do; 
if flag1 = false then 
   M1 := M; 
end if; 
RETURN(M1); 
end proc:
################################################################################################
with(LinearAlgebra):
with(Groebner): 
with(PolynomialIdeals):
############################################### Zero Dim-Check ################################################OK
Zdim := proc (E, f, R) 
local ns, rv, poly, d, M, EE, Inf, Rr, EG,ms,sms; 

LinearAlgebra; 
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); 
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; 

E := [op(EE)]; 
Rr:=tdeg(op(R));
LE:=LeadingMonomial(E,Rr);
V:=Groebner:-MaximalIndependentSet(LE,{op(R)}); 
d := nops(V); 
abar := [seq(i-1/3, i = 1 .. d)];

EV := seq(subs(V[i] = abar[i], E), i = 1 .. nops(V)); 
spE := Groebner:-Basis([EV][nops([EV])], Rr); 
   if IsZeroDimensional(`<,>`(op(spE))) then 
      if Zdim(spE, f, R) then 
         RETURN(true); 
      end if; 
   end if; 
RETURN(false); 
end proc:
########################################################## ICheck sub-algorithm ######################################################
Icheck := proc (E, f, R) 
local loop, p, i, s, M, m,Rr; 
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 1 Kapur tests(Zdim, Ccheck, Icheck) for a polynomial################
consist2 := proc (E, f, R) 
local EE, Inf, Rr, flag1, flag2; 

if nops(E) = 0 then 
   RETURN(true) 
end if; 
if IdealMembership(f, `<,>`(op(E))) then 
   RETURN(false) 
end if; 
EE := `<,>`(op(E)); 
Inf := IdealInfo[Variables](EE); 
Rr:=op(0,R)(op(Inf));

if IsZeroDimensional(EE) then 
   RETURN(Zdim(E, f, R)); 
end if; 
flag1 := Ccheck(E, f, R); 
if flag1 then  
   RETURN(true); 
end if; 
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;
f:=mul(i,i in N);
RETURN(consist2(EE,f,R)); 
end proc:
#############################################################
with(LinearAlgebra):
#########################################################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; 
if nops(LL) = 0 then 
   RETURN([true, LL, M]); 
end if; 
N:=LL;
if M = [] or M={} or M=[1] then 
   if Groebner:-Basis(LL,R)<>[1] then
      RETURN([true, LL, M]); 
   fi;
end if;
 
W := M; 
WW := seq(NormalForm(W[i], N, R), i = 1 .. nops(W)); 
WWW := seq(factor([WW][i]), i = 1 .. nops([WW])); 
test := true; 
t := true; 
NN := N; 
i := 1; 
flag := consist1(N, M, R); 
if flag = false then 
   RETURN([flag, N, [WWW]]) 
else 
   while t do 
      t := false; 
      facN := [seq([op(i)], `in`(i, [seq(fac([NN[i]]), i = 1 .. nops(NN))]))]; 
      facW := [op(fac([WWW], R))]; 
      NNN := []; 
      L := NULL; 
      for i to nops(facN) do 
         for x in facN[i] do 
            flag := true; 
            for w in facW while flag do 
               if divide(w, x) then 
                  flag := false; 
                  L := L, x; 
               end if; 
            end do; 
         end do; 
         NnN[i] := [op(`minus`({op(facN[i])}, {L}))]; 
         n[i] := simplify(expand(product(NnN[i][j], j = 1 .. nops(NnN[i])))); 
      end do; 
      NNNN := [seq(n[i], i = 1 .. nops(facN))]; 
     
      if flag=false then
         t := true; 
         NN := Groebner:-Basis(NNNN, R); 
         WW := seq(NormalForm([WWW][i], NN, R), i = 1 .. nops([WWW])); 
         WWW := seq(factor([WW][i]), i = 1 .. nops([WW])); 
      end if; 
   end do; 
end if; 
ww := NULL; 
for i to nops([WWW]) do 
   if type(WWW[i], 'constant') = true then  
   else 
      ww := ww, [WWW][i]; 
   end if; 
end do; 
WWW := ww; 
MM := [op(fac([WWW], R))]; 
WWW := [op({seq(i/LeadingTerm(i, R)[1], `in`(i, MM))})]; 
www := [seq(nrm([ee]), `in`(ee, WWW))]; 
if www <> [] then 
   www := [op(www)]; 
   www := [seq(op(jj), `in`(jj, www))]; 
end if; 
RETURN([flag, NN, www]); 
end proc:
######################################
#####################################################input: parametric polynomial f, output: variables appeared in f###############################
Variables := proc (e, p::(({list, rtable, set})(name))) 
`minus`(indets(e, And(name, Not(constant))), convert(p, set)):
end proc:
#####################################################input: parametric poly ideal F, output: parameters appeared in F ###############################
variables:=proc(order)
	local l;
	l:=[op(order)];
	if type(l[1],`list`)=true then
		RETURN(l[2]);
	else
		RETURN(l);
	end if:
end:

pars := proc (F, T) 
local p, Ps; 
p:=variables(T);

Ps := `union`(seq(Variables(F[i], p), i = 1 .. nops(F))); 
RETURN(Ps); 
end proc:
##################################################################Return true if a is constant#########################################
isconstant := proc (a) RETURN(type(a, 'constant')); end proc:
##############################################Return true if polynomial f is w.r.t variables 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 polynomial f is w.r.t parameters 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 for product of h_i's###################################################
cunion := proc (H, E, N) 
local EE, NN, k, List, i, h; 
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:
##################################################the correction of consistent checking in Kapur(when check N*Gr)#######################
mycheck := proc (E, N, Gr) 
local i;
global ordp;
for i to nops(Gr) do 
   if ccc(E, zarb(N, [Gr[i]]), ordp)[1] = true then 
      RETURN([true,E,zarb(N,[Gr[i]])]) ;
   end if; 
end do; 
RETURN([false,1]); 
end proc:
#############################################################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;

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));
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(type, H, 'constant');
   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)#########################
PGBMAIN := proc(E,N,F,R,T)
local AA;
global PGB,LIST;
PGB:=NULL;
AA:=pgbmainc(E,N,F,R,T);	
end:

#*******************************************************************************
#*******************************************************************************
#*****************Parametric Generic Groebner Walk Algorithm********************
#*******************************************************************************
#*******************************************************************************

with(Groebner):
with(LinearAlgebra):
with(RegularChains):

####################### This function returns list of variables of a monomial order ###############################

variables:=proc(order)
	local l;
	l:=[op(order)];
	if type(l[1],`list`)=true then
		RETURN(l[2]);
	else
		RETURN(l);
	end if:
end:

####################### This function returns list of parametric coefficients of a polynomial #######################
#######################  'f' is parametric polynomilal and 'var' is  'var' is variables appeared in  f ###############
Pcoeff:=proc(f,var) 
	local c,p,N,i;
	c:=[coeffs(simplify(f),var,'e')];
	p:=[e];
	N:={seq([c[i],p[i]],i=1..nops(c))};
	RETURN(N);
end:

####################### This function eliminates zero terms of a parametric polynomial ########################################
####################### 'f' is parametric polynomial and 'par' is list of parameters f and 'N' is null condition set for f######

ValidTerms:=proc(f,par,N) 
      
        RETURN(Groebner[NormalForm](f,N,plex(op(par))));
end:

nonZeroTerms:=proc(G,var,par,N)
	local r,i,k;
	r:=NULL:
	for i from 1 to nops(G) do
		k:=ValidTerms(G[i],par,N);
		if k <> 0 then
			r:=r,collect(k,var,distributed);
		fi:
	od:
	RETURN([r]);
end:

####################### This function returns vector of degrees of terms of a polynomial ###################################
####################### Input:'f' is parametric polynomilal and 'var' is  'var' is variables appeared in  f ################
support := proc (f, var) 
        local L,N,i,a,n,j,v,pc;
        pc:=Pcoeff(f,var);
        L:=[seq(pc[i][2],i=1..nops(pc))]; 
        n:=[seq([seq(degree(L[i],var[j]),j=1..nops(var))],i=1..nops(L))];
	N:={seq([Vector(n[i]),pc[i][1]*pc[i][2]],i=1..nops(n))};
	RETURN(N); 
end proc:

########################  Implementation of Facet preorder defined in page 8 of paper #######################################
######################## Input:'u' and 'v' are two vectors, 'order1' and 'order2' are two term orders #######################
FacetPreOrder :=proc(u,v,order1,order2) # 
		
	local var,n,T,tv,tu,Tuv,Tvu,t,a,term,i;
	var:=variables(order2);
	n:=nops(var);
	T:=MatrixOrder(order2);
	T:=convert(T,Matrix);
	tv:=Transpose(v);
	tu:=Transpose(u);                                           
	Tuv:=T.u.tv;
	Tvu:=T.v.tu;
	t:=Tvu-Tuv;
	
	for i from 1 to n do
		a:=Row(t,i);
		term:=mul(var[j]^a[j],j=1..n);
		
		if term <> 1 then
			
			if TestOrder(1,term,order1)=false then 
				RETURN(false);
			else
				RETURN(true);
			fi:
			
		fi:
	od:
	RETURN(true);
end:

######################## Compute leading part ideal of G w.r.t. vector w defined in page 8 of paper #########################
######################## Input: 'g' is a parametric ideal, 'w' is a vector, 'order1' and 'order2' are two term orders #######
initial_g := proc (g,w,order1,order2) 
	
	local var,supp,L,ltsupp,lt,N,i,terms,init,k,ltSupp;
	var:=variables(order2);
	supp:=support(g[2],var);
	L:=[seq(supp[i][1],i=1..nops(supp))];
	ltSupp:=support(g[1],var);
	lt:=ltSupp[1][1];
	N:=NULL:
	
        for i from 1 to nops(L) do
        	if Equal(L[i],lt)=false then
        		k:=lt-L[i];
        		if FacetPreOrder(k,w,order1,order2)=true and FacetPreOrder(w,k,order1,order2)=true then
        			N:=N,supp[i][2];
        		fi:
        	fi:
         od:
         N:={N};
         init:=add(N[i],i=1..nops(N));
         RETURN([g[1],collect(g[1]+init,var)]);
end:

initial_G:=proc(G,w,order1,order2)
	RETURN({seq(initial_g(G[i],w,order1,order2),i=1..nops(G))});
end:

############################## Compute Gamma(F) defined in page 8 of paper ##########################################
############################# Input : 'f' is parametric polynomial, 'var' is list of variables of polynomial f ######
Rond_f:=proc(f,var) 
	
	local L,supp,N,s,lt,i;
	L:=support(f[2],var);
	supp:={seq(L[i][1],i=1..nops(L))};
	N:=NULL:
	s:=support(f[1],var);
	lt:=s[1][1];
	for i from 1 to nops(supp) do
		if Equal(lt,supp[i])=false then
			N:=N,lt-supp[i];
		fi:
	od:
	RETURN(N);
end:

######################## Input :'F' is list of polynomilas, 'var' is variables appeared in  F #######################
Rond_F:=proc(F,var) 
	RETURN({seq(Rond_f(F[i],var),i=1..nops(F))});                           
end:

############################# Compute marked Groebner basis ##########################################################
############################# Input : 'G' is a Groebner basis, 'order' is a term order ###############################

mark:=proc(G,order) 
	RETURN([seq([LeadingCoefficient(G[i],order)*LeadingMonomial(G[i],order),G[i]],i=1..nops(G))]);
end:

########### If vector v is in Delta(G) defined in page 8 of paper, this function returns true, else returns false ####
########## Input : 'v' is a vector, 'order1' and 'order2' are two term orders ########################################

IsInIntersect:=proc(v,order1 , order2) 
	
	local var,term;
	var:=variables(order2);
	term:=mul(var[j]^v[j],j=1..nops(var));
	if TestOrder(1,term,order1)=true and TestOrder(term,1,order2)=true then
		RETURN(true);
	else
		RETURN(false);
	fi:
end:

############################## Find minimal gamma in Gamma(G) in PGGW algorithm of section 4 of paper ###############
############################## Input : 'G' is a Groebner basis, 'order1' and 'order2' are two term orders ###########
Compute_last_w:=proc(G,order1,order2)  
	
	local var,rond,m,n,i;
	var:=variables(order2);
	rond:=Rond_F(G,var);
	m:=-1;
	n:=NULL:
	
	for i from 1 to nops(rond) do
		if IsInIntersect(rond[i],order1,order2)=true then
			n:=n,rond[i];
		fi:
	od;
	n:=[n];
	if n <> [] then
		m:=n[1];
	fi:
	
	for i from 2 to nops(n) do
		if FacetPreOrder(n[i],m,order1,order2)=true then
			m:=n[i];
		fi:
	od:
	RETURN(m);
end:
	
############################## PNormalForm Algorithm defined in page 8 of paper ######################################


lt:=proc(f,T)

	RETURN(LeadingTerm(f,T)[1]*LeadingTerm(f,T)[2]):
end:

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));
	od:
	RETURN(L);
end:

poly:=proc(f,var)

	local g,fl;
	g:=collect(f,var);
	fl:=FL(g,plex(op(var)));
	RETURN(fl);
end:


Division:=proc(f,inl,L::list,var)
	
	local B,p,flag,i,r,ff,pl,ord,lm,lc;
	B:=[seq(0,i=1..nops(L))];
	ff:=collect(f,var);
	p:=poly(ff,var);
	r:=0;
	ord:=plex(op(var));
	while p <> [] do
		pl:=p[-1];
		i:=1;
		flag:=false;
		while i<= nops(L) and flag=false do
			lm:=LeadingMonomial(inl[i],ord);
			lc:=LeadingCoefficient(inl[i],ord);
			if divide(pl,lm,'q') then 
				B[i]:=collect(simplify(expand(B[i]+q/lc)),var);
				ff:=collect(simplify(expand(ff-simplify(expand((q/lc*L[i]))))),var);
				p:=poly(ff,var);
				flag:=true;
			else
				i:=i+1;
			fi:
		od:
		if flag=false then
			r:=r+pl;
			ff:=collect(simplify(expand(ff-pl)),var);
			p:=subsop(-1=NULL,p);
		fi:
	od:
	RETURN(B,r);
end:
################################# sub-algorithms used in the H_Prime algorithm #################################
coeffOfLT:=proc(lt,f,var)

	local c1,c2,i,t;
	c1:=Pcoeff(lt,var)[1];
	c2:=Pcoeff(f,var);
	for i from 1 to nops(c2) do
		t:=c2[i][2];
		if t=c1[2] then
			RETURN(c2[i][1]*c2[i][2]);
		fi:
	od:
end:

simCoeff:=proc(f,var)
	local c,s,lc,sim_f;
	c:=[coeffs(f,var)];
	s:=seq(expand(denom(c[i])),i=1..nops(c));
	lc:=lcm(s);
	sim_f:=collect(f*lc,var,distributed);
	RETURN(simplify(sim_f));
end:
	
################################# H'={f-f^G | f in H} where H=GB(in_w(I)) over <_2 and G is Grobner basis of last cone ##########
H_Prime:=proc(G,H,var)

	local N,i,k,sim,lt,inl,L,s;
	N:=NULL:
	inl:=[seq(G[i][1],i=1..nops(G))];
	L:=[seq(G[i][2],i=1..nops(G))];
	for i from 1 to nops(H) do
		k:=[Division(H[i][2],inl,L,var)];
		sim:=collect(H[i][2]-k[2],var,distributed);
		s:=collect(simCoeff(sim,var),var,distributed);
		if sim <> 0 then
			lt:=coeffOfLT(H[i][1],s,var);
			N:=N,[lt,s];
		fi:
	od:
	RETURN({N});
end:

##################################  Interreduce function in PGGW algorithm of section 4 of paper #################################
################################## Input :  'G' is a Groebner basis, 'var' is variables appeared in G  ###########################

interReduced:=proc(G,var) # 
	
	local i,a,r,L,g,lt,inl,k,result,l;
	inl:=[seq(G[i][1],i=1..nops(G))];
	g:=[seq(G[i][2],i=1..nops(G))];
	if nops(G)=1 then
		RETURN([[inl[1],g[1]]]);
	fi:
	for i from 1 to nops(g) do
		a:=g[i];
		L:=subsop(i=NULL,g);
		l:=subsop(i=NULL,inl);
		k:=[Division(a,l,L,var)];
		
		r:=collect(simCoeff(k[2],var),var,distributed);
		if r <> 0 then
			lt:=coeffOfLT(inl[i],r,var);
			inl:=subsop(i=lt,inl);
			g:=subsop(i=r,g);
		fi:
	od:
	
	result:=[seq([inl[i],g[i]],i=1..nops(g))];
	RETURN(result);
end:	

################################## Parametric Generic Groebner Walk (PGGW) algorithm of section 4 of paper ######################################
################################## Input : 'G' is first Groebner basis, 'N' is null condition set, 'W' is non-null condition set,################
##################################  'order1' is first term order, 'order2' is second term order, 'POrder' is parametric order ###################
PGGW:=proc(G,N,W,order1,order2,POrder) 
				       																
	
	local gb,reduced,var,w,g,initial_w,initial,H,HP,ma,sim,par,m,i,j,r,Hmark,sim_reduce,result,pgb,d,flag,h,si,rr,Result,Rh,te,temp,lts;
	var:=variables(order2);
	w:=1;
	flag:=true;
	par:=variables(POrder);	
	sim:=nonZeroTerms(G,var,par,N);
	ma:=mark(sim,order1);
	rr:=interReduced(ma,var);
	g:=[[N,W,rr]];
	Result:=NULL:
	while g <> [] do
		gb:=g[1];
		g:=subsop(1=NULL,g);
		w:=Compute_last_w(gb[3],order1,order2);
		
		
		if w <> -1 then
			initial_w:=initial_G(gb[3],w,order1,order2);
			initial:=[seq(initial_w[z][2],z=1..nops(initial_w))];
			H:={PGBMAIN(gb[1],gb[2],initial,POrder,order2)};
			
			for j from 1 to nops(H) do
				h:=H[j];
				si:=nonZeroTerms(h[3],var,par,h[1]);
				Rh:=Groebner[Basis](si,order2);
				Hmark:=mark(Rh,order2);
				HP:=H_Prime(gb[3],Hmark,var);
				reduced:=interReduced(HP,var);
				r:=[seq(reduced[z][2],z=1..nops(reduced))];
				sim_reduce:=nonZeroTerms(r,var,par,h[1]);
				lts:=NULL;
				lts:=[seq(coeffOfLT(reduced[kk][1],sim_reduce[kk],var),kk=1..nops(reduced))];
				result:=[seq([lts[z],sim_reduce[z]],z=1..nops(reduced))];
				g:=[op(g),[h[1],h[2],result]];
			
			od:
		else
			te:=gb[3];
			temp:=[seq(te[i][2],i=1..nops(te))];	
			Result:=Result,[gb[1],gb[2],temp];
		fi:
		  
	od:
	RETURN({Result});
end:

#################################### sub-algorithm used in the PGSC algorithm ###########################################
Pwalks:=proc(result,order1,order2,porder)
        
	local N;
	N:={seq(op(PGGW(result[i][3],result[i][1],result[i][2],order1,order2,porder)),i=1..nops(result))} minus {[{1}, [], [0]]};
	RETURN(N);
end:
#################################### Parametric Groebner System Conversion (PGSC) algorithm in page 11 of paper ##########

PGSC := proc ( ideal, N, W, porder, order1, order2 )
	
 	local G,P;
	G := [PGBMAIN(N, W, ideal, porder, order1)];
	P := Pwalks(G, order1, order2, porder);
	RETURN(P);
end:
PGSC( [a*x-b*y, y^2+c*x*y-1] , [], [1], plex(a, b, c, d, m, n , r, t), tdeg(x, y, z, u), plex(x, y, z, u)):	

#*******************************************************************************
#*******************************************************************************
#****************** Check correctness of PGSC Algorithm ************************
#*******************************************************************************
#*******************************************************************************

CHECK := proc (G,N,W, GS1, R, T, m, n, number) 
	
	local P, flag, rd, f, RD, i, G1, G2, GG, GS, GS2, FLAG,test,can;
	FLAG := NULL; 
	P := [op(R)]; 
	GS2 := NULL; 
	test:=NULL:
	for i to nops(GS1) do
		can:=canspec(GS1[i][1], GS1[i][2], plex(op(R)));
		GS2 := GS2, [can[2], can[3], [seq(expand(g/LeadingCoefficient(g, T)), `in`(g, GS1[i][3]))]]; 
	end do; 
	
	GS := [GS2]; 
	f := rand(m .. n); 
	RD := [seq([seq(f(), j = 1 .. nops(P))], i = 1 .. number)]; 
	for rd in RD do 
		flag := false; 
		if  Groebner:-Basis(subs(seq(P[j] = rd[j], j = 1 .. nops(P)), G), T) <> [1] and subs(seq(P[j] = rd[j], j = 1 .. nops(P)), [op(N)]) = [seq(0, i = 1 .. nops(N))] and evalb(`in`(0, subs(seq(P[j] = rd[j], j = 1 .. nops(P)), W)))=false then
		
			
		for i to nops(GS) while flag = false do 
			if subs(seq(P[j] = rd[j], j = 1 .. nops(P)), [op(GS[i][1])]) <> [seq(0, i = 1 .. nops(GS[i][1]))] or evalb(`in`(0, subs(seq(P[j] = rd[j], j = 1 .. nops(P)), GS[i][2]))) then  
			else 
				flag := true; 
				GG := Groebner:-Basis(subs(seq(P[j] = rd[j], j = 1 .. nops(P)), G), T); 
				G1 := [seq(expand(g/LeadingCoefficient(g, T)), `in`(g, GG))]; 
				G2 := subs(seq(P[j] = rd[j], j = 1 .. nops(P)), GS[i][3]); 
				if {op(G1)} = {op(G2)} then 
					FLAG := FLAG, true; 
				else 
					FLAG := FLAG, false; 
					test:=test,[rd,GG,G1,G2];
				end if; 
			end if; 
		end do;
		else
			flag:=true;
		fi:
	if flag = false then 
		RETURN(ERROR, rd); 
	end if; 
	end do; 
	RETURN([FLAG],[test]); 
end proc:

#*******************************************************************************
#*******************************************************************************
#******************************* Examples **************************************
#*******************************************************************************
#*******************************************************************************

print("*************************EXAMPLES Page 13 of paper (Section 5)**************************");

EX1 := [b*z^3+c*y^2, m*y^3-r*x*y+z^2, n*x^3+a*y^2];
EX2 := [c*y+d+x, b*y^2+c*x^2+2*c*x*y+2*d*x+2*m*y+n, b*y+m*x-1];
EX3 := [a*y^2+a*x+c*z, 2*a*y-x^2-x, -a*b*y^2-a*z^3+3*b*x^2];
EX4 := [b*z^2+c*y, a*z+c*x-x*y+b+y, b*x^2+a*y];
EX5 := [(m*n-a)*x+y^2+a, c*x+z^2, x^2+(a*d-c)*y*x+a*b, u^2+(a^4-4)*x];
EX6 := [m*x^2+b, c*y+d+x-1, b*y^2+2*c*x*y+2*d*x+x+2];
EX7 := [r*x^5+(a*b-c)*z-n, c*y^3+a*c*x+d*n, z^3-(c-t)*y];
EX8 := [a*b*x*y+a*y^3-c*z+1, a*x^2+a*x+c*u, t*u^3+t*u, b*z^3+m*n*x-b*z];
EX9 := [(a*b-c)*z+r*x^4-n, y^2+(a*c-b)*x+d, z^2-(b*c-a)*y+t];
EX10 := [u^3+(a-1)*x*b, a*x*y+x^3+m*y-1, c*y^2+n*x+a, a*z^2+c*x];


#PGSC(EX1, [], [1], plex(a, b, c, m, n, r), tdeg(x, y, z), plex(x, y, z));
#PGSC(EX2, [], [1], plex(b, c, d, m, n), tdeg(x, y), plex(x, y));
#PGSC(EX3, [], [1], plex(a, b, c), tdeg(x, y, z), plex(x, y, z));
#PGSC(EX4, [], [1], plex(a, b, c), tdeg(x, y, z), plex(x, y, z));
#PGSC(EX5, [], [1], plex(a, b, c, d, m, n), tdeg(x, y, z, u), plex(x, y, z, u));
#PGSC(EX6, [], [1], plex(b, c, d, m), tdeg(x, y), plex(x, y));
#PGSC(EX7, [], [1], plex(a, b, c, d, n, r, t), tdeg(x, y, z), plex(x, y, z));
#PGSC(EX8, [], [1], plex(a, b, c, m, n, t), tdeg(x, y, z, u), plex(x, y, z, u));
#PGSC(EX9, [], [1], plex(a, b, c, d, n, r, t), tdeg(x, y, z), plex(x, y, z));
#PGSC(EX10, [], [1], plex(a, b, c, m, n), tdeg(x, y, z, u), plex(x, y, z, u));