#########################################
########			  ####### 
########  F5 + Montes algorithms  #######
########                          #######
#########################################

with(Groebner):
with(Algebraic):
with(PolynomialIdeals):
############################Select Algorithm###################### 
sff := proc (x) 
evalb(type(x, 'constant')): 
end proc:
################
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; 
#option trace; 
if F = [] then 
RETURN(F); 
end if; 
KK := F; 
if type(denom(KK[1]), monomial) then 
h := LeadingMonomial(denom(KK[1]), plex(op(indets(denom(KK[1]))))); 
else 
h := denom(KK[1]); 
end if; 
for i from 2 to nops(KK) do 
if type(denom(KK[i]), monomial) then 
h := lcm(h, LeadingMonomial(denom(KK[i]), plex(op(indets(denom(KK[i])))))); 
else 
h := lcm(h, denom(KK[i])); 
end if; 
end do; 
NM := simplify(expand(h*KK)); 
RETURN(NM); 
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; 
#option trace;
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:

################## New Condition #################
new := proc (N, W, f, R, T) 
local ff, test, NN, WW, cd, KK, ww, i, cc, ccc, k, M, www,fff; 
#option trace:
ff := simplify(expand(f)); 
test := true; 
NN := N; 
KK := []; 
while test and NN <> [1] and ff <> 0 do 
    while type(LeadingCoefficient(ff, T), 'constant') = true and ff <> 0 do 
        KK := [op(KK), lt(ff, T)]; 
        ff := simplify(expand(ff-simplify(expand(lt(ff, T))))) 
    end do; 
    if RadicalMembership(LeadingCoefficient(simplify(ff), T), simplify(`<,>`(NN))) = true then 
        NN := [op({op(NN), LeadingCoefficient(ff, T)})]; 
        ff := simplify(ff-simplify(expand(lt(ff, T)))): 
    else 
        test := false:
    end if :
end do; 
NN := Basis(NN, R); 
WW := [op({seq(NormalForm(W[i], NN, R), i = 1 .. nops(W))})]; 
ff := op(nrm([NormalForm(ff, NN, R)])); 
if WW <> [] then 
    #`WW&Assign;nrm`(WW);
     WW:=nrm(WW):
    ww := NULL; 
    for i to nops(WW) do 
        if not type(WW[i], 'constant') = true then 
            ww := ww, WW[i] 
        end if 
    end do; 
    WW := [ww] 
end if; 
fff := ff; 
while type(LeadingCoefficient(fff, T), 'constant') = true and fff <> 0 do 
    fff := simplify(expand(fff-simplify(expand(lt(fff, T))))) 
end do; 
cc := fac([LeadingCoefficient(fff, T)]); 
ccc := {seq(op(nrm([k/LeadingCoefficient(k, R)])), `in`(k, cc))};
cd := `minus`({seq(expand(kk), `in`(kk, ccc))}, {-op(WW), op(WW)}); 
k := add(KK[j], j = 1 .. nops(KK)); 
M := [op(fac(WW, R))]; 
WW := [op({seq(ii/LeadingCoefficient(ii, R), `in`(ii, M))})]; 
www := [seq(nrm([ee]), `in`(ee, WW))]; 
if www <> [] then 
    www := [op(www)]; 
www := [seq(op(jj), `in`(jj, www))] 
end if; 
if nops(cd)<>0 then
     RETURN([{op([nrm([op(cd minus select(sff,cd))][1])])}, ff+k, NN, www]); 
else
     RETURN([{}, ff+k, NN, www]);
fi:
end proc:
##################
sf := proc (x) 
evalb(x <> 0): 
end proc: 

##############################
##       G2V Algorithm       ##
##############################
LM:=proc(f)
global tord;
if f<>0 then
    RETURN(LeadingMonomial(op(nrm([f])),tord));
fi:
RETURN(0);
end:
LC:=proc(f)
global tord;
RETURN(LeadingCoefficient(op(nrm([f])),tord));
end:
#############################
refine:=proc(LL,Jsig)
local u,f:
u:=Jsig:
f:=proc(Q)
global u:
if not divide(LM(Q[1]),u) then
    RETURN(true);
fi:
RETURN(false);
end:
RETURN(select(f,LL)):
end:
#############################
JPair2:=proc(P,Q,LMH)
#option trace;
global tord,Polys,Deg,NumOfFirstBuch,m,DDD;
local t,t1,t2,f,L,T:
####
f:=proc(r)
if r<>0 then 
    RETURN(true);
fi:
RETURN(false);
end:
####
T:=tord;
L:=select(f,LMH):
t:=lcm(P[3],Q[3]):
DDD:=max(DDD,degree(t));
t1,t2:=simplify(t/P[3]) , simplify(t/Q[3]):
if TestOrder(t1*LM(P[1]) , t2*LM(Q[1]),T) and NormalForm(t2*LM(Q[1]),L,T)<>0  then
    RETURN(expand([LeadingCoefficient(Polys[P[2]],T)*t2*Q[1],Q[2],t2*Q[3],t2*Q[4],P[2]]));
elif TestOrder(t2*LM(Q[1]) , t1*LM(P[1]),T) and NormalForm(t1*LM(P[1]),L,T)<>0 then 
    RETURN(expand([LeadingCoefficient(Polys[Q[2]],T)*t1*P[1],P[2],t1*P[3],t1*P[4],Q[2]]));
fi:
RETURN(NULL);
end:
#############################
tidy:=proc(P,Q)
global tord;
local T:
T:=tord;
if TestOrder(LM(P[1]),LM(Q[1]),T) then
    RETURN(true);
fi:
RETURN(false);
end:
#############################
RegRed:=proc(P)
#option trace:
local u1,v1,i,flag,Q,u2,v2,t,c,cmpt,T:
global tord,r,Grob,Polys,NumOfSup,R,GG,tp:
T:=tord;
u1,v1:=P[1],expand(P[4]*Polys[P[2]]):
i:=1:
cmpt:=0;
while i<=nops(R) and v1<>0 do
    flag:=false:
    Q:=R[i]:
    u2,v2:=Q[1],expand(Q[4]*Polys[Q[2]]):
    t:=LM(v1)/LM(v2):
    if divide(LM(v1),LM(v2)) and TestOrder(LM(expand(t*u2)) , LM(u1),T) and (P[2]<>Q[2] or cmpt=1) then# and Q[1]<>t*P[1] then
        cmpt:=1;
        c:=LeadingCoefficient(op(nrm([v1])),T)/LeadingCoefficient(op(nrm([v2])),T):
        if LM(u1)=t*LM(u2) then 
	       RETURN(1); 
	    fi;
        u1:=op(nrm([simplify(u1-c*t*u2)])): 
        v1:=SPolynomial(v1, v2, T);
        if LM(expand(t*u2)) = LM(u1) and c<>1 then
            u1:=op(nrm([u1/(1-c)])):
            v1:=op(nrm([v1/(1-c)])):
        fi:
        u1:=op(nrm([NormalForm(u1,Grob,T)])):
        v1:=op(nrm([NormalForm(v1,Grob,T)])):
        i:=1:
        flag:=true:
    else
        i:=i+1:
    fi:
