with(PolynomialIdeals):
with(DifferentialAlgebra):with(Tools):
R:= DifferentialRing ( derivations = [x,y,z], blocks = [[u,v,w,m]] ):
H:= DifferentialRing ( derivations = [x,y,z], blocks = [m,w,v,u] ):
#divide 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:
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:
#......................................................................................................
#add 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 {seq(is(M[i]=s[j]),j=1..nops(s))}={false} then s:=[op(s),M[i]]:fi:
end:
return(s):
fi:
end:


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


#simplify the coeffs of f
SIM:=proc(f)
if f=0 then return(f):else return(f/GCD([coeffs(expand(f))])):fi:
end:
#SIM(-20*(u-1)^2):simplify(-20*(u-1)^2):

delta_leader:=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:

autopr := proc (A, H, R)
local B, L, a, u, b, HA, HB, H1, K:
#option trace:
 B := {NULL}:
 L := SortByRank(A,ascending,R):
 for a in L do 
  u := LeadingDerivative(a, R):
  b := [DifferentialPrem(a,[op(B)],R)][2]:
  if LeadingRank(expand(b), R) = LeadingRank(a, R) then 
   B := B union {b}:
  else 
  return({}):
  end if:
 end do:
 HA := {seq(Initial(a, R), a in A), seq(Separant(a, R), a in A)}:
 HB := {seq(Initial(a, R), a in B), seq(Separant(a, R), a in B)}:
 H1 := {op(H)} minus HA:
 K := HB union {seq(expand([DifferentialPrem(h, [op(B)],reduction=partial,R)][2]), h in H1)}:
 
  if member(0, K) then
  return({}):
 else 
 K := remove(type, K, 'realcons'):
  return({[[op(B)], [op(K)]]}):
 end if:
 end proc:

#autopr([-u*v[x, x]+2*u[y]*v[y]+v[x, x], -u, v[y]^2], [2*u[x], -u+1, 2*v[y]],R):

Test:=proc(f)
local H:
#option trace:
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:
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:

#DifferentialPrem(-u*v[x, x]+2*u[y]*v[y]+v[x, x], [2],R):
#second criterian of bukh.............................................................................
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, y],u[x]*v[y], v*u[x, x]^2+v*u[x, x]+u[x]],R):
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:

 
complete := proc (T, f, R)
 local a, uu, t,GA, i, A1, G1, D1, A, uuu, H1,DD,inn,countD,C,Bnew,D,E,s,hd,k,B,d,ldk1,ldk2:
#option trace:
C:=[]:
 uu :=LeadingDerivative(f, R):
 GA := []:
 for i in T[3] do
   if delta_leader(i, uu, R) <> false then 
     GA := ADD(GA,[i]):
   end if:
 end do:
 A1 := Extract(T[3],GA):
 d:=[FactorDerivative(SIM(LeadingDerivative(f,R)),R)]:
G1:=T[1]:
 D1 := T[2]:
for i from 1 to nops(T[3]) do
if [FactorDerivative(SIM(LeadingDerivative(T[3][i],R)),R)][2]=d[2] then
C:=ADD(C,[{SIM(T[3][i]),f}]):
fi:
od:
t:=Test(f,R):
E:=[]:
while C<>[] do
  s:=C[1]:C:=Extract(C,[C[1]]):
  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:=T[2]:i:=1:
while i<=nops(B) do
  k:=B[i]:B:=Extract(B,[k]):SC([k[1],f,k[2]],R):
  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)]:
inn:=remove(type,[Initial(f, R), Separant(f, R)],constant):
 H1 := ADD(T[4],inn):
 return([[op(G1)], Bnew, [op(A1), f], [op(H1)]]):
end proc:

RR := proc (X)
 local L,a,b,B,A,m,C:
#option trace:
 L := NULL:
 for a in X do
  b := [op(a[1]), op(a[2])]:
  B := SortByRank(indets(b),descending,R):
  m := nops(B):
  A := subs(seq(B[j] = N[j], j = 1 .. m), a):
  C := [op(A[1]),seq(g[i]*A[2][i]-1,i=1..nops(a[2]))]:
 if not IdealMembership(1, EliminationIdeal(<(op(C))>, {seq(N[i],i=1...m)})) then
   L := L, a:
 end if:
 end do:
  return(L):
 end proc:

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:



RosGrobim := proc (F, K, R)
 # option trace:
 local S, A, T, S1, p, G1, D1, p1, p2,L, p3, GD,count,countD,D3,Au,non,TT,de,BBB,BB,b,a,Bb,fe,gf,dee,fee:
 S := {[F, [], [], K]}: 
 A := {}:countD:=0:non:=[]:
 while S <> {} do 
  T := S[1]:
  S1 := S minus {T}:
  if T[1]= [] and T[2]=[] then
    Au:=autopr(T[3], T[4], R): 
    if Au={} then non:=[op(non),T]: else
    A := A union Au:
    fi: 
  else
   L:=[]:
 for a in T[4] do
  if type(a,'constant')=false then L:=[op(L),[DifferentialPrem(a,T[3],R)][2]]:fi:
  end do:
  if member(0,L) then non:=[op(non),T]:
  else
   G1:=T[1]:
   D1:=T[2]:
   if T[1]<>[] then 
    p := T[1][1]:G1:= Extract(G1,[p]):
   else
    p:=DeltaPolynomial(op(D1[1]),R):countD:=countD+1:
    D1:=Extract(D1,[D1[1]]):
    fi:
    TT:=RZ(T[3]):
    p1 := SIM(expand([DifferentialPrem(SIM(p), [op(TT)],R)][2])): 
    p1:=DIVIDE(p1,T[4]):
    if p1 = 0 then
      S1 := S1 union {[G1, D1, T[3], T[4]]}:count:=count+1:
    elif type(p1, realcons) then non:=[op(non),T]:
    else
     p2 := SIM(expand(p1-Initial(p1, R)*LeadingDerivative(p1,R)^degree(p1,LeadingDerivative(p1,R)))):
     p3 := SIM(expand(degree(p1, LeadingDerivative(p1, R))*p1-LeadingDerivative(p1, R)*Separant(p1, R))):
     D3:=[complete([G1, D1, T[3], T[4]], p1, R)]:
     S1 := S1 union {D3[1]}:
     if type(Separant(p1, R), realcons) = false then
      if type(Initial(p1, R), realcons) = false then
        S1 := S1 union {[ADD(G1,[SIM(Separant(p1, R)), p3]), D1, T[3], [op(T[4]), SIM(Initial(p1, R))]]}:
      else 
        S1 := S1 union {[ADD(G1,[SIM(Separant(p1, R)), p3]), D1, T[3], [op(T[4])]]}:
        count:=count+1:
      end if
     end if:
     if type(Initial(p1, R), realcons) = false then
       S1 := S1 union {[[op(G1), p2, SIM(Initial(p1, R))], D1, T[3], T[4]]}:count:=count+1:
     end if
    end if:
   end if: 
   end if:
  S := S1:

 end do:
fe:=[]:
for a in A do
 gf:=[seq(SIM(DIVIDE(a[1][i],a[2])),i=1..nops(a[1]))]:
if member(1,gf)=false then
 fe:=[op(fe),[gf,a[2]]]:
fi:
end do:
dee:=[seq([fe[i][1],[op({seq(op([SIM(Initial(fe[i][1][j],R)),SIM(Separant(fe[i][1][j],R))]),j=1..nops(fe[i][1]))} minus {1})]],i=1..nops(fe))]:
fee:=dee:
for a in fe do
 for b in Extract(fe,[a]) do
   if {op(a[1])} intersect {op(b[1])}={op(a[1])} then 
    fee:=Extract(fee,[b]):
   fi:
 end do:
end do:
de:=[RR(fee)]:
#print(nops(non)+nops(A)-nops(de),countD):
 return(de):
 end proc: