restart:
with(PolynomialIdeals):
with(DifferentialAlgebra):with(Tools):

R := DifferentialRing( derivations = [x,y,z], blocks = [w,v,u,m] ):
vars:=[x,y,z]:
parms:=[a,b,c,d]:
F:=DifferentialRing(derivations=vars,blocks=[w,v,u,m],arbitrary=parms):
F1:=DifferentialRing(derivations=vars,blocks=[[m,u,v,w]],arbitrary=parms):
#..............................................................................................
#gcd of memebers of L
GCD:=proc(L)
local g,i:
#option trace:
g:=L[1]:
if nops(L)>1 then
for i from 1 to nops(L) do
g:=gcd(g,L[i]):
od:
fi:
return(g):
end:
#...........................................................................................................
#it divides f by u^m if u belongs to U 
DIVIDE:=proc(f,U)
local i,q,sw,g:
if f=0 then return(0):else
q:=f:
for i to nops(U) do
sw:=true:
while sw do
g:=gcd(q,U[i]):
if g<>1 then q:=expand(simplify(q/g)): else sw:=false:fi:
end do:
od:
return(q):
fi:
end:
#....................................................................................................
#it removes {0} from L
RZ:=proc(L)
local i,S:
S:=[]:
for i to nops(L) do
if L[i]<>0 then S:=[op(S),L[i]]:fi:
end do:
return(S):
end:
RZ([]):

#..................................................................................................
#it simplifies the coeffs of f
SIM2:=proc(f)
local v:
v:=indets(f) minus {op(parms)}:
if f=0 then return(f):
elif type(f,'monomial') then return(f/DIVIDE(f,[op(v)])):
else return(expand((f/DIVIDE(GCD([op(expand(f))]),[op(v)])))):fi:
end:
#..........................................................................................................
SIM:=proc(f)
if f=0 then return(f):else return(f/GCD([coeffs(expand(f))])):fi:
end:

#......................................................................................................
#it adds the members of M to L without repeated members
ADD:=proc(L::list,M::list)
local s,i:
if L=[] and M=[] then return([]):
else
if L=[] then s:=[M[1]]:else
s:=L:fi:
for i to nops(M) do
if {member(M[i],s),member(-M[i],s)}={false} then s:=[op(s),M[i]]:fi:
end:
return(s):
fi:
end:
#..................................................................................................
#it extracts the members of T from L
Extract:=proc(L::list,T::list)
local S,i:
S:=[]:
if T=[] then return(L):
else
for i to nops(L) do
if {seq(is(L[i]=T[j]),j=1..nops(T))}={false} then S:=[op(S),L[i]]:fi:
end:
return(S):
fi:
end:

#...................................................................................................................
ParmDiffRed:=proc(p,L,r,R)
local H:
#option trace:
H:=[DifferentialPrem(p,[op(L),seq(seq(parms[j][vars[i]],i=1...nops(vars)),j=1...nops(parms))],r,R)]:
return(H[2]):
end:
#.............................................................................................................
#if p is reduceble w.r.t q then it returns true else returns false
isreduceble:=proc(p,q,R)
local dp,dq,dpp,dqq,H,i,sw:
#option trace:
H:=indets(p):
dq:=LeadingDerivative(q,R):
dqq:=[FactorDerivative(dq,R)]: degree(q,dq):
i:=1:
sw:=false:
while i<=nops(H) and sw=false do 
dp:=H[i]:
dpp:=[FactorDerivative(dp,R)]:
 if dp=dq and degree(q,dq)<=degree(p,dp) then sw:=true:
elif dp=dq and degree(q,dq)>degree(p,dp) then i:=i+1:
elif dpp[2]=dqq[2] and divide(dpp[1],dqq[1]) then sw:=true:
else i:=i+1:
fi:
end do:
if sw=true then return(true): else return(false):fi:
end:
#.................................................................................................................
#if the leader of p is reduceble w.r.t leader q then it returns true else returns false
leaderreduced:=proc(p,q,R)
local dp,dq,dpp,dqq:
#option trace:
dp:=LeadingDerivative(p,R):
dq:=LeadingDerivative(q,R):
dqq:=[FactorDerivative(dq,R)]:
dpp:=[FactorDerivative(dp,R)]:
if dp=dq and degree(q,dq)<=degree(p,dp) then return(true):
elif dp=dq and degree(q,dq)>degree(p,dp) then return(false):
elif dpp[2]=dqq[2] and divide(dpp[1],dqq[1]) then return(true):
else return(false):
fi:
end:
#first criterion of buchberger.......................................................................................