od:
RETURN([op(nrm([u1])),op(nrm([v1]))]);
end:
############################
g2v:=proc(fff,GGG,N,W,RR,T,CPP,RRR)
#option trace:
local n,U,V,LMH,v,u,JP,NF,UV,Pm,g,flag,nr,m,P,k,Y,c,S,test,CP,f,fun,i,funi,G;
global tp,Grob,L,Polys,R,NumOfZero,NumOfF5,H,tord,null,notnull,fflag,iii,DDD,GG,rdz:
f:=fff[2];
CP:=CPP:
GG:=NULL:
tp:=RR;
R:=NULL:
fun:=proc (p, i) if i in {p[2],p[-1]} then RETURN(true) else RETURN(false) end if end proc;
for i from 1 to nops(GGG) do
    g:=NormalForm(GGG[i],N,RR):
    if g=0 then 
    CP:=remove(fun,CP,i);
    funi:=proc(p,j) if p<j then RETURN(p); else RETURN(p-1); fi;end:  
	CP:={seq([p[1],funi(p[2],i),p[3],p[4],funi(p[5],i)], p in CP)}:
    else
	   GG:=GG,g;
	   R:=R,[RRR[i][1],RRR[i][2],LM(g),RRR[i][4]];
    fi;
od;
G:=[GG];if N=[1] then RETURN(true,G,N,W,0,{},{}); ;fi;
R:=[R];
tord:=T;
null, notnull:=N,W:
Polys:=G:
m:=nops(Polys);
LMH:=[seq(LM(Grob[j]), j=1..nops(Grob))]:
g:=expand(NormalForm(f,N,RR)):
v:=op(nrm([NormalForm(g,Grob,T)])):
if v=0 and CP={} then 
    RETURN(true,G,N,W,0,{},{}):
elif v=0 then
    JP:=CP:
else
    
    u:=fff[1];
    Y := new(null, notnull, op(nrm([v])), RR, T);
                    c[dec] := Y[1];
                    S:=Y[2]: 
                    null := Y[3]; 
                    notnull := Y[4]; 
                    if S=0 then RETURN(true,G,N,W,0,{},{}): fi:
                    if c[dec] = {} then 
                        if S <> 0 then 
			                 v:=simplify(op(nrm([S])));
			                 nr:=nops(R):
			                 R:=[op(R),[u,nr+1,LM(v),1]]:
			                 Polys:=[op(Polys),v]:
			                 JP:=[op({op(CP),seq(JPair2(R[-1],R[j],LMH),j=1..nr)})]:
			                 JP:=sort(JP,tidy);
			             end if: 
                    elif S <> 0 then 
                        RETURN(false,Polys,null,notnull,[u,S],CP,R);
                    end if: 
fi:
while nops(JP)<>0 do #print(nops(JP));
        P:=JP[1]:
        JP:=JP[2..-1]:
        NF:=RegRed(P):
        if NF<>1 then
            u:=NF[1]:
            v:=NF[2]:
            if v=0 then
                #H:=[op(H),u]:
                LMH:=[op(LMH),LM(u)]:
                JP:=refine(JP,LM(u)):
                rdz:=rdz+1;
            else
                    v:=NormalForm(v,N,RR);
                    ##print(null, notnull, op(nrm([v])), RR, T);
                    Y := new(null, notnull, op(nrm([v])), RR, T);
                    c[dec] := Y[1];
                    S:=Y[2]: 
                    null := Y[3]; 
                    notnull := Y[4]; 
                    if c[dec] = {} then 
                        if S <> 0 then 
			                 v:=simplify(op(nrm([S])));
			                 nr:=nops(R):
			                 R:=[op(R),[u,nr+1,LM(v),1]]:
			                 Polys:=[op(Polys),v]:
			                 JP:=[op({op(JP),seq(JPair2(R[-1],R[j],LMH),j=1..nr)})]:
			                 JP:=sort(JP,tidy);
			             end if: 
                    elif S <> 0 then 
                        RETURN(false,Polys,null,notnull,[u,S],JP,R);
                    end if: 
            fi:
        fi:
od:
RETURN(true,Polys,null,notnull,0,{},{});
end:
#######################################
#######################################
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,NnN,n; 
#option trace;
N := Basis(LL,R); 
W := M; 
WW := seq(NormalForm(W[i], N, R), i = 1 .. nops(W));
WWW := seq(factor([WW][i]), i = 1 .. nops([WW])); 
test := true; 
t:= true; 
NN := N; 
h := product(W[i], i = 1 .. nops(W)); 
if RadicalMembership(h, `<,>`(NN)) = true then
    test := false;
    NN := {1};
    RETURN([test, NN,[WW]]):
end if; 
while t do
    t := false;
    facN := [seq([op(i)], `in`(i, [seq(fac([NN[i]]), i = 1 .. nops(NN))]))];
    facW := [op(fac([WWW], R))];
    NNN := [];
    L := NULL;
    for i from 1 to nops(facN) do
        for x in facN[i] do
           flag := true;
           for w in facW while flag do
             if divide(w, x) then
                flag := false;
                L := L, x:
             end if:
           end do:
           
        end do;
        NnN[i]:= [op(`minus`({op(facN[i])}, {L}))];
        n[i]:=simplify(expand(product(NnN[i][j], j = 1 .. nops(NnN[i]))));
     end do;   
     NNNN :=[seq(n[i],i=1..nops(facN))];
    if {op(NNNN)} <> {op(NN)} then 
        t := true;
        NN := Basis(NNNN, R);
        WW := seq(NormalForm([WWW][i], NN, R), i = 1 .. nops([WWW]));
        WWW := seq(factor([WW][i]), i = 1 .. nops([WW]))
    end if:
end do; 
ww := NULL; 
for i to nops([WWW]) do
    if type(WWW[i], 'constant') = true then
    else
        ww := ww, [WWW][i]:
    end if
end do; 
WWW := ww; 
MM := [op(fac([WWW], R))]; 
WWW :=[op({seq(i/LeadingTerm(i, R)[1], `in`(i, MM))})]; 
www :=[seq(nrm([ee]), `in`(ee, WWW))]; 
if www <> [] then
    www := [op(www)];
    www := [seq(op(jj), `in`(jj, www))]:
end if; 
RETURN([test, NN, www]): 
end proc:

##########################################
branch := proc (v, ffff,BBB, NNN, WWW, R, T,CPP,RR) 
local cd,l,ll,i,ff,Y,X,BB,pivot,test,TTT,Bb,g,TT,vv,LL,gg,B,CN,N,W,NN,WW,f,CP,RRRR;
global LIST, Grob,termord1, termord2, termord3, iii, DDD,ind,rdz; 
#option trace:
f:=ffff[2];
CN := canspec(NNN, WWW, R); 
B:=BBB:
N:=CN[2];
W:=CN[3];
cd := []; 
Bb:=B: 
Y := new(N, W, f, R, T);
cd := [op(Y[1])];
ff := simplify(expand(Y[2])); 
NN := Y[3]; 
WW := Y[4];  
if v = [] then 
    vv := NULL: 
