with(PolynomialIdeals): with(diffalg): R := differential_ring(ranking = [u, v, w], derivations = [x, y]): sorth := proc (a, b, R) if greater(a, b, R) = true then return(true): elif greater(b, a, R) = true then return(false): else return(sorth(initial(a, R), initial(b, R), R)); end if; end proc: f1 := proc (p, R, A) local q,c,flag,i,d: q := denote(p, jet, R); if type(q, realcons) = true then return(false): else c := denote(leader(q, R), jet, R); flag := false: i := 1: while i <= nops(A) and flag = false do d := delta_leader(c, denote(leader(A[i], R), jet, R), R); if d <> FAIL and d <> NULL then flag := [truepar, i, [d], A[i], c]: elif d = NULL and degree(A[i], c) <= degree(q, c) then flag := [truealg, i]: else i := i+1: end if: end do: if flag = false then return(f1(initial(q, R), R, A)): else return(flag): end if: end if: end proc: f4 := proc (p, R, k) local q, v, s, d, a3, T, W, n, hh, N, F, h, i; q := p; if type(q, realcons) = true then return(q): end if; v := k[5]; s := separant(k[4], R); d := degree(q, denote(v, jet, R)); a3 := differentiate(k[4], op(k[3]), R); T := s*v-a3; h := 0; W := [seq(coeff(q, v, i), i = 0 .. d)]; return add(T^i*W[i+1]*s^(d-i), i = 0 .. d); end proc: pdred := proc (p, R, A) #option trace; local q, r, F, k; q := denote(p, jet, R); r := 0; if type(expand(q), `+`) = true then F := [op(expand(q))]: else F := [q]: end if; F := sort(F, (a,b)->sorth(a,b,R) ); while q <> 0 do k := f1(F[1], R, A); #print(k); if k[1] = truepar then q := f4(q, R, k): else r := r+F[1]; q := q-F[1]: end if; if type(expand(denote(q, jet, R)), `+`) = true then F := [op(expand(q))] else F := [q] end if; F := sort(F, (a,b)->sorth(a,b,R)): end do; return(expand(denote(r, jet, R))) end proc: dred := proc (p, R, A) #option trace; local q, r, F, k, i; q := denote(p, jet, R); r := 0; if type(expand(q), `+`) = true then F := [op(expand(q))]: else F := [q]: end if; #print(salam1); #F := sort(F, (a,b)->sorth(a,b,R)); while q <> 0 do #print(salam2); k := f1(F[1], R, A); if k<>false then #print(k); fi: if k[1] = truepar then q := f4(q, R, k): elif k[1] = truealg then i := initial(A[k[2]], R); q := simplify(i*q-simplify(denote(A[k[2]]*F[1], jet, R)/rank(A[k[2]], R))): else r := r+F[1]; q := q-F[1]: end if; if type(expand(denote(q, jet, R)), `+`) = true then F := [op(expand(q))]: else F := [q]: end if; end do; return(expand(denote(r, jet, R))): end proc: autopr := proc (A, H, R) local B, L, a, u, b, HA, HB, H1, K; B := {NULL}; L := sort(A, (a, b)->greater(b,a,R)): for a in L do u := leader(a, R); b := pdred(a, R, [op(B)]); if rank(b, R) = rank(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(pdred(h, R, [op(B)]), h in H1)}; K := remove(type, K, 'realcons'); if member(0, K) then return({}): else return({[[op(B)], [op(K)]]}): end if: end proc: lcd := proc (p, q, R) #option trace; local u1, u2, u3, a1, a2, a,N; u1 := leader(p, R); u2 := leader(q, R); u3 := op(0, u1); a1 := [delta_leader(p, u3, R)]; a2 := [delta_leader(q, u3, R)]; N:=NULL; for a in [a1] do if member(a,a2,'k') then N:=N,a: a2:=subsop(k=NULL,a2): fi: od: if N=NULL then return(1): else return(differentiate(u3, N , R)); end if: end proc: lh := proc (q, R) #option trace; local u, L, a; if degree(q)<>1 then RETURN(false); fi: u := op(0, leader(q, R)); L := NULL; if type(q, `+`) then L := [op(q)]; else L := [q]; end if; for a in L do if op(0, simplify(a/coeffs(a))) <> u then return(false): end if; end do; return(true): end proc: b3 := proc (p, R) #option trace: local pp, i,p1; p1:=factor(p); pp := [op(p1)]; if not lh(p1,R) and not type(p1,`*`) then return(false): fi: for i in pp do if type(i, `^`) = true and not type(a,`realcons`) then if lh(op(i)[1], R) = false then return(false): fi: elif type(i, `+`) = true and not type(a,'realcons') then if lh(i, R) = false then return(false): fi: fi: end do; return(true): end proc: b4:=proc(p,a,R) #option trace; if b3(p,R) and b3(a,R) and lcd(leader(p,R),leader(a,R),R)=1 then return(true): else return(false): fi: end: update := proc (T, p, R) local a, uu, GA, i, A1, G1, D1, A, uuu, H1; uu := leader(p, R); GA := {}; for i in T[3] do if delta_leader(i, uu, R) <> FAIL then GA := GA union {i}: end if: end do; A1 := {op(T[3])} minus GA; G1 := {op(T[1])} union GA; D1 := {op(T[2])}; for a in A1 do uuu := leader(a, R); if op(0, uu) = op(0, uuu) and not b4(p,a,R) then D1 := D1 union {delta_polynomial(p, a, R)}: end if: end do; D1 := D1 minus {0}; H1 := {op(T[4])} union {initial(p, R), separant(p, R)}; return([[op(G1)], [op(D1)], [op(A1), p], [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 := indets(b); m := nops(B); A := subs(seq(B[i] = t[i], i = 1 .. m), a); C := [op(A[1]), seq(1-T[i]*A[2][i], i = 1 .. nops(A[2]))]; if not IdealMembership(1, EliminationIdeal(<(op(C))>, {seq(t[i], i = 1 .. m)})) then L := L, a: end if: end do; return(L): end proc: RosGrobim := proc (F, K, R) #option trace; local S, A, T, S1, p, G1, D1, p1, p2, p3, GD; S := {[F, [], [], K]}; A := {}; while S <> {} do T := S[1]; S1 := S minus {T}; GD := {op(T[1])} union {op(T[2])}; if GD = {} then A := A union autopr(T[3], T[4], R): else p := GD[1]; G1, D1 := {op(T[1])} minus {p}, {op(T[2])} minus {p}; p1 := differential_sprem(p, T[3],R); if p1 = 0 then S1 := S1 union {[[op(G1)], [op(D1)], T[3], T[4]]}: elif type(p1, realcons) = false then p2 := simplify(p1-initial(p1, R)*rank(p1, R)); p3 := simplify(degree(p1, leader(p1, R))*p1-leader(p1, R)*separant(p1, R)); S1 := S1 union {update([[op(G1)], [op(D1)], T[3], T[4]], p1, R)}; #print(sa); if type(separant(p1, R), realcons) = false then if type(initial(p1, R), realcons) = false then S1 := S1 union {[[op(G1 union {separant(p1, R), p3})], [op(D1)], T[3], [op(T[4]), initial(p1, R)]]}: #print(sal); else S1 := S1 union {[[op(G1 union {separant(p1, R), p3})], [op(D1)], T[3], [op(T[4])]]}: #print(sala); end if end if; if type(initial(p1, R), realcons) = false then S1 := S1 union {[[op(G1), p2, initial(p1, R)], [op(D1)], T[3], T[4]]}: #print(salam); end if end if; end if; S := S1: end do; return(RR(A)); end proc: #RosGrobim([(u[x, x, x, x]-u[])^2*(u[x, x]+u[x])^3, (u[y, y, y]-u[y, y])^2*(u[y]+u[])^3], [], R)