##############################################
########			       ####### 
########  Parametric Involutive Bases  #######
########                               #######
##############################################

######Needed packages
with(Groebner):
with(Algebraic):
with(PolynomialIdeals):
with(ArrayTools):
######Select Algorithm
sff := proc (x) 
evalb(type(x, 'constant')): 
end proc:

######Computing the leading term of a poly.
lt := proc (f, T) 
RETURN(LeadingTerm(f, T)[1]*LeadingTerm(f, T)[2]): 
end proc: 

######Normalizing procedure
nrm := proc(F) 
local h, i; 
#option trace; 
if F = [] then 
   RETURN(F); 
fi:
if type(denom(F[1]), monomial) then 
   h := LeadingMonomial(denom(F[1]), plex(op(indets(denom(F[1]))))); 
else 
     h := denom(F[1]); 
fi: 
for i from 2 to nops(F) do 
    if type(denom(F[i]), monomial) then 
        h := lcm(h, LeadingMonomial(denom(F[i]), plex(op(indets(denom(F[i])))))); 
    else 
    	h := lcm(h, denom(F[i])); 
    fi: 
od: 
RETURN(simplify(expand(h*F))): 
end:

nrm:=proc(F)
local i,f,A;
#option tarce;
f:=proc(u)
if type(u, monomial) then
    LeadingMonomial(u, plex(op(indets(denom(u)))));
else
    u;
fi:
end:
A:=map(f,F):
RETURN(simplify(expand(lcm(seq(denom(A[i]),i=1..nops(A)))*F))): 
end:


#####Convert a list of polynomial to the set "FACVAR"
fac := proc(L) 
local AA, AAA, N, P, p, B, C, i; 
#option trace; 
AA := [seq(factor(i), i in L)]; 
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]: 
    fi: 
od: 
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]: 
    fi: 
od: 
RETURN({N}): 
end:


fac2:=proc(L)
local A,i,f,g;
g:=proc(v)
if irreduc(v) then
     v
else
     op(convert(factor(v), list));;
fi:
end:
A := {op(map(g,L))}; 
f:=proc(u)
	not type(u, 'constant');
end:
select(f,A);
end:

###### New cond algorithm
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))))) 
    od:
    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:
    fi:
od: 
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] 
        fi: 
    od: 
    WW := [ww] 
fi: 
fff := ff; 
while type(LeadingCoefficient(fff, T), 'constant') = true and fff <> 0 do 
    fff := simplify(expand(fff-simplify(expand(lt(fff, T))))) 
od; 
cc := fac([LeadingCoefficient(fff, T)]); 
ccc := {seq(op(nrm([k/LeadingCoefficient(k, R)])), k in 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), jj in  www)] 
fi:
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:


#####Canspec algorithm
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:


#####Janet Division
Janet:=proc(u,U,Vars)
local nn,m,d,L,i,j,dd,V,v,Mult;
#option trace;
if degree(u)=0 then RETURN(Vars); fi;
nn:=nops(Vars);
m:=nops(U);
d:=[seq(max(seq(degree(U[j],Vars[i]),j=1..m)),i=1..nn)];
Mult:=NULL;
    if degree(u,Vars[1])=d[1] then
       Mult:=Mult,Vars[1];
    fi:
for j from 2 to nn do
    	dd:=[seq(degree(u,Vars[l]),l=1..j-1)];
    	V:=NULL:
    	for v in U do
    	    if [seq(degree(v,Vars[l]),l=1..j-1)]=dd then
	       V:=V,v:
	    fi:
	od:
	if degree(u,Vars[j])=max(seq(degree(v,Vars[j]), v in [V])) then
	   Mult:=Mult,Vars[j];
	fi:
od:
RETURN([Mult]);
end:

############################
#         Reduction        # 
############################
###### This proc checks involutive divisibility ###### 
InvDivide2:=proc(u,v,L)
global tord,Vars;
local InvDiv,uu,M;
#option trace;
if v=0 then RETURN(false); fi;
if divide(u,v) then
   uu:=u/v;
   M:=Janet(v,L,Vars);
   if indets(uu) subset {op(M)} then
      RETURN(true);
   else
      RETURN(false);
   fi;
else
   RETURN(false);
fi;
RETURN(false);
end:
################ AutoReduced Algorithm #################
AutoReduction:=proc(f,F,TermOrder)
#option trace;
local LMF,Vars,p,r,flag,g,u,Mult,FF,i;
global tord;
LMF:=[seq(LeadingMonomial(v,tord), v in F)];
Vars:=[op(TermOrder)];
p:=f:
r:=0;
FF:=F:#{op(F)} minus {0};
while p<>0 do
      flag:=false:
      i:=1;
      while not flag and i<= nops(FF) do
          g:=FF[i]:
      	  if divide(LeadingMonomial(p,TermOrder),LeadingMonomial(g,TermOrder)) then
	     u:=LeadingMonomial(p,TermOrder)/LeadingMonomial(g,TermOrder);
	     Mult:=Janet(LeadingMonomial(g,TermOrder),LMF,Vars);
	     if indets(u) subset {op(Mult)} then
	     	flag:=true:
		p:=simplify(p-(LeadingCoefficient(p,TermOrder)/LeadingCoefficient(g,TermOrder))*u*g);
	     else 
                i:=i+1;
	     fi:
	     else 
	        i:=i+1;
	  fi:
       od:
       if not flag  then 
	  r:=simplify(r+LeadingCoefficient(p,TermOrder)*LeadingMonomial(p,TermOrder));
          p:=simplify(p-LeadingCoefficient(p,TermOrder)*LeadingMonomial(p,TermOrder));
       fi:
od:
r:=op(nrm([r])):
RETURN(r);
end:


#########################
tidy:=proc(p,q)
global Polys,tord;
if LM(Polys[p[1]])<>LM(Polys[q[1]]) then
   RETURN(TestOrder(LM(Polys[p[1]]),LM(Polys[q[1]]),tord));
else
   RETURN(evalb(p[2]<q[2]));
fi;
end:
#########################
LM:=proc(f)
global tord;
if f<>0 then
    RETURN(LeadingMonomial(f,tord));
fi:
RETURN(0);
end:
#################
LC:=proc(f)
global tord;
RETURN(LeadingCoefficient(f,tord));
end:    
#####################
Criteria2:=proc(p,gg)
global Tpol,tord,Vars,Q,Polys,c1,c2,LL,ALM;
local InvDiv,flag,t,i,NM,y,g;
#option trace;
g:=LM(Polys[gg[2]]);
    if LM(Polys[p[2]])*g=LM(Polys[p[1]]) then 
    c1:=c1+1;
    RETURN(true);
elif  divide(LM(Polys[p[1]]),lcm(LM(Polys[p[2]]),g)) and LM(Polys[p[1]])<>lcm(LM(Polys[p[2]]),g) then 
    c2:=c2+1;
    RETURN(true);
fi;
RETURN(false);
end:
#####################
#################
LM:=proc(f)
global tord;
if f<>0 then
    RETURN(LeadingMonomial(f,tord));
fi:
RETURN(0);
end:
#################
LC:=proc(f)
global tord;
RETURN(LeadingCoefficient(f,tord));
end:    
#####################
#####################
HeadNormalForm:=proc(p,BB,NN,WW,Rr, Tt ,CPP,RRRR)
global Tpol,tord,Q,Polys,Vars;
local h,G,flag,i,divisor,g,cd,NNN,WWW,Y;
#option trace;
Tpol:=BB: 
h:=Polys[p[1]];
G:=[seq(LM(Polys[q[1]]), q in Tpol)];
flag:=false:
for i from 1 to nops(Tpol) while not flag do
    if InvDivide2(LM(h),LM(Polys[Tpol[i][1]]),G) then 
	divisor:=i;
	flag:=true;
    fi:
od:
if not flag then
   RETURN(Polys[p[1]]);
fi;
if p[1]<>p[2] then 
     if Criteria2(p,Tpol[divisor]) then
        RETURN(c1c2c3c4);
     else    
       i:=1;
       while h<>0 and i<=nops(Tpol) do 
   	 g:=Polys(Tpol[i][1]);
	 if InvDivide2(LM(h),LM(g),G) then
	     h:=SPolynomial(h,g,tord);
             h:=NormalForm(h,NN,Rr):
	     if h<>0 then
	     	Y := new(NN, WW, h, Rr, Tt);
		cd := [op(Y[1])];
		h := simplify(expand(Y[2]));
		NNN := Y[3]; 
		WWW := Y[4]; 
	     	if cd <> [] then
                         RETURN(false,h);
                fi	     	
	     fi:
             i:=1;
         else
	     i:=i+1;
         fi;
       od;
    fi:   
else
   i:=1;
   while h<>0 and i<=nops(Tpol) do 
   	 g:=Polys(Tpol[i][1]);
	 if InvDivide2(LM(h),LM(g),G) then
	     h:=SPolynomial(h,g,tord);
             h:=NormalForm(h,NN,Rr):
	     if h<>0 then
	     	Y := new(NN, WW, h, Rr, Tt);
		cd := [op(Y[1])];
		h := simplify(expand(Y[2]));
		NNN := Y[3]; 
		WWW := Y[4]; 
	     	if cd <> [] then
                          RETURN(false,h);
                fi:
              fi:	     	
             i:=1;
         else
	     i:=i+1;
         fi;
    od;
fi;
RETURN(h);
end:
#####################
HeadReduce:=proc(BB,NN,WW,Rr, Tt ,CPP,RRRR)
#option trace;
global Vars,Tpol,tord,Polys,rdz,isr,n;
local S,Q,p,h,SS,s,RR;
S:=CPP:
Q:=NULL;
while nops(S)<>0 do
     p:=S[1];
     S:=S[2..-1];
     h:=HeadNormalForm(p,BB,NN,WW,Rr, Tt ,CPP,RRRR); 
     if nops([h])>1 then 
          if h[1]=false then 
               RR:=RRRR: 
               RR[p[1]]:=h[2]: 
	       RETURN(false,p,BB,NN,WW,CPP,RR); 
          fi:
     fi:
     if h<>0 and h<>c1c2c3c4  then
     	if LM(Polys[p[1]])<>LM(h) then
	     n:=n+1:
	     Polys(n):=h;
             Q:=Q,[n,n,{}];
	else
             Polys(p[1]):=h;
	     Q:=Q,p;
        fi;
     elif h=0  then 
        rdz:=rdz+1:
        if p[1]=p[2] then
	   SS:=NULL;
           for s in S do
	       if s[2]<>p[1] then
	       	  SS:=SS,s;
	       else
	          isr:=isr+1:
               fi;
           od:
           S:=[SS];
         fi:
      fi:
od:
RETURN([Q]);
end:                 

#####################
Gerdt:=proc(ffff,BB, NN, WW, Rr, Tt ,CPP,RRRR)
#option trace;
#global tord,Tpol,Q,Polys,Vars,c1,c2,rdz,ALM,isr,n;
global Polys,Tpol,Vars,tord,c1,c2,rdz,n,deg,isr;
local InvDiv,F,p,TT,q,h,NM,np,j,firsttime,firstbytes,secondtime,secondbytes,test,Q,yy,f,i;  
tord:=Tt;  
Vars:=[op(Tt)];
test := true;
F:=BB;
Q:=CPP:
if Q=[] then #and BB=[] then 
    c1:=0:
    c2:=0:
    rdz:=0;
    deg:=0:
    isr:=0:
    f:=proc(u)
      u<>0:
    end:
    F:=convert(select(f,convert(RRRR,list)),Array):
    n:=NumElems(F):if n=0 then RETURN([true,0, [[1]], NN, WW,[],[0]]);fi;
    F:=sort(F,(a,b)->TestOrder(a,b,tord));
    Polys:=F:
    Tpol:=[[1,1,{}]];
    Q:=[seq([i,i,{}],i=2..n)];
else 
    Tpol:=BB: 
    Polys:=RRRR:
    Q:=CPP:
fi:
Q:=sort(Q,tidy);
Q:=HeadReduce(Tpol,NN,WW,Rr, Tt ,Q,Polys); 
if Q=[] then  RETURN(true,0,Tpol,NN,WW,[],Polys); fi:
if nops([Q])>1 then if Q[1]=false  then RETURN(Q); fi: fi: 
while  nops(Q)<>0 do
       p:=Q[1]; 
       Q:=Q[2..-1]; 
       if p[1]=p[2]   then
       	  TT:=NULL;
          for q in Tpol do
	           if divide(LM(Polys[q[1]]), LM(Polys[p[1]])) and LM(Polys[p[1]]) <> LM(Polys[q[1]]) then
	      	        Q:=[op(Q),q]; 
                else
		        TT:=TT,q;
                fi:
          od:
          Tpol:=[TT];
        fi;
	#deg:=max(deg,degree(Polys[p[1]],{op(Vars)})):
        h:=AutoReduction(Polys[p[1]],[seq(Polys[q[1]],q in Tpol)],tord); 
	Polys[p[1]]:=h;
        if h<>0 then
             Tpol:=[op(Tpol),p]; 
             for i from 1 to nops(Tpol) do 
       	         q:=Tpol[i];
       	         NM:={op(Vars)} minus {op(Janet(LM(Polys[q[1]]),[seq(LM(Polys[w[1]]), w in Tpol)],Vars))};
	             for yy in NM minus q[3] do
		       n:=n+1: 
	               Polys(n):=expand(yy*Polys[q[1]]);
	               Q:=[op(Q),[n,q[2],{}]];
	   	     od:
                 Tpol[i][3]:=(NM minus q[3]) union (NM intersect q[3]);
             od:
             Q:=sort(Q,tidy);
             Q:=HeadReduce(Tpol,NN,WW,Rr, Tt ,Q,Polys);
             if Q=[] then  RETURN(true,0,Tpol,NN,WW,[],Polys); fi:
             if nops([Q])>1 then if Q[1]=false  then RETURN(Q); fi: fi: 
             
        fi:
od:
#Tpol:=[seq(Polys[tt[1]], tt in Tpol)];
RETURN(true,0,Tpol,NN,WW,[],Polys);
end:	   

#####Branch algorithm
branch := proc (ffff,BBB, NNN, WWW, R, T,CPP,RR) 
local cd,l,ll,i,ff,Y,X,BB,pivot,test,TTT,Bb,g,TT,LL,gg,B,CN,N,W,NN,WW,f,CP,RRRR;
global LIST, Grob,termord1, termord2, termord3, iii, DDD,ind,rdz,inputsize,RR1; 
#option trace:
RRRR:=RR:
#CN := canspec(NNN, WWW, R);
B:=BBB:
N:=NNN:#CN[2];
W:=WWW:#CN[3];
cd := []; 
Bb:=B: 
f:=RRRR[ffff[1]]:
Y := new(N, W, f, R, T);
cd := [op(Y[1])];
ff := simplify(expand(Y[2]));
NN := Y[3]; 
WW := Y[4]; 
#RRRR[ffff[1]]:=ff:
if ind<inputsize and cd <> [] then 
    RR1:=RRRR:
    #RR1[ffff[1]]:=NormalForm(ff, Basis([op(NN),op(cd)],R),R):
    RETURN(branch(ffff,Bb, NN, [op(WW),op(cd)], R, T,CPP,RRRR),branch(ffff,Bb, [op(NN),op(cd)], WW, R, T,CPP,RR1)):
elif ind<inputsize and cd = [] and ff<>0  then 
    RETURN([ffff,Bb, NN, WW, R, T,CPP,convert(NormalForm(convert(RRRR,list),NN,R),Array)]):
elif ind<inputsize and cd = [] and ff=0  then 
    RETURN([[],Bb, NN, WW, R, T,CPP,convert(NormalForm(convert(RRRR,list),NN,R),Array)]):
fi:
if cd = [] then 
    X:=Gerdt(ffff,Bb, NN, WW, R, T,CPP,RRRR);
    test := X[1]; 
    ff := X[2]; 
    BB := X[3]; 
    NN := X[4]; 
    WW:=X[5];
    CP:=X[6];
    RRRR:=X[7]:
    if test  then 
	if CP=[] then 
	    LIST := LIST, [seq(RRRR[tt[1]], tt in BB), NN, WW]:
        
	fi;
    else 
        branch(ff, BB, NN, WW, R, T,CP,RRRR); 
    end if: 
else
    branch(ffff,Bb, NN, [op(WW),op(cd)], R, T,CPP,RR);  
    branch(ffff,Bb, [op(NN),op(cd)], WW, R, T,CPP,convert(NormalForm(convert(RRRR,list),[op(NN),op(cd)],R),Array)); 
   
fi:
end:

################## Gerdt-Montes algorithm
Montes:=proc(FF,R,T)
local f,t1,t2,b1,b2,F;
global LIST,DDD,ind,rdz,Grob,termord1, termord2,inputsize;
local A,G,i;
#option trace;
t1,b1:=kernelopts(cputime,bytesused);
termord1, termord2:=R,T:
inputsize:=nops(FF);
F:=Array(1..inputsize,i->FF[i]);
rdz:=0;
DDD:=0;
LIST:=NULL;
ind:=1:
LIST:=NULL:
G:=[branch([1,1,{}],[], [], [], R, T,[],F)];
for i from 2 to inputsize do 
    ind:=ind+1;
    G:=[seq(branch([i,i,{}],op(A[2..8])), A in G)];
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):
#lprint(LIST);
RETURN(LIST);
end:

#print("********B2*******");
F:=[a*x^2*y-y^3,b*x+y^2+c]:
Montes(F, plex(a, b,c), plex(x, y)); 
print("********B3*******");
F:=[a*y^2+b*x^2+c*x^3+d, 2*a*y, 3*b*x+3*c*x^2]:
#Montes(F, plex(a, b, c,d), plex(x, y));    
#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("********B4*******");
F:=[x+c*y+b*z+a, c*x+y+a*z+b, b*x+a*y+z+c];
#Montes(F, plex(a, b, c), plex(x, y, z));
#Montes([x^3-a*x*y, x^2*y-2*y^2+x], plex(a), plex(x, y));
F:=[x+a*y^2+b*y+c,d*y^2+e*y+f]:

#Montes(F, plex(a, b, c,d,e,f), tdeg(x, y));
#Montes(F, plex(a, b, c,d,e,f), tdeg(y,x));