else 
    vv := op(v): 
end if;  
TT[vv] := [v, Bb, NN, WW]; 
if cd = [] then 
    X:= g2v([ffff[1],ff],Bb, NN, WW, R, T,CPP,RR); 
    test := X[1]; 
    BB := X[2]; 
    NN := X[3]; 
    WW := X[4]; 
    ff:=X[5];
    CP:=X[6];
    RRRR:=X[7]:
    if test  then 
        #TT[vv] := [[vv], nrm(InterReduce(BB,T)), NN, WW]; 
        TT[vv] := [[vv], nrm(BB), NN, WW]; 
	if CP={} then
	    LIST := LIST, TT[vv]:
        
	fi;
    else 
        branch([vv],ff, BB, NN, WW, R, T,CP,RRRR); 
    end if: 
else 
    branch([0,vv], [ffff[1],ff],Bb, [op(N),op(cd)], W, R, T,CPP,RR); 
    branch([1,vv], [ffff[1],ff],Bb, N, [op(W),op(cd)], R, T,CPP,RR);
end if; 
end proc:
################## Incremental Montes Algorithm ##############
IncMontes:=proc(f,Grobsys,R,T)
local L,i,G,N,W,ff,CPP;
global LIST, DDD,rdz,Grob;
#option trace;
L:=NULL:
for i from 1 to nops(Grobsys) do 
    CPP:={};
    LIST:=NULL:
    G:=Grobsys[i];
    Grob:=G[2]:
    N:=G[3]:
    W:=G[4]:
    ff:=NormalForm(f,N,R);
    if ff=0 then
        LIST:=LIST,[G[1],Grob,G[3],G[4]]: 
    else 
	 branch([], [1,ff],Grob, N, W, R, T, {},[seq([0,i,LM(Grob[i]),1],i=1..nops(Grob))]);
    fi:
    L:=L,LIST;
od:
LIST:=L;
RETURN(LIST);
end:

#######################
selpoly1 := proc (f, g) 
local ZPf, ZPg;
global termord1, termord2; 
ZPf := LeadingCoefficient(f, termord2); 
ZPg := LeadingCoefficient(g, termord2); 
RETURN(TestOrder(ZPg, ZPf, termord1)): 
end proc:
#######################
################## F5Montes Algorithm ##############
Montes:=proc(FF,R,T)
local f,t1,t2,b1,b2,F;
global LIST,DDD,ind,rdz,Grob,termord1, termord2;
local IndPolys,POLYS;
#option trace;
t1,b1:=kernelopts(cputime,bytesused);
termord1, termord2:=R,T:
IndPolys:=[seq(i,i=1..nops(FF))];
F:=InterReduce(FF,prod(T,R));
POLYS:=FF;
rdz:=0;
DDD:=0;
LIST:=NULL;
f:=F[1];
Grob:=[];
ind:=nops(F)-1;
branch([], [1,f],[], [], [], R, T,{},[]);
for f in F[2..-1] do
    ind:=ind-1;print(ind);
    IncMontes(f,[LIST],R,T); 
od;
t2,b2:=kernelopts(cputime,bytesused);
printf("%-1s %1s %1s%1s    : %3a %3a\n",The, cpu, time, is,t2-t1,(sec)):
 printf("%-1s %1s %1s   : %3a %3a\n",The,used,memory,b2-b1,(bytes)):
 printf("%-1s %1s %1s        : %3g\n",Num,of,Rds,rdz):
 printf("%-1s %1s %1s        : %3a\n",Num,of,sys,nops([LIST])):
 printf("%-1s %1s %1s %1s    : %3a\n",The, max, deg, is,DDD):
RETURN(LIST);
end:




#############P3P from Kapur paper#########
#F:=subs(b=1,a=1,[(1 - a)*y^2 - a*x^2 - p*y + a*r*x*y + 1, (1 - b)*x^2 - b*y^2 - q*x + b*r*x*y + 1]);
#Montes(F,plex(p,q,r,a,b),plex(x, y));
F:=[x^2 + y^2 - p^2,
(a - 1)*y + b*(1 - x),
-a*y + b*(x + p),
(1 - x)^2 + y^2 - (1 - r)^2,
m*y - x*n,
(1 - m)*y + (x + r - 2)*n,
a^2 + b^2 - (1 - m)^2 + n^2];
                      #############################Examples#################################### 
print("********B2*******");
Montes([a*y^2+b*x^3+c, 2*a*y, 3*b*x^2], plex(a, b, c), plex(x, y)):   
print("********B3*******");
Montes([a*y^2+b*x^2+c*x^3+d, 2*a*y, 3*b*x+3*c*x^2], plex(a, b, c,d), plex(x, y)):                       
print("********B4*******");
Montes([x+c*y+b*z+a, c*x+y+a*z+b, b*x+a*y+z+c], plex(a, b, c), plex(x, y, z)):       
print("********B5*******");
Montes([(a-1)*y[2]-b*(x[2]-1),(a-1)*(x[2]+1)+b*y[2],(a+1)*y[3]-b*(x[3]+1),(a+1)*(x[3]-1)+b*y[3],(x[3]-a)^2+y[3]^2-(x[2]-a)^2-y[2]^2 ], plex(a,b), plex(x[2],x[3],y[2],y[3])):
print("********ex3*******");
Montes([x^2+b*y^2+2*c*x*y+2*d*x+2*e*y+f, x+c*y+d, b*y+c*x+e], plex(b, c, d, e, f), plex(x, y)):
print("********F4*******");
Montes([a*x^3*y+c*x*y^2, x^4*y+3*d*y, c*x^2+b*x*y, x^2*y^2+a*x^2, x^5+y^5], plex(a, b, c, d), plex(x, y)):
print("********B6*******");
Montes([a*x^2*y^2+b*x*z^2, b*x*y^2+c*x^2+2, c*x^2+b*y^2*z], plex(a, b, c), plex(x, y, z)):
print("********B7*******");
Montes([x^3-a*x*y, x^2*y-2*y^2+x], plex(a), plex(x, y)):
print("********S7*******");
Montes([y^2-z*x*y+x^2+z-1, x*y+z^2-1, y^2+x^2+z^2-r^2], plex(r), plex(x, y, z)):
print("********ex1*******");
Montes([a*x^2*y+b*x*y, b*x+y^2], plex(a, b), plex(x, y)):
print("********ex2*******");
Montes([2599*e2^2-400*e2*e3+80*e2*f3-80*e3*f2+2599*f2^2-400*f2*f3-2200*e2-240*f2-20*Q2, 14-12*e2-110*f2-2*e3-10*f3-P1, 2397-2200*e2+240*f2-200*e3+40*f3-20*Q1, 16*e2^2-4*e2*e3-20*e2*f3+20*e3*f2+16*f2^2-4*f2*f3-12*e2+110*f2-P2], plex(P1, Q1, P2, Q2), plex(e2, f2, e3, f3)):
print("********ex3*******");
Montes([x^2+b*y^2+2*c*x*y+2*d*x+2*e*y+f, x+c*y+d, b*y+c*x+e], plex(b, c, d, e, f), plex(x, y)):
print("********ex4*******");
Montes([r-c1+l*(s1*s2-c1*c2), z-s1-l*(s1*c2+s2*c1), s1^2+c1^2-1, s2^2+c2^2-1], plex(r, z, l), plex(c1, s1, c2, s2)):
print("********ex6*******");
Montes([a*x*z+b*x*z+a, (a^2+a)*x*y], plex(a, b), plex(x, y, z)):
print("********ex7*******");
Montes([a+b^2, y+1, b*x+1, a*x-b], plex(a, b), plex(x, y)):
print("********ex8*******");
Montes([x+c*y+b*z+a, c*x+y+a*z+b, b*x+a*y+z+c], plex(a, b, c), plex(x, y, z)):
print("********F1*******");
Montes([a*x^4*y+x*y^2+b*x, x^3+2*x*y, b*x^2+x^2*y], plex(a, b), plex(x, y)):
print("********F2*******");
Montes([a*x^2*y^3+b*y+y, x^2*y^2+x*y+2, a*x^2+b*y+2], plex(a, b), plex(x, y)):
print("********F3*******");
Montes([b*x^3+x^2+2, c*x^2+d*x, a*x^4+c*x^2+b], plex(a, b, c, d), plex(x)):
print("********F4*******");
Montes([a*x^3*y+c*x*y^2, x^4*y+3*d*y, c*x^2+b*x*y, x^2*y^2+a*x^2, x^5+y^5], plex(a, b, c, d), plex(x, y)):
print("********F5*******");
#Montes([a*x^2*y+b*x+y^3, a*x^2*y+b*x*y, y^2+b*x^2*y+c*x*y], plex(a, b, c), plex(x, y)):
print("********F6*******");
Montes([a*x^2*y^2+b*x*z^2, b*x*y^2+c*x^2+2, c*x^2+b*y^2*z], plex(a, b, c), plex(x, y, z)):
print("********F7*******");
Montes([x^4+a*x^3+b*x^2+c*x+d, 4*x^3+3*a*x^2+2*b*x+c], plex(a, b, c, d), plex(x)):
print("********S1*******");                         
Montes([a*(x+y), b*(x+y), x^2+a*x], plex(a, b), plex(x, y)):
print("********S2*******");
Montes([x1^2, x1*x2, x1*x3^2, x1*a+x2, x2*x3-x3^2, x2*a, x3^3, x3^2*a, a^2], plex(a), plex(x1, x2, x3)):
print("********S3*******");
Montes([x^3-a*x*y, x^2*y-2*y^2+x], plex(a), plex(x, y)):
print("********S4*******");
Montes([a*x+y-1, b*x+y-2, 2*x+a*y, b*x+a*y+1], plex(a, b), plex(x, y)):
print("********S5*******");
Montes([x4-a4+a2, x1+x2+x3+x4-a1-a3-a4, x1*x3*x4-a1*a3*a4, x1*x3+x1*x4+x2*x3+x3*x4-a1*a4-a1*a3-a3*a4], plex(a1, a2, a3, a4), plex(x1, x2, x3, x4)):
print("********S6*******");
Montes([v*x*y+u*x^2+x, u*y^2+x^2], plex(u, v), plex(x, y)):
print("********S7*******");
Montes([y^2-z*x*y+x^2+z-1, x*y+z^2-1, y^2+x^2+z^2-r^2], plex(r), plex(x, y, z)):
print("********S8*******");
Montes([a-b+(x*y*a-x^2*y*b-3*a)^3+(x*y*b-3*x*b-5*b)^4, x*y*a-x^2*y*b-3*a, x*y*b-3*x*b-5*b], plex(a, b), plex(x, y)):
print("********S9*******");
Montes([x+c*y+b*z+a, c*x+y+a*z+b, x*b+a*y+z+c], plex(a, b, c), plex(x, y, z)):
print("********S10*******");
Montes([a-l3*c3-l2*c1, b-l3*s3-l2*s1, c1^2+s1^2-1, c3^2+s3^2-1], plex(a, b, c), plex(c1, c3, l1, l2, l3, s1, s3)):
print("********S15*******");
Montes([a+d*s1, b-d*c1, l2*c2+l3*c3-d, l2*s2+l3*s3-c, c1^2+s1^2-1, c3^2+s3^2-1], plex(a, b, c, d), plex(c1, c2, c3, l1, l2, l3, s1, s2, s3)):
print("********S16*******");
Montes([w12+w14, w12+w13, w12+w15, w12+w23+w25-w26*x+w26, w12+w25-w26*y+w26, w12+w23-w26*z+w26, w23+w34+x*w36, w13+w34-w36*y+w36, w23+z*w36, w14+w34+w45-w46*x+w46, w34+y*w46, w45+z*w56, w15+w45-z*w56+w56, -w26+w26*x+x*w36-w46+w46*x+w56*x, -w26+w26*y-w36+w36*y+y*w46+w56*y, -w26+w26*z+z*w36+w46*z-w56+z*w56], plex(x, y, z), plex(w12, w13, w14, w15, w23, w25, w34, w45, w26, w36, w46, w56)):
print("********S13*******");
Montes([a*x^2*y+a+3*b^2, a*(b-c)*x*y+a*b*x+5*c], plex(a, b, c), plex(x, y)):                                      
print("********B1*******");
#Montes([a*x^4+b*x^2*y^2-2*c*x^2*y-d*x*y^2+e*y^2+f, 4*a*x^3+2*b*x*y^2-4*c*x*y-d*y^2, 2*b*x^2*y-2*c*x^2-2*d*x*y+2*e*y], plex(a, b, c, d, e, f), plex(x, y)): 
####################################################
########			             ####### 
########    Improved Montes's algorithm      #######
########                                     #######
####################################################

															     															     
#########Division algorithm ############
Div:=proc(f,L,T)															     
p:=f;
r:=0;
while p<>0 do
    i:=1:
    flag:=false:
    while i<=nops(L) and not flag do
        if divide(LeadingMonomial(p,T),LeadingMonomial(L[i],T)) then
            p:=SPolynomial(p,L[i],T);
            flag:=true;
        else
            i:=i+1;
        fi:
    od:
    if not flag then
        r:=r+LeadingCoefficient(p,T)*LeadingMonomial(p,T);
        p:=simplify(p-LeadingCoefficient(p,T)*LeadingMonomial(p,T));
    fi:
od:
RETURN(r);
end:            
###PRIMEPART Algorithm###################
ppart := proc (K, T)
local h, i, ti, o, oo, ooo;
if {op(K)} <> {0}and {op(K)} <> {} then
    h := LeadingTerm(K[1], T)[1];
    for i from 2to nops(K) do
        ti := LeadingTerm(K[i], T)[1];
        h := gcd(h, ti)
    end do;
    RETURN(simplify(expand(K/h)))
else
    RETURN(K)
end if
end proc:
####lt- Algorithm######################
lt := proc (f, T) 
RETURN(LeadingTerm(f, T)[1]*LeadingTerm(f, T)[2]): 
end proc:                      
####Normalize Algorithm################
nrm := proc (F) 
local KK, h, i, NM: 
if F = [] then 
    RETURN(F) 
end if; 
KK := F; 
h := denom(KK[1]); 
for i from 2 to nops(KK) do 
    h := simplify(lcm(denom(KK[i]), h)) 
end do; 
NM := simplify(expand(h*KK)); 
RETURN(NM): 
end proc:
#####convert a polynomial to the list of it terms############                                  
FL := proc (f, T) 
local L, p; 
L := []; 
p := f; 
while p <> 0 do 
    L := [op(L), lt(p, T)]; 
    p := simplify(p-lt(p, T)) 
end do; 
RETURN(L) 
end proc:
####PRIMITIVATION one list 'G' of polynomials################  
plp := proc (G, T) 
local FFFF, FFF, FF, FFi, Gi, i, j, gj, k; 
FFFF := []; 
FFF := []; 
for i to nops(G) do 
    Gi := FL(G[i], T); 
    FFi := ppart(Gi, T); 
    FFF := [op(FFF), FFi] 
end do; 
for j to nops(FFF) do 
    gj := 0; 
    for k to nops(FFF[j]) do 
        gj := gj+FFF[j][k] 
    end do; 
    FFFF := [op(FFFF), gj] 
end do; 
RETURN(FFFF) 
end proc:
####GGE ALGORITM###############################
gge := proc (F, T) 
RETURN(InterReduce(F, T)): 
end proc:
####convert a list of polynomial to the set "FACVAR"##########
fac := proc (L, T) 
local A, AA, AAA, N, P, p, B, C, i; 
A := L; 
AA := [seq(factor(i), `in`(i, A))]; 
AAA := NULL; 
B := []; 
for i to nops(AA) do 
    p := AA[i]; 
    if irreduc(p) = true then 
        AAA := AAA, p: 
    else 
        AAA := AAA, op(convert(p, list)): 
    end if 
end do; 
P := NULL; 
for i to nops([AAA]) do 
    if `not`(type(AAA[i], 'constant')) then 
        P := P, [AAA][i]: 
    end if 
end do; 
N := NULL; 
for i to nops([P]) do 
    if irreduc([P][i]) = false and type([P][i], 'constant') <> true then 
        N := N, op([P][i])[1]: 
    else 
        N := N, [P][i]: 
    end if 
end do; 
RETURN({N}): 
end proc:
####CANSPEC ALGORITM############################  
canspec := proc (LL, M, R) 
local N, W, WW, WWW, NN, NNN, test, i,h, t, facN, facW, x, NNNN, L, w, flag, ww, www, MM; 
#option trace;
N := Basis(LL,R); 
W := M; 
WW := seq(NormalForm(W[i], N, R), i = 1 .. nops(W));
WWW := seq(factor([WW][i]), i = 1 .. nops([WW])); 
test := true; 
t:= true; 
NN := N; 
h := product(W[i], i = 1 .. nops(W)); 
if RadicalMembership(h, `<,>`(NN)) = true then
    test := false;
    NN := {1};
    RETURN(test, NN):
end if; 
while t do
    t := false;    
    facN := [seq([op(i)], `in`(i, [seq(fac([NN[i]]), i = 1 .. nops(NN))]))];####
    facW := [op(fac([WWW], R))];
    NNN := [];
    L := NULL;
    for i from 1 to nops(facN) do
        for x in facN[i] do
           flag := true;
           for w in facW while flag do
             if divide(w, x) then
                flag := false;
                L := L, x:
             end if:
           end do:
        end do;
        NnN[i]:= [op(`minus`({op(facN[i])}, {L}))];
        n[i]:=simplify(expand(product(NnN[i][j], j = 1 .. nops(NnN[i]))));
     end do;   
     NNNN :=[seq(n[i],i=1..nops(facN))];
     if {op(NNNN)} <> {op(NN)} then
        t := true;
        NN := Basis(NNNN, R);
        WW := seq(NormalForm([WWW][i], NN, R), i = 1 .. nops([WWW]));
        WWW := seq(factor([WW][i]), i = 1 .. nops([WW]))
    end if:
end do; 
ww := NULL; 
for i to nops([WWW]) do
    if type(WWW[i], 'constant') = true then
    else
        ww := ww, [WWW][i]:
    end if
end do; 
WWW := ww; 
MM := [op(fac([WWW], R))]; 
WWW :=[op({seq(i/LeadingTerm(i, R)[1], `in`(i, MM))})]; 
www :=[seq(nrm([ee]), `in`(ee, WWW))]; 
if www <> [] then
    www := [op(www)];
    www := [seq(op(jj), `in`(jj, www))]:
end if; 
RETURN([test, NN, www]): 
end proc:
####NEW COND ALGORITM##########################
new := proc (N, W, f, R, T) 
local ff, test, NN, WW, cd, KK, ww, i, cc, ccc, k, M, www,fff; 
#option trace:
ff := simplify(expand(f)); 
test := true; 
NN := N; 
KK := []; 
while test and NN <> [1] and ff <> 0 do 
    while type(LeadingCoefficient(ff, T), 'constant') = true and ff <> 0 do 
        KK := [op(KK), lt(ff, T)]; 
        ff := simplify(expand(ff-simplify(expand(lt(ff, T))))) 
    end do; 
    if RadicalMembership(LeadingCoefficient(simplify(ff), T), simplify(`<,>`(NN))) = true then 
        NN := [op({op(NN), LeadingCoefficient(ff, T)})]; 
        ff := simplify(ff-simplify(expand(lt(ff, T)))): 
    else 
        test := false:
    end if :
end do; 
NN := Basis(NN, R); 
WW := [op({seq(NormalForm(W[i], NN, R), i = 1 .. nops(W))})]; 
ff := op(nrm([NormalForm(ff, NN, R)])); 
if WW <> [] then 
     WW:=nrm(WW):
    ww := NULL; 
    for i to nops(WW) do 
        if not type(WW[i], 'constant') = true then 
            ww := ww, WW[i] 
        end if 
    end do; 
    WW := [ww] 
end if; 
fff := ff; 
while type(LeadingCoefficient(fff, T), 'constant') = true and fff <> 0 do 
    fff := simplify(expand(fff-simplify(expand(lt(fff, T))))) 
end do; 
cc := fac([LeadingCoefficient(fff, T)]); 
ccc := {seq(op(nrm([k/LeadingCoefficient(k, R)])), `in`(k, cc))};
cd := `minus`({seq(expand(kk), `in`(kk, ccc))}, {-op(WW), op(WW)}); 
k := sum(KK[j], j = 1 .. nops(KK)); 
M := [op(fac(WW, R))]; 
WW := [op({seq(ii/LeadingCoefficient(ii, R), `in`(ii, M))})]; 
www := [seq(nrm([ee]), `in`(ee, WW))]; 
if www <> [] then 
    www := [op(www)]; 
www := [seq(op(jj), `in`(jj, www))] 
end if; 
if nops(cd)<>0 then
     RETURN([{op([nrm([op(cd minus select(sff,cd))][1])])}, ff+k, NN, www]); 
else
     RETURN([{}, ff+k, NN, www]);
fi:
end proc:
#####Select Algoritm###################### 
sff := proc (x) 
evalb(type(x, 'constant')): 
end proc:
#####Select Algoritm###################### 
sf := proc (x) 
evalb(x <> 0): 
end proc: 
######Normal Sterategy######################
f2 := proc (p, q) 
local t; 
global setB, Tord; 
t := Tord; 
RETURN(TestOrder(lcm(lt(setB[p[1]], t), lt(setB[p[2]], t)), lcm(lt(setB[q[1]], t), lt(setB[q[2]], t)), t)) 
end proc:
######Update Algoritm########################
update := proc (h, LL, BB) 
local L1, L, L2, L3, L4, L5, i, lth, flag, j, p, E,nn; 
global setB, Tord; 
#option trace;
lth := LeadingTerm(setB[h], Tord)[2]; 
L := LL; 
L1 := NULL; 
L2 := NULL; 
for i to nops(L) do 
    if gcd(LeadingTerm(setB[i], Tord)[2], lth) = 1 then 
        L1 := L1, i: 
    else 
        L2 := L2, i: 
    end if: 
end do; 
L1 := [L1]; 
nn:=nops(L1):
L2 := [L2]; 
for i in L2 do 
    flag := false;
    member(i,L2,'zz');
    L2:=subsop(zz=NULL,L2); 
    for j in L1 while not flag do 
        if divide(lcm(lth, LeadingTerm(setB[i], Tord)[2]), LeadingTerm(setB[j], Tord)[2]) then 
            flag := true: 
        end if: 
    end do; 
    for j in L2 while not flag do 
        if divide(lcm(lth, LeadingTerm(setB[i], Tord)[2]), LeadingTerm(setB[j], Tord)[2]) then 
            flag := true: 
        end if: 
    end do; 
    if not flag then 
        L1 := [op(L1), i]: 
    end if: 
end do; 
L5 := NULL; 
for p in BB do 
    if not divide(lcm(LeadingTerm(setB[p[1]], Tord)[2], LeadingTerm(setB[p[2]], Tord)[2]), lth) or lcm(LeadingTerm(setB[p[1]], Tord)[2], LeadingTerm(setB[p[2]], Tord)[2]) = lcm(LeadingTerm(setB[p[1]], Tord)[2], lth) or lcm(LeadingTerm(setB[p[1]], Tord)[2], LeadingTerm(setB[p[2]], Tord)[2]) = lcm(LeadingTerm(setB[p[2]], Tord)[2], lth) then 
        L5 := L5, p 
    end if: 
end do; 
if BB=[] and nops(LL)=1 then
    RETURN([{h,op(LL)}]);
fi;
E := [seq({i, h}, `in`(i, L1[nn+1..nops(L1)])), L5]; 
RETURN(sort(E, f2)): 
end proc:
#####CondPGB algorithm ####################
condpgb := proc (B, N, W, R, T,CPP)
local test, NN, WW, BB, i, j, t, s, J, x, c, Y, S, ii;
global setB, Tord, DDD, iii,inde;
#option trace;
felag:=true;
CP:=CPP:
inde:=true;
test := true;
s := nops(B);
if s=0 then RETURN([true, [], N, W,1,[]]);fi;
setB := B;
Tord := T;
BB := B;
NN := N;
WW := W;
t := s;
if CP={} then  
    for tt from 1 to nops(B) do
	CP:=update(tt, [seq(i, i = 1 .. tt-1)], CP):
   od:
else
    CP:=update(t, [seq(i, i = 1 .. t-1)], CP):
fi:    
J := sort(CP, f2);
while nops(J) <> 0 and test do
    x := J[1];
    J := J[2 .. -1];
    i := x[1];
    j := x[2];
    S := expand(simplify(Div(op(nrm([SPolynomial(BB[i], BB[j], T)])), select(sf, BB), T)));
    S := simplify(NormalForm(op(nrm([S])), NN, R));
    if S<>0 then 
       felag:=false;
    fi;
    DDD := max(DDD, degree(lcm(BB[i], BB[j]), {op(T)}));
    if S <> 0 then
        Y := new(NN, WW, S, R, T);
        c[dec] := Y[1];
        S := Y[2];
        NN := Y[3];
        WW := Y[4];
        if c[dec] = {} then
            if S <> 0 then
                t := t+1;
                BB := [op(BB), S];
                setB := BB;
                J := update(t, [seq(i, i = 1 .. t-1)], J): 
            end if
        elif S <> 0 then
            test := false;
            BB := [op(BB), S]:
        end if:
    else
        iii := iii+1:
    end if:
end do;
if test then 
    BB := InterReduce(BB, prod(T, R)):
end if;
if felag=true then 
   RETURN([test, B, NN, WW,S,J]):
else 
    RETURN([test, BB, NN, WW,S,J]):
fi;  
end proc:
#####Selection strategy for polynomials############
selpoly1 := proc (f, g) 
local ZPf, ZPg;
global termord1, termord2; 
ZPf := LeadingCoefficient(f, termord2); 
ZPg := LeadingCoefficient(g, termord2); 
RETURN(TestOrder(ZPg, ZPf, termord1)): 
end proc:
######Branch algorithm ###############################
branch := proc (v, BBB, N, W, R, T,CPP) 
local cd, l, ll, i, f, ff, Y, X, BB, NN, WW, pivot, test, TTT, Bb, g, TT, vv, LL,gg,B; 
global LIST, termord1, termord2, termord3, iii, DDD,inde; 
#option trace:
cd := []; 
termord1 := R; 
termord2 := T; 
termord3 := prod(termord2, termord1); 
B:=BBB;
Bb := B; 
if inde then 
    Bb:=nrm([seq(NormalForm(h,N,R), h in Bb)]);
    CP:=CPP;
    f := B[-1];
    Y := new(N, W, f, R, T);
    if Y[2]=0 then 
	    fun:=proc (p, j) if j in p then RETURN(true) else RETURN(false) end if end proc;
	    CP:=remove(fun,CP,nops(B));
	    funi:=proc(p,j)  if p<j then RETURN(p); else RETURN(p-1); fi;end:  
	    CP:={seq({funi(p[1],nops(B)), funi(p[2],nops(B))}, p in CP)}:
    fi:
	cd := [op(Y[1])]; 
	ff := simplify(expand(Y[2])); 
	NN := Y[3]; 
	WW := Y[4]; 
	Bb[-1] := Y[2]; 
	if cd <> [] then 
	    pivot := nops(B): 
	end if: 
else
    CP:=CPP;
    for i from nops(B) to 1 by -1 while cd = [] do 
	f := B[i];
	Y := new(N, W, f, R, T);
	if Y[2]=0 then 
	    fun:=proc (p, j) if j in p then RETURN(true) else RETURN(false) end if end proc;
	    CP:=remove(fun,CP,i);
	    funi:=proc(p,j)  if p<j then RETURN(p); else RETURN(p-1); fi;end:  
	    CP:={seq({funi(p[1],i), funi(p[2],i)}, p in CP)}:
	fi:
	cd := [op(Y[1])]; 
	ff := simplify(expand(Y[2])); 
	NN := Y[3]; 
	WW := Y[4]; 
	Bb[i] := Y[2]; 
	if cd <> [] then 
	    pivot := i: 
	end if: 
    end do;
fi: 
if nops(B)=0 then
    NN := N; 
    WW := W;
fi:
gg := Bb[pivot]; 
LL := NULL; 
for i to nops(Bb) do 
    if Bb[i] <> 0 then 
        LL := LL, Bb[i]: 
	fun:=proc (p, j) if j in p then RETURN(true) else RETURN(false) end if end proc;
	    CP:=remove(fun,CP,i);
	    funi:=proc(p,j)  if p<j then RETURN(p); else RETURN(p-1); fi;end:  
	    CP:={seq({funi(p[1],i), funi(p[2],i)}, p in CP)}:	
    end if: 
end do; 
Bb := [LL];
if v = [] then 
    vv := NULL: 
else 
    vv := op(v): 
end if; 
TT[vv] := [v, Bb, NN, WW]; 
if cd = [] then 
    X := condpgb(Bb, NN, WW, R, T,CP); 
    test := X[1]; 
    BB := X[2]; 
    NN := X[3]; 
    WW := X[4]; 
    CP:=X[6];
    if test then 
        TT[vv] := [[vv], BB, NN, WW]; 
        LIST := LIST, TT[vv]: 
    else 
        branch([vv], BB, NN, WW, R, T,CP) 
    end if: 
else 
    newvertex(1, [vv], cd, Bb, NN, WW, pivot, R, T, gg,CP); 
    newvertex(0, [vv], cd, Bb, NN, WW, pivot, R, T, gg,CP): 
end if; 
RETURN([LIST]): 
end proc: 
#######NewVertex algorithm ##################
newvertex := proc (n, u, cd, B, N, W, pivot, R, T, gg,CPP) 
local i, r,v, WW, NN, BB, cond, CN, TT, vv; 
#option trace;
if u = 0 then
    vv := [NULL]
else
    vv := op(u), n
end if; 
BB := simplify(B); 
r := product([op(cd)][i], i = 1 ..nops(cd)); 
if n = 0 then
    cond := r;
    WW := W;
    NN := Basis([op(cd), op(N)], R);
else
    cond := r;
    WW := {op(cd), op(W)};
    
    NN := N;
    BB := B;
    if new(NN, WW, lt(gg, T), R, T)[1]={} then 
    for i to nops(B) do
         if B[i] <> gg and divide(LeadingMonomial(B[i], T), LeadingMonomial(gg, T)) = true then
           BB := subs(B[i] = expand(simplify(SPolynomial(B[i], gg, T))), BB)
         end if
    end do
    end if;
end if;
CN := canspec(NN, WW, R); 
CP:=CPP:
if CN[1] then
LLL:=NULL:
    TT[vv] := [cond, BB, NN, nrm(WW)]; 
    if vv = NULL then
        v := 0;
          branch(v, BB, CN[2], CN[3], R, T,CP):#######
    else
         branch([vv], BB, CN[2], CN[3], R, T,CP):##########
    end if:
else
    RETURN([vv, BB, {1}, WW])
end if: 
end proc:

######Impeoved DisPGB algorithm###########
DISPG := proc (F, R, T) 
local B,t,q; 
global LIST, iii, DDD,NR,NF,termord1,termord2,termord3,inde;
t1,b1:=kernelopts(cputime,bytesused);
NR:=0:
NF:=0: 
LIST := NULL; 
iii := 0; 
DDD := 0;
inde:=false;
termord1 := R; 
termord2 := T; 
termord3 := prod(termord2, termord1);  
B := sort(InterReduce(F, prod(T, R)), selpoly1);
J:={};
branch([], B, [], [], R, T, J):
t2,b2:=kernelopts(cputime,bytesused);
printf("%-1s %1s %1s%1s    : %3a %3a\n",The, cpu, time, is,t2-t1,(sec)):
 printf("%-1s %1s %1s   : %3a %3a\n",The,used,memory,b2-b1,(bytes)):
 printf("%-1s %1s %1s %1s    : %3a\n",The, num, red, is,iii):
 printf("%-1s %1s %1s %1s    : %3a\n",The, num, sys, is,nops([LIST])):
  printf("%-1s %1s %1s %1s    : %3a\n",The, max, deg, is,DDD):
RETURN(LIST): 
end proc:    

##################
##################
                      #############################Examples#################################### 
print("********B2*******");
DISPG([a*y^2+b*x^3+c, 2*a*y, 3*b*x^2], plex(a, b, c), plex(x, y)):   
print("********B3*******");
DISPG([a*y^2+b*x^2+c*x^3+d, 2*a*y, 3*b*x+3*c*x^2], plex(a, b, c,d), plex(x, y)):                       
print("********B4*******");
DISPG([x+c*y+b*z+a, c*x+y+a*z+b, b*x+a*y+z+c], plex(a, b, c), plex(x, y, z)):       
print("********B5*******");
DISPG([(a-1)*y[2]-b*(x[2]-1),(a-1)*(x[2]+1)+b*y[2],(a+1)*y[3]-b*(x[3]+1),(a+1)*(x[3]-1)+b*y[3],(x[3]-a)^2+y[3]^2-(x[2]-a)^2-y[2]^2 ], plex(a,b), plex(x[2],x[3],y[2],y[3])):
print("********ex3*******");
DISPG([x^2+b*y^2+2*c*x*y+2*d*x+2*e*y+f, x+c*y+d, b*y+c*x+e], plex(b, c, d, e, f), plex(x, y)):
print("********F4*******");
DISPG([a*x^3*y+c*x*y^2, x^4*y+3*d*y, c*x^2+b*x*y, x^2*y^2+a*x^2, x^5+y^5], plex(a, b, c, d), plex(x, y)):
print("********B6*******");
DISPG([a*x^2*y^2+b*x*z^2, b*x*y^2+c*x^2+2, c*x^2+b*y^2*z], plex(a, b, c), plex(x, y, z)):
print("********B7*******");
DISPG([x^3-a*x*y, x^2*y-2*y^2+x], plex(a), plex(x, y)):
print("********S7*******");
DISPG([y^2-z*x*y+x^2+z-1, x*y+z^2-1, y^2+x^2+z^2-r^2], plex(r), plex(x, y, z)):
print("********ex1*******");
DISPG([a*x^2*y+b*x*y, b*x+y^2], plex(a, b), plex(x, y)):
print("********ex2*******");
DISPG([2599*e2^2-400*e2*e3+80*e2*f3-80*e3*f2+2599*f2^2-400*f2*f3-2200*e2-240*f2-20*Q2, 14-12*e2-110*f2-2*e3-10*f3-P1, 2397-2200*e2+240*f2-200*e3+40*f3-20*Q1, 16*e2^2-4*e2*e3-20*e2*f3+20*e3*f2+16*f2^2-4*f2*f3-12*e2+110*f2-P2], plex(P1, Q1, P2, Q2), plex(e2, f2, e3, f3)):
print("********ex3*******");
DISPG([x^2+b*y^2+2*c*x*y+2*d*x+2*e*y+f, x+c*y+d, b*y+c*x+e], plex(b, c, d, e, f), plex(x, y)):
print("********ex4*******");
DISPG([r-c1+l*(s1*s2-c1*c2), z-s1-l*(s1*c2+s2*c1), s1^2+c1^2-1, s2^2+c2^2-1], plex(r, z, l), plex(c1, s1, c2, s2)):
print("********ex6*******");
DISPG([a*x*z+b*x*z+a, (a^2+a)*x*y], plex(a, b), plex(x, y, z)):
print("********ex7*******");
DISPG([a+b^2, y+1, b*x+1, a*x-b], plex(a, b), plex(x, y)):
print("********ex8*******");
DISPG([x+c*y+b*z+a, c*x+y+a*z+b, b*x+a*y+z+c], plex(a, b, c), plex(x, y, z)):
print("********F1*******");
DISPG([a*x^4*y+x*y^2+b*x, x^3+2*x*y, b*x^2+x^2*y], plex(a, b), plex(x, y)):
print("********F2*******");
DISPG([a*x^2*y^3+b*y+y, x^2*y^2+x*y+2, a*x^2+b*y+2], plex(a, b), plex(x, y)):
print("********F3*******");
DISPG([b*x^3+x^2+2, c*x^2+d*x, a*x^4+c*x^2+b], plex(a, b, c, d), plex(x)):
print("********F4*******");
DISPG([a*x^3*y+c*x*y^2, x^4*y+3*d*y, c*x^2+b*x*y, x^2*y^2+a*x^2, x^5+y^5], plex(a, b, c, d), plex(x, y)):
print("********F5*******");
#DISPG([a*x^2*y+b*x+y^3, a*x^2*y+b*x*y, y^2+b*x^2*y+c*x*y], plex(a, b, c), plex(x, y)):
print("********F6*******");
DISPG([a*x^2*y^2+b*x*z^2, b*x*y^2+c*x^2+2, c*x^2+b*y^2*z], plex(a, b, c), plex(x, y, z)):
print("********F7*******");
#DISPG([x^4+a*x^3+b*x^2+c*x+d, 4*x^3+3*a*x^2+2*b*x+c], plex(a, b, c, d), plex(x)):
print("********S1*******");                         
DISPG([a*(x+y), b*(x+y), x^2+a*x], plex(a, b), plex(x, y)):
print("********S2*******");
DISPG([x1^2, x1*x2, x1*x3^2, x1*a+x2, x2*x3-x3^2, x2*a, x3^3, x3^2*a, a^2], plex(a), plex(x1, x2, x3)):
print("********S3*******");
DISPG([x^3-a*x*y, x^2*y-2*y^2+x], plex(a), plex(x, y)):
print("********S4*******");
DISPG([a*x+y-1, b*x+y-2, 2*x+a*y, b*x+a*y+1], plex(a, b), plex(x, y)):
print("********S5*******");
DISPG([x4-a4+a2, x1+x2+x3+x4-a1-a3-a4, x1*x3*x4-a1*a3*a4, x1*x3+x1*x4+x2*x3+x3*x4-a1*a4-a1*a3-a3*a4], plex(a1, a2, a3, a4), plex(x1, x2, x3, x4)):
print("********S6*******");
DISPG([v*x*y+u*x^2+x, u*y^2+x^2], plex(u, v), plex(x, y)):
print("********S7*******");
DISPG([y^2-z*x*y+x^2+z-1, x*y+z^2-1, y^2+x^2+z^2-r^2], plex(r), plex(x, y, z)):
print("********S8*******");
DISPG([a-b+(x*y*a-x^2*y*b-3*a)^3+(x*y*b-3*x*b-5*b)^4, x*y*a-x^2*y*b-3*a, x*y*b-3*x*b-5*b], plex(a, b), plex(x, y)):
print("********S9*******");
DISPG([x+c*y+b*z+a, c*x+y+a*z+b, x*b+a*y+z+c], plex(a, b, c), plex(x, y, z)):
print("********S10*******");
#DISPG([a-l3*c3-l2*c1, b-l3*s3-l2*s1, c1^2+s1^2-1, c3^2+s3^2-1], plex(a, b, c), plex(c1, c3, l1, l2, l3, s1, s3)):
print("********S15*******");
DISPG([a+d*s1, b-d*c1, l2*c2+l3*c3-d, l2*s2+l3*s3-c, c1^2+s1^2-1, c3^2+s3^2-1], plex(a, b, c, d), plex(c1, c2, c3, l1, l2, l3, s1, s2, s3)):
print("********S16*******");
DISPG([w12+w14, w12+w13, w12+w15, w12+w23+w25-w26*x+w26, w12+w25-w26*y+w26, w12+w23-w26*z+w26, w23+w34+x*w36, w13+w34-w36*y+w36, w23+z*w36, w14+w34+w45-w46*x+w46, w34+y*w46, w45+z*w56, w15+w45-z*w56+w56, -w26+w26*x+x*w36-w46+w46*x+w56*x, -w26+w26*y-w36+w36*y+y*w46+w56*y, -w26+w26*z+z*w36+w46*z-w56+z*w56], plex(x, y, z), plex(w12, w13, w14, w15, w23, w25, w34, w45, w26, w36, w46, w56)):
print("********S13*******");
DISPG([a*x^2*y+a+3*b^2, a*(b-c)*x*y+a*b*x+5*c], plex(a, b, c), plex(x, y)):                                      
print("********B1*******");
#DISPG([a*x^4+b*x^2*y^2-2*c*x^2*y-d*x*y^2+e*y^2+f, 4*a*x^3+2*b*x*y^2-4*c*x*y-d*y^2, 2*b*x^2*y-2*c*x^2-2*d*x*y+2*e*y], plex(a, b, c, d, e, f), plex(x, y)):