Test:=proc(f,R)
local H:
#option trace:
if nops(indets(f))<>1 then return(false) else
if type(f,`+`) then H:=[op(f)]: else H:=[f]:fi:
if {seq(type(H[i],'constant'),i=1..nops(H))}<>{false} then return(false) else
 if {seq(type(Initial(H[i],R),'constant'),i=1..nops(H))}<>{true} then return(false) else
   if nops({seq(degree(FactorDerivative(SIM(H[i]), R)[1]),i=1..nops(H))})<>1 then return(false): else
     return(true):fi:
   fi:
 fi:
fi:
end:
 firstcri:=proc(p,a,R)
 #option trace:
 if Test(p,R) and Test(a,R) and gcd([FactorDerivative(LeadingDerivative(p,R),R)][1],[FactorDerivative(LeadingDerivative(a,R),R)][1])=1 then
   return(true):
 else
   return(false):
 fi:
 end:
#second criterion of buchberger.............................................................................
SC:=proc(P::list,R)
local H,S:
#option trace:
H:=[seq(FactorDerivative(SIM(LeadingDerivative(P[i],R)), R)[1],i=1..3)]:
if nops(convert(H,set))<>3 then return(false): else
  if GCD(H)=1 then return(false): else
   if divide(lcm(H[1],H[3]),H[2])=false then return(false): else 
     if {seq(seq(divide(H[i],H[j]),i=1..3),j in Extract([1,2,3],[i]))}={false} then return(true): else
       S:=SortByRank(P,ascending,R):
       if S[2]=P[2] then return(true):
       elif degree(S[2],LeadingDerivative(S[2],R))=1 then return(true):
       else return(false):
       fi:
     fi:
   fi:
 fi:
fi:
end: 
SC([u[x],v[y], b*u[x, x]*v],F):
Sec_Crit:=proc(f,g,ZD,R)
local i,sw,n,j:
#option trace:
if ZD=[] then return(false) else
sw:=true:
i:=1:
while sw and i<=nops(ZD) do
 n:=num(f,ZD[i]):
 if n=false then i:=i+1:
 else
  if n=1 then j:=2:else j:=1:fi:
  if member({g,ZD[i][j]},ZD) then
   if SC([f,ZD[i][j],g],R)=true then  sw:=false: else i:=i+1:fi:
  else i:=i+1:
  fi:
 fi:
end do:
if sw=true then return(false) else return(true):fi:
fi:
end:


# update subalgorithm..............................................................................................
update:=proc(h,Gold,Bold,R)
local i,C,D,E,d,T,s,hd,Bnew,Gnew,B,k,Ex,f,K,DD,ldk1,ldk2:
#option trace:
#C:={seq({f,Gold[i]},i=1..nops(Gold))}:
D:=[]:
f:=SIM(h):
E:=[]:C:={}:K:=[]:
d:=[FactorDerivative(SIM(LeadingDerivative(f,R)),R)]:
i:=1:Gnew:=Gold:Ex:=[]:
while i<=nops(Gold) do
  if leaderreduced(Gold[i],f,R)=true then Gnew:=Extract(Gnew,[Gold[i]]):K:=[op(K),Gold[i]]:
  elif isreduceble(Gold[i],f,R)=true then Ex:=[op(Ex),Gold[i]]: Gnew:=Extract(Gnew,[Gold[i]]): fi:
  i:=i+1:
end do:
DD:=ADD(Gnew,K):
for i from 1 to nops(DD) do
if [FactorDerivative(SIM(LeadingDerivative(DD[i],R)),R)][2]=d[2] then
C:=C union {{SIM(DD[i]),f}}:
fi:
od:
T:=Test(f,R):
while C<>{} do
  s:=C[1]:C:=C minus {s}:
  hd:=[FactorDerivative(SIM(LeadingDerivative(s[2],R)),R)]:
  if T=true and Test(s[2],R)=true and gcd(hd[1],d[1])=1 then
    D:=[op(D),s]:
  elif Sec_Crit(f,s[2],C,R)=false or Sec_Crit(f,s[2],D,R)=false then
   D:=[op(D),s]:E:=[op(E),s]:
  fi:
end do:
Bnew:=[]:
B:=Bold:i:=1:
while i<=nops(B) do
  k:=B[i]:B:=Extract(B,[k]):
  ldk1:=[FactorDerivative(SIM(LeadingDerivative(k[1],R)),R)][1]:
  ldk2:=[FactorDerivative(SIM(LeadingDerivative(k[2],R)),R)][1]:
  if SC([k[1],f,k[2]],R)=false or gcd(ldk1,ldk2)=gcd(ldk1,d[1]) or gcd(ldk1,ldk2)=gcd(ldk2,d[1]) then Bnew:=[op(Bnew),k]:fi:
end do:
Bnew:=[op(Bnew),op(E)]:
Gnew:=[op(Gnew),f]:
return(Gnew,Ex,Bnew):
end:
# Find-New-Condition sub algorithm.................................................................................
#f is a new polynomial and OC is the list of old conditions,NZP is the non zero parametrs.
NewCondition:=proc(f,oc,p,NZP,R)
local init,inn,s,sep,NC,E,F,e,OC,OC1:
#option trace:
OC1:=ADD(oc,p):
OC:=ADD(OC1,NZP):
E:=[]:F:=p:
inn:=SIM(Initial(f,R)):
s:=SIM(Separant(f,R)):
init:=DIVIDE(inn,OC):
if type(init,'constant')=false then
 if  member(0,[seq(ParmDiffRed(OC1[i],[init],reduction=full,R),i=1..nops(OC1))])=false then 
E:=[init]:OC1:=ADD(OC1,[init]): e:=i:else F:=[op(F),init]:fi:
fi:
sep:=DIVIDE(expand(s),OC):
if type(sep,'constant')=false then
 if member(0,[seq(ParmDiffRed(OC1[i],[sep],reduction=full,R),i=1..nops(OC1))])=false then E:=ADD(E,[sep]):e:=se 
elif init<>sep then F:=[op(F),sep]:print(yes):fi:
fi:
 if E=[] then return([false,F]):
 else return([true,E,F,e]):
 fi:
end:
#Make-Tree subalgorithm ..........................................................................................
#P=[processed-poly,poly,delta,non-processed-polys,nzc] is a vertex,f is a polynomial with new conditions and C is  its #condition
Branch:=proc(P,p,R)
local NewP,f1,g,delta,T,TR,Re,NC,Z,i,P4,P42,Z2,Z3,New,O,h,count,D2,B,f,o:
#option trace:
B:=P[3]:
f:=factor(p):
if type(f,`^`) then f:=[op(f)][1]:fi:
if type(f,`+`)=false then
h:=SIM(f):
 if type(h,`*`)=true then O:=[op(h)]:New:=[]:
  for i from 1 to nops(O) do
   if type(O[i],'constant')=false then
   if type(O[i],`^`) then o:=[op(O[i])][1]:else o:=O[i]:fi:
    New:=[op(New),[P[1],ADD([o],P[2]),P[3],P[4],P[5],P[6],P[7]]]:
   fi:
  end do:
  NewP:=op(New):
 else
   Re:=update(f,P[1],B,R):
   TR:=Re[1]:
    B:=Re[3]:
   NewP:=[RZ(TR),ADD(Re[2],Extract(P[2],[f,-f])),B,P[4],P[5],P[6],P[7]]:
 fi:
else
NC:=NewCondition(f,P[4],P[5],P[7],R):Z:=[]:Z2:=[]:P42:=P[4]:Z3:=[]:
 if NC[1]=false then 
    Z:=NC[2]:
    Re:=update(f,P[1],B,R):
    TR:=[seq(DIVIDE(Re[1][i],NC[2]),i=1..nops(Re[1]))]:
    B:=Re[3]:
   NewP:=[RZ(TR),ADD(Re[2],Extract(P[2],[f,-f])),B,P[4],Z,P[6],P[7]]:
 elif nops(NC[2])=1  then
 P4:=P[4]:
 #for i from 1 to nops(P[4]) do
  #if ParmDiffRed(NC[2][1],[P[4][i]],reduction=full,R)=0 then P4:=Extract(P4,[P[4][i]]):Z3:=ADD(Z3,[P[4][i]]):    fi:
 #od:
   Re:=update(f,P[1],B,R):
   TR:=[seq(DIVIDE(Re[1][i],NC[2]),i=1..nops(Re[1]))]:
   Z:=ADD(Z,NC[3]):Z3:=ADD(Z3,NC[3]):
   B:=Re[3]:
  if NC[4]=i then
     f1:=expand(f-LeadingDerivative(f,R)^degree(f,LeadingDerivative(f,R))*Initial(f,R)):
  else 
     f1:=expand(degree(f,LeadingDerivative(f,R))*f-LeadingDerivative(f,R)*Separant(f,R)):
  fi:
   NewP:=[TR,ADD(Re[2],Extract(P[2],[f,-f])),B,ADD(P4,[NC[2][1]]),Z3,P[6],P[7]],[P[1],RZ(ADD(Extract(P[2],[f,-f]),[NC[2][1],SIM(f1)])),Extract(P[3],[f,-f]),P[4],Z,P[6],P[7]]:
 else P4:=P[4]:
  for i from 1 to nops(P[4]) do
  if ParmDiffRed(NC[2][1],[P[4][i]],reduction=full,R)=0 then P4:=Extract(P4,[P[4][i]]):Z:=ADD(Z,[P[4][i]]):      P42:=Extract(P42,[P[4][i]]):Z2:=ADD(Z,[P[4][i]]):
  elif  ParmDiffRed(NC[2][2],[P[4][i]],reduction=full,R)=0 then P4:=Extract(P4,[P[4][i]]):Z:=ADD(Z,[P[4][i]]):
  fi:
 od:
 f1:=expand(f-LeadingDerivative(f,R)^degree(f,LeadingDerivative(f,R))*Initial(f,R)):
  g:=expand(degree(f,LeadingDerivative(f,R))*f-LeadingDerivative(f,R)*Separant(f,R)): 
   Re:=update(f,P[1],B,R):
   TR:=seq(DIVIDE(Re[1],NC[2]),i=1..nops(Re[1])):
   Z:=ADD(Z,NC[3]):Z2:=ADD(Z2,NC[3]):
   B:=Re[3]:
  NewP:=[RZ(TR),ADD(Re[2],Extract(P[2],[f,-f])),B,ADD(P4,NC[2]),Z,P[6],P[7]],[P[1],RZ(ADD(Extract(P[2],[f,-f]),[ NC[2][1],SIM(f1)])),Extract(P[3],[f,-f]),P[4],Z,P[6],P[7]],[P[1],RZ(ADD(Extract(P[2],[f,-f]),[NC[2][2],SIM(g)])),Extract(P[3],[f,-f]),ADD(P42,[NC[2][1]]),Z2,P[6],P[7]]:
   
 fi:
fi:
return(NewP):
end:
#.....................................................................................................................
#it selects the components of f which just involve with parametrs.
FACTOR:=proc(f,parm::set)
local v,h,tt,kg:
#option trace:
v:=indets(f):
if v intersect parm=v then
  return(f):
else
  h:=factor(f):
  if type(h,`*`)=false then 
    return(false):
  else
    h:=[op(h)]:kg:=1:
    for tt in h do
     if indets(tt) minus parm={} then kg:=kg*tt:fi:
    od:
    if type(kg,'constant') then return(false): else return(kg):fi:
  fi:
fi:
end:
#....................................................................................................................
#p=polynomial.....N=list zero parametrs...W=list of non zero parametrs
NewPCondition:=proc(p,N,W,R)
local i,par,Ini,In,Se,See,v,v1,w,w1,r:
#option trace:
par:=[]:
if indets(p) minus {op(parms)}={} then par:=p:
else
 Se:=SIM(Separant(p,R)):
 See:=FACTOR(Se,{op(parms)}):
 if See<>false then 
  if N=[] or RadicalMembership(See,<op(N)>)=false then
    See:=DIVIDE(See,W):
    if type(See,'constant')=false then
      r:=s: par:=[See]:
    fi:
  else par:=false:
  fi:
 fi:
if par<>false then
 In:=SIM(Initial(p,R)):
  Ini:=FACTOR(In,{op(parms)}):
  if Ini<>false then 
    Ini:=DIVIDE(Ini,W):
    if type(Ini,'constant')=false then
     r:=i: 
     if type(DIVIDE(Ini,par),'constant')=false then
       par:=ADD([Ini],par):
     fi:
    fi:
  fi:
 fi:
fi:
if par=p then return (false,p):
elif par=[] then return(false,[]):
elif par=false then return(false,[See]):
elif nops(par)=1 then return(true,par,r):
else return(true,par):
fi:
end:
#................................................................................................................
#it outputs the prim components of p
PComponent:=proc(p)
local f,s1,h:
f:=factor(p):
h:=[]:
if type(f,`^`) then return([[op(f)][1]]):
elif type(f,`*`)=false then
  return([p]):
else
  f:=[op(f)]:
  for s1 in f do
   if type(s1,'constant')=false then
    if type(s1,`^`)=false then h:=[op(h),s1]:
    else
      h:=[op(h),[op(s1)][1]]:
    fi:
   fi:
  od:
fi:
return(h):
end:
#....................................................................................................................
MakeTree:=proc(p,P,R)
local NC,P1,P12,q,q1,q2,NewP,PC,New,New1,PC2,i,P2,j:
#option trace:
New:=[]:New1:=[]:P2:=P[2]:P1:=P[1]:
NC:=[NewPCondition(p,P[6],P[7],R)]:
if NC[1]=false and NC[2]=p then
PC:=PComponent(p):
  for i from 1 to nops(PC) do
   for j from 1 to nops(P[1]) do
    if Groebner[NormalForm](P[1][j],[PC[i]],plex(op(parms)))<>P[1][j] then
     P1:=Extract(P1,[P[1][j]]): P2:=ADD([P[1][j]],P2):fi:
   od:
   #P1:=[seq(DIVIDE(Groebner[NormalForm](P[1][j],[PC[i]],plex(op(parms))),ADD(P[4],P[5])),j=1..nops(P[1]))]:
   
   NewP:=[P1,P2,P[3],P[4],P[5],ADD(P[6],[PC[i]]),P[7]]:
 od:
elif NC[1]=false and NC[2]=[] then NewP:=Branch(P,p,R):
elif NC[1]=false then print("see"):
  q:=expand(degree(p,LeadingDerivative(p,R))*p-LeadingDerivative(p,R)*Separant(p,R)):
  NewP:=[P[1],ADD([q],P[2]),P[3],P[4],P[5],P[6],P[7]]:
elif NC[1]=true and nops(NC[2])=1 then 
  if NC[3]=i then 
    q:=expand(p-LeadingDerivative(p,R)^degree(p,LeadingDerivative(p,R))*Initial(p,R)):
  else
    q:=expand(degree(p,LeadingDerivative(p,R))*p-LeadingDerivative(p,R)*Separant(p,R)):
  fi:
  PC:=PComponent(NC[2][1]):
  for i from 1 to nops(PC) do
    for j from 1 to nops(P[1]) do
    if Groebner[NormalForm](P[1][j],[PC[i]],plex(op(parms)))<>P[1][j] then
     P1:=Extract(P1,[P[1][j]]): P2:=ADD([P[1][j]],P2):fi:
   od:
   #P1:=[seq(Groebner[NormalForm](P[1][j],[PC[i]],plex(op(parms))),j=1..nops(P[1]))]:
   New:=[op(New),[P1,ADD([q],P2),P[3],P[4],P[5],ADD(P[6],[PC[i]]),P[7]]]:
 od:
 NewP:=op(New),Branch([P[1],P[2],P[3],P[4],P[5],P[6],ADD(P[7],PC)],DIVIDE(p,PC),R):
else
  q1:=expand(p-LeadingDerivative(p,R)^degree(p,LeadingDerivative(p,R))*Initial(p,R)):
  q2:=expand(degree(p,LeadingDerivative(p,R))*p-LeadingDerivative(p,R)*Separant(p,R)):
  PC:=PComponent(NC[2][1]):
  PC2:=PComponent(NC[2][2]):
  for i from 1 to nops(PC1) do
   for j from 1 to nops(P[1]) do
    if Groebner[NormalForm](P[1][j],[PC[i]],plex(op(parms)))<>P[1][j] then
     P1:=Extract(P1,[P[1][j]]): P2:=ADD([P[1][j]],P2):fi:
   od:
   #P1:=[seq(Groebner[NormalForm](P[1][j],[PC[i]],plex(op(parms))),j=1..nops(P[1]))]:
   New:=[op(New),[P1,ADD([q1],P2),P[3],P[4],P[5],ADD(P[6],[PC[i]]),P[7]]]:
  od:
  for i to nops(PC2) do
  P12:=[seq(Groebner[NormalForm](P[1][j],[PC2[i]],plex(op(parms))),j=1..nops(P[1]))]:
  New1:=[op(New2),[P12,ADD([q2],P[2]),P[3],P[4],P[5],ADD(P[6],[PC2[i]]),ADD(P[7],PC1)]]:
  od:
  NewP:=op(New),op(New1),Branch([P[1],P[2],P[3],P[4],P[5],P[6],ADD(ADD(P[7],PC1),PC2)],DIVIDE(p,[op(PC),op(PC2)]),R):
fi:
return(NewP):
end:
#.....................................................................................................................
#pol is the set of proccesd polynomials and zp,nzp are the set of zero and non zero parms
CheckParm:=proc(Pol,Zp,Nzp)
local N,h,ZP,dd:
ZP:=[]:
for dd in Pol do
  if type(dd,'constant') then return(false):
  elif indets(dd) minus {op(parms)}={} then ZP:=ADD([dd],ZP):fi:
od:
N:=Groebner[Basis](ADD(ZP,Zp),plex(op(parms))):
if N=[1] then return(false): else
  h:=product(Nzp[i],i=1..nops(Nzp)):
  if Nzp=[] or IdealMembership(h,<op(N)>)=false then 
    return(N,Nzp,ZP):
  else
    return(false):
  fi:
fi:
end:

#Differential-consistency subalgorithm............................................................................
CheckBranch:=proc(P,R)
local i,sw,j,L,S,L1,Z,F,C,B,P1,P2,g:
#option trace:
C:=CheckParm(P[1],P[6],P[7]):
if C=false then return(false):
else
P1:=Extract(P[1],C[3]):
P2:=RZ(P[2]):
  S:={seq(type(P2[i],'constant'),i=1..nops(P2))}:
  if member(true,S) then return(false):
  elif P[4]=[] and P[5]=[] then return([P1,P2,P[3],P[4],P[5],C[1],C[2]]):
  else
    Z:=[seq(Groebner[NormalForm](P[4][i],C[1],plex(op(parms))),i=1..nops(P[4]))]:
    L:=[seq(DIVIDE(SIM(expand(ParmDiffRed(Z[i],P1,reduction=full,R))),[seq(Z[j],j=i+1..nops(Z)),op(P[7])]),i=1..nops(P[4]))]:
    if member(0,L) then return(false):
    else L1:=L:
     B:=[]:
     for i from 1 to nops(L) do
       if type(L[i],'constant') then L1:=Extract(L1,[L[i]]):
       elif indets(L[i]) minus {op(parms)}={} then L1:=Extract(L1,[L[i]]): B:=ADD([L[i]],B):
       else
         g:=FACTOR(L[i],{op(parms)}):
         if g<>false then L1:=ADD(Extract(L1,[L[i]]),[expand(L[i]/g)]):B:=ADD(PComponent(g),B):fi:
       fi:
     od:
     #F:=convert({op(P[4])} minus {op(L1)},list):
     return([P1,P2,P[3],L1,P[5],C[1],ADD(C[2],B)]):fi:
  fi:
fi:
end:


#Rosenfeld's lemma.............................................................................
#it detects consistent regular systems by rosenfeld lemma
with(PolynomialIdeals):
checkregular := proc (X,R)
 local L,a1,b1,B,A,m,C,A1:
#option trace:
 L := NULL:
 for a1 in X do
  b1 := [op(a1[1]), op(a1[2]),op(a1[3]),op(a1[4])]:
  B := SortByRank(indets(b1),descending,R):
  m := nops(B):
  A := subs(seq(B[j] = N[j], j = 1 .. m), a1):A1:=ADD(A[2],A[4]):
  C := [op(A[1]),op(A[3]),seq(g[i]*A1[i]-1,i=1..nops(A1))]:
 if not IdealMembership(1, EliminationIdeal(<(op(C))>, {seq(N[i],i=1...m)})) then
   L := L, a1:
 end if:
 end do:
  return(L):
 end proc:
#MainProc.........................................................................................................

MainProc:=proc(P,Q,R)
local Decom,de,q,chBr,T,NP,Br,chBr1,chBr2,Br4:
#option trace:
Decom:=[]:
NP:=[]:
Br:=[[[],P,[],Q,[],[],[]]]:
 while Br<>[] do
  chBr2:=Br[1]:Br:=Extract(Br,[Br[1]]):T:=CheckBranch(chBr2,R):
  if T=false then 
NP:=[op(NP),chBr2]:
  else chBr:=T:
   if chBr[3]<>[] then
     q:=DeltaPolynomial(op(chBr[3][1]),R):
     q:=[DifferentialPrem(q,[seq(seq(parms[i][vars[j]],j=1..nops(vars)),i=1..nops(parms))],F)][2]:
     
     q:=expand(ParmDiffRed(q,chBr[1],reduction=full,R)):
     if chBr[6]<>[] then
      q:=Groebner[NormalForm](q,chBr[6],plex(op(parms))):
     fi:
     q:=SIM(q):
     if type(q,'constant') and q<>0 then NP:=[op(NP),chBr]:
     elif q=0 then Br:=[[chBr[1],chBr[2],Extract(chBr[3],[chBr[3][1]]),chBr[4],chBr[5],chBr[6],chBr[7]],op(Br)]:
     else 
      q:=DIVIDE(q,[op(ADD(chBr[5],chBr[4])),op(chBr[7])]):
      if type(q,'constant') then NP:=[op(NP),chBr]:
      else
       chBr1:=[chBr[1],chBr[2],Extract(chBr[3],[chBr[3][1]]),chBr[4],chBr[5],chBr[6],chBr[7]]:
       Br4:=[MakeTree(q,chBr1,R)]:
       Br:=[op(Br4),op(Br)]:
      fi:
    fi:
  elif chBr[2]<>[] then
    q:=chBr[2][1]:
    if chBr[6]<>[] then
      q:=Groebner[NormalForm](q,chBr[6],plex(op(parms))):
    fi:
    q:=expand(ParmDiffRed(q,chBr[1],reduction=full,R)):q:=SIM(q):
    if type(q,'constant')=true and q<>0  then NP:=[op(NP),chBr]:
    elif q=0 then Br:=[[chBr[1],Extract(chBr[2],[chBr[2][1]]),chBr[3],chBr[4],chBr[5],chBr[6],chBr[7]],op(Br)]:
    else
      q:=DIVIDE(q,[op(ADD(chBr[5],chBr[4])),op(chBr[7])]):
      if type(q,'constant') then NP:=[op(NP),chBr]:
      else
       chBr1:=[chBr[1],Extract(chBr[2],[chBr[2][1]]),chBr[3],chBr[4],chBr[5],chBr[6],chBr[7]]:
       Br4:=[MakeTree(q,chBr1,R)]:
       Br:=[op(Br4),op(Br)]:
     fi:
   fi:
 else
    Decom:=[op(Decom),chBr]: 
 fi:
fi:
 end do:

de:=[checkregular([seq([Decom[i][1],ADD(Decom[i][4],Decom[i][5]),Decom[i][6],Decom[i][7]],i=1..nops(Decom))],R)]:
return(de):
end: