#*******************************************************************************
#*******************************************************************************
#***************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 
with(Ore_algebra):
func:=proc(order1,order2,l,par) 
  global term1:=order1;  
  global term2:=order2;  
  global var:=l;   
  global li:=[op(par),op(l)];
  aa := poly_algebra(op(li)):
  pp:=proc(t1,t2)
    global li,term1,term2,var; 
    co1:=coeffs(t1,var,'e');    
    te1:=e;      
    co2:=coeffs(t2,var,'e');   
    te2:=e;
    `or`(`and`(TestOrder(te1, te2, term1) = true, te1 <> te2), `and`(te1 = te2, TestOrder(co1, co2, term2) = true))
    end:
    T := MonomialOrder(aa, ('user')(pp, li));
     RETURN(T);  
end:



############################
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;
var:=variables(T);
par:=variables(R);
TT:=func(T,R,var,par);
ordp:=R;
ordv:=T;
E:=Groebner:-Basis(EE,R);#print(E, N, R);
CC := ccc(E, N, R); 
if CC[1] = false then 
   RETURN([]); 
else 
   #G:=Groebner:-Basis([op(F),op(CC[2])],prod(T,R));
    G:=Groebner:-Basis([op(F),op(CC[2])],TT);
end if; 
if pars(F,T)={} then 
   RETURN([[E,N,F]] );
fi;
Gr:=select(issubsetp,G);
 
GG:=remove(issubsetp,G,'constant');
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,'constant'); 
   print("hi");
   h:=lcm(op(H));
   Zh := zarb(N, [h]); 
   CCCC := ccc(Gr, [op(Zh)], R); 
   if CCCC[1]=true then 
      PGB := PGB, [CCCC[2], factor(CCCC[3]), Gm];
      CU:=[cunion(H,CCC[2],N)];print(injaaaaaaaaaaaaaaaaaaaaaaa);
   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:
#####################
pgbmainc2 := 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;
var:=variables(T);
par:=variables(R);
TT:=func(T,R,var,par);
ordp:=R;
ordv:=T;
E:=Groebner:-Basis(EE,R);#print(E, N, R);
CC := ccc(E, N, R); 
if CC[1] = false then 
   RETURN([]); 
else 
   G:=Groebner:-Basis([op(F),op(CC[2])],prod(T,R));
   # G:=Groebner:-Basis([op(F),op(CC[2])],TT);
end if; 
if pars(F,T)={} then 
   RETURN([[E,N,F]] );
fi;
Gr:=select(issubsetp,G);
 
GG:=remove(issubsetp,G,'constant');
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,'constant'); 
   print("hi");
   h:=lcm(op(H));
   Zh := zarb(N, [h]); 
   CCCC := ccc(Gr, [op(Zh)], R); 
   if CCCC[1]=true then 
      PGB := PGB, [CCCC[2], factor(CCCC[3]), Gm];
      CU:=[cunion(H,CCC[2],N)];print(injaaaaaaaaaaaaaaaaaaaaaaa);
   else
      CU:=[cunion(H,CCC[2],CCC[3])];
   end if; 
   LIST:=[seq(pgbmainc2(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)#########################
PGBMAIN2 := proc(E,N,F,R,T)
local AA;
global PGB,LIST;
PGB:=NULL;
AA:=pgbmainc2(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) 
        print(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])],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 #######


############################## Compute Gamma(F) defined in page 8 of paper ##########################################
############################# Input : 'f' is parametric polynomial, 'var' is list of variables of polynomial f ######


######################## 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;
  #option trace;
	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;
  #option trace;
  print (G,H,var);
	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:
	

#*******************************************************************************
#*******************************************************************************
#****************** 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:

##############################
checkPoint:=proc(GS,rd,par)
    #option trace;
    flag:=false;
    k:=0;
    print(nops(GS));
    for i to nops(GS) while flag = false do 
      print("hi");
      if subs(seq(par[j] = rd[j], j = 1 .. nops(par)), [op(GS[i][1])]) <> [seq(0, j = 1 .. nops(GS[i][1]))] or evalb(`in`(0, subs(seq(par[j] = rd[j], j = 1 .. nops(par)), GS[i][2]))) then  
        print("if");
      else  
        prin("else");  
        flag:=true;
        k:=i;
      fi:
    
    od:
    if flag=true then
      RETURN(k);
    else
      RETURN(false);
    fi:
end:
checkConditions:=proc(N,W,rd,par)
  if subs(seq(par[j] = rd[j], j = 1 .. nops(par)), [op(N)]) <> [seq(0, j = 1 .. nops(N))] or evalb(`in`(0, subs(seq(par[j] = rd[j], j = 1 .. nops(par)), W))) then 
    RETURN(false);
  else
    RETURN(true);
  fi:
end:

getIdeal:=proc(ideal,rd,par)
    RETURN(subs(seq(par[j] = rd[j], j = 1 .. nops(par)), ideal));
end:

getPoint:=proc(m,n,number,par)
    f := rand(m .. n); 
    RD := [seq([seq(f(), j = 1 .. nops(par))], i = 1 .. number)];
    RETURN(RD);
end:

CheckIdeal:=proc(ideal,rd,par,order)
option trace;
  ii:=getIdeal(ideal,rd,par);
  gb:=Groebner[Basis](ii,order);
  if gb=[1] then
    RETURN(false);
  else
    RETURN(true);
  fi:

end:
conditions:=proc(N,W)
  local L:
  L:=NULL:
  L:=seq(N[i]=0,i=1..nops(N));
  L:=L,seq(W[i] <> 0 , i=1..nops(W));
  RETURN([L]);
end:

samplepoint:=proc(inequality,var)
  option trace;
  local R,n,c,i,sp; 
  R:=PolynomialRing(var);
  sp:=SamplePoints(inequality,R,output=record);
  if  sp <> NULL then
    c:=convert([sp][1], list);
    n:=NULL:
    for i from 1 to nops(c) do

      if type(op(2,op(1,c[i])),'list') then
        cc:=op(2,op(1,c[i]))[1];
      else
        cc:=op(2,op(1,c[i]));
      fi:
      n:=n,cc;
    od:
    RETURN([n]);
  else
    RETURN(-1);
  fi:
end:

#*******************************************************************************
#*******************************************************************************
#******************************* Examples **************************************
#*******************************************************************************
#*******************************************************************************
IsInSystem := proc (G, w, l) 
	s := [seq(l[i] = w[i], i = 1 .. nops(l))];
  count := 0;
	for i to nops(G) do 
		u := eval(G[i], s);
    uu:=op(2,G[i]);
    if eval(uu,s)=0 then
      count:=count+1;
    fi:
		if  evalb(u) = false then 
			RETURN(false);
		end if:
	end do;
  if count<2 then 
	   RETURN(true);
  else
    RETURN(false);
  fi:
end proc:
#u <> (0 < 0) and
############################
NumberOfValidPoint:=proc(GS, w, l,facet)
	count:=0;
	for i from 1 to nops(GS) do
    a := GS[i];
		if IsInSystem(a[2],w,l) = true then
      flag := true;
      
          if a[1]=1 then
            flag :=false;
          else
            for j from 1 to nops(l) do
              if a[1][j] <> facet[j] then
                flag := false;
              fi:
            od:
          fi:
      if flag = false then
			 count:=count+1;
      fi:
			
	  fi:	
  od:
	RETURN(count);
	
end:
		
###############################
isValid:=proc(v1,v2,var,order1)
  varr:=[op(order1)];
	n:=nops(var);
	term1:=mul(varr[j]^v1[j],j=1..n);
	term2:=mul(varr[j]^v2[j],j=1..n);
	if TestOrder(term2,term1,order1)=false then 
		RETURN(false);
	else
		RETURN(true);
	fi:
end:

validFacet:=proc(vect,ine,v,facet,n,order)
	if isValid(vect[1],vect[2],n,order) and NumberOfValidPoint(ine,v,n,facet) < 2 then
		RETURN([isValid(vect[1],vect[2],n,order),NumberOfValidPoint(ine,v,n,facet) < 2,0]);
	else
		RETURN([isValid(vect[1],vect[2],n,order),NumberOfValidPoint(ine,v,n,facet) < 2,1]);
	fi:
end:
##############################
GB:=proc(i,order)
	RETURN(Groebner[Basis](i,order));
end:

reducing:=proc(i,order)
  RETURN(Groebner[InterReduce](i,order));
end:

elim:=proc(l,i)
	RETURN(subsop(i=NULL,l));
end:

init := proc (G, v,l) 
  option trace;
	N := NULL; 
	for i to nops(G) 
		do N := N, InitialForm(G[i][2], v, l) 
	end do; 
	RETURN([N]); 
end proc:

singularOrder:=proc(l)
	RETURN(wp(op(l)));
end:

vec:=proc(v)
	RETURN('v');
end:

ListToSet:=proc(l)
	RETURN([op({op(l)})]);
end:

mapleList:=proc(l)
	temp:=NULL:
	for i from 1 to nops(l) do
		temp:=temp,op(l[i]);
	od:
	RETURN([temp]);
end:

#########################


flip:=proc(G,v,al,var,par)
    option trace;
    local ini,alp,gg,final,si,gb,Hmark,HP,reduced,r,consist,sim_reduce,lts,result;
    ini := init(G[3],v,var);
    alp:= -1*al;
    gg := PGBMAIN(G[1], G[2], ini, plex(op(par)), wdeg(alp, var));
    nonZeroGB:= PGBMAIN([], par, ini, plex(op(par)), wdeg(alp, var)); 
    final := NULL;
    for j in gg do
        
        si := nonZeroTerms(j[3],var,par,j[1]);
        #gb :=Groebner[InterReduce](si,wdeg(alp,var));
        gb :=Groebner[Basis](si,wdeg(alp,var));

        Hmark := mark(gb,wdeg(alp,var));
        HP := H_Prime(G[3],Hmark,var);
        #return(HP);
        reduced := interReduced(HP,var);
        r:=[seq(reduced[z][2],z=1..nops(reduced))];
        consist:=canspec(j[1],j[2],plex(op(par)));
        if consist[1]=true then
          if consist[3]=[] then
            c3:=[1];
          else
            c3:=consist[3];
          fi:
          sim_reduce:=nonZeroTerms(r,var,par,consist[2]);
          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))];
          final := final , [consist[2],c3,result];
        fi:
    od:
    RETURN([final]);
end:

Gflip:=proc(G,v,al,var,par)
    option trace;
    local ini,alp,gg,final,si,gb,Hmark,HP,reduced,r,consist,sim_reduce,lts,result;
    ini := init(G[3],v,var);
    alp:= -1*al;
    
    nonZeroGB:= PGBMAIN([], par, ini, plex(op(par)), wdeg(alp, var)); 
    final := NULL;
    for j in nonZeroGB do
        
        si := nonZeroTerms(j[3],var,par,j[1]);
        gb :=Groebner[Basis](si,wdeg(alp,var));

        Hmark := mark(gb,wdeg(alp,var));
        HP := H_Prime(G[3],Hmark,var);
        #return(HP);
        reduced := interReduced(HP,var);
        r:=[seq(reduced[z][2],z=1..nops(reduced))];
        #consist:=canspec(j[1],j[2],plex(op(par)));
        if consist[1]=true then
          if consist[3]=[] then
            c3:=[1];
          else
            c3:=consist[3];
          fi:
          sim_reduce:=nonZeroTerms(r,var,par,consist[2]);
          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))];
          final := final , result;
        fi:
    od:
    RETURN(final);
end:

################################# cartprodrt product
with(combinat, cartprod):

cart := proc( list2)
  local T,N,RESULT;
  T := cartprod(list2);
  N := NULL:
  while not T[finished] do N:=N,T[nextvalue]() end do:
  RESULT := NULL:
  N:=[N];
  #for i from 1 to nops(N) do
  #    RESULT := RESULT , [op(list1),op(N[i])];
  #od:    
  RETURN(N):
end:

vec:=proc(v1,v2)
  if evalb(v1=v2)=true then
    RETURN(true);
  else
    RETURN(false);
  fi:
end:
###################################
SimplifyOutput:=proc(PGF,par)
  local finalN,finalW;
  option trace;
  Final:=NULL:
  for i from 1 to nops(PGF) do
      ele:=PGF[i];
      finalN:=NULL:
      finalW:=NULL:
      output:=NULL:
      for j from 1 to nops(ele) do
          a:=ele[j];
          finalN:=finalN,op(a[1]);
          finalW:=finalW,op(a[2]); 
          finalG:=NULL:
          for k from 1 to nops(a[3]) do
            finalG:=finalG,a[3][k][2];

          od: 
          output:=output,[finalG];
      od:
      can:=ccc([finalN],[finalW],plex(op(par)));
      if can[1]=true then
        if can[3]=[] then
          c3:=[1];
        else
          c3:=can[3];
        fi:
        Final:=Final,[can[2],c3,[output]] ;
      fi:
  od:
  RETURN([Final]);
end:
####################################
ch :=proc(bb,par)
  N:=NULL:
  W:=NULL:
  for i from 1 to nops(bb) do
    ele := bb[i][2];
    N:=N,op(ele[1]);
    W:=W,op(ele[2]);
  od:
  can:=ccc([N],[W],plex(op(par)));
  if can[1]=true then
    return(true);
  else
    return(false);
  fi:
end:
######################################
simCoeff1:=proc(f,var)
  #option trace;
  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);
  gcdd:=igcd(coeffs(sim_f,var));

  RETURN(simplify(sim_f/gcdd));
end:
#######################################
sim:=proc(G,var)
  N:=NULL:
  for i from 1 to nops(G) do
    #a:=simCoeff1(G[i],var);
    a:=Groebner[InterReduce]([G[i]],plex(op(var)));
    N:=N,a;
  od:
  RETURN({N});
end:
#######################################
SetOfLT:=proc(set,var)
  N:=NULL:
  for i from 1 to nops(set) do
      #p:=coeffs(set[i][1],var,'e');
      N:=N,set[i][2];
  od:
  RETURN({N});
end:
#######################################
equalSet:=proc(set1,set2)
  if evalb(set1=set2) then
      RETURN(true);
  else
      RETURN(false);
  fi:
end:
########################################
UniounCondition := proc(N1,W1,N2,W2,par)
  NewN:=[op(N1),op(N2)];
  NewW:=[op(W1),op(W2)];
  can := canspec(NewN,NewW,plex(op(par)));
  if can[1]=true then
    return([can[2],can[3]])
  else
    return(false)
  fi:
end:
#########################################
EqualityOfTwoList := proc(list1,list2)
  if nops(list1)<> nops(list2) then
    return false;
  else
    for i from 1 to nops(list1) do
        if member(list1[i],list2)=false then
          return false;
        fi:
    od:
    return true;
  fi:
  end:
##########################################
  EqualityOfTwoList2 := proc(list1,list2)
  option trace;
  if nops(list1[3])<> nops(list2[3]) then
    return false;
  else
    for i from 1 to nops(list1[3  ]) do
        if member(list1[3][i],list2[3])=false then
          return false;
        fi:
    od:
    return true;
  fi:
  end:
############################################  
MemberList:=proc(l,List)
option trace;
for i from 1 to nops(List) do
  a:=List[i][2];
  if EqualityOfTwoList2(a,l)=true then
    return true;
  fi:
od:
return false;
end:
############################################
  Permute := proc (n) 
     result := NULL;
     for i to n do
       l := seq(1, k = 1 .. i);
       for j to n-i do 
          l := l, 0 ;
       end do: 
       p := combinat:-permute([l], n);
       result := result, op(p) 
     end do:
    for i to n do
      l := seq(0, k = 1 .. i);
      for j to n-i do
         l := l, 1; 
      end do:
      p := combinat:-permute([l], n); 
      result := result, op(p);
    end do: 
    return(result);
 end proc:
#############################################
getPoint2:=proc(m,n,number,par)
    option trace;
    f := rand(m .. n); 
    p := Permute(nops(par));
    RD := [p,seq([seq(f(), j = 1 .. nops(par))], i = 1 .. number)];
    RETURN(RD);
end: