restart: 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] ): num:=proc(f,L) if member(f,L,'q') then return(q): fi: return(false): 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: #................................................................................................ FACTOR:=proc(f) local i,L,g: if type(f,`+`) then L:=[op(f)]: g:=GCD(L): return(g*(add(L[i]/g,i=1..nops(L)))): else return(f): fi: end: #.................................................................................................. #it simplifies the coeffs of f SIM:=proc(f) if f=0 then return(f):else return(f/GCD([coeffs(expand(f))])):fi: 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 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: #...................................................................................................... #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 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: #............................................................................................................. #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, 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: # 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]):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)]: Gnew:=[op(Gnew),f]: return(Gnew,Ex,Bnew): end: CH:=proc(f,S,R) local i,s: i:=SIM(Initial(f,R)): s:=SIM(Separant(f,R)): if {type(i,'constant'),type(s,'constant')}={true} then return(true): elif type(i,'constant')=true and Extract(S,[s,-s])<>S then return(true): elif type(s,'constant')=true and Extract(S,[i,-i])<>S then return(true): elif {type(i,'constant'),type(s,'constant')}={false} and Extract(S,[s,-s])<>S and Extract(S,[i,-i])<>S then return(true): else return(false): fi; end: # Find-New-Condition sub algorithm................................................................................. #f is a new polynomial and OC is the list of old conditions NewCondition:=proc(f,oc,p,R) local init,inn,s,sep,NC,E,OC,F,e: #option trace: E:=[]:F:=p: OC:=ADD(oc,p): inn:=Initial(f,R): s:=Separant(f,R): init:=DIVIDE(SIM(inn),OC): if type(init,'constant')=false then if {member(init,OC),member(-init,OC)}={false} then if CH(init,OC,R)=false or member(0,[seq([DifferentialPrem(OC[i],[init],reduction=full,R)][2],i=1..nops(OC))])=false then E:=[init]:OC:=ADD(OC,[init]): e:=i:else F:=[op(F),init]:fi: fi: fi: sep:=DIVIDE(expand(SIM(s)),OC): if type(sep,'constant')=false then if {member(sep,OC),member(-sep,OC)}={false} then if CH(sep,OC,R)=false or member(0,[seq([DifferentialPrem(OC[i],[sep],reduction=full,R)][2],i=1..nops(OC))])=false then E:=ADD(E,[sep]):e:=se else F:=[op(F),sep]:fi: 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 #Re:=update(O[i],P[1],B,R): #TR:=Re[1]: #B:=Re[3]: 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]]]: 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]]: fi: else NC:=NewCondition(f,P[4],P[5],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]: elif nops(NC[2])=1 then P4:=P[4]: 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[1],RZ(ADD(Extract(P[2],[f,-f]),[NC[2][1],SIM(f1)])),Extract(P[3],[f,-f]),P[4],Z]: else P4:=P[4]: for i from 1 to nops(P[4]) do if [DifferentialPrem(NC[2][1],[P[4][i]],reduction=full,R)][2]=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 [DifferentialPrem(NC[2][2],[P[4][i]],reduction=full,R)][2]=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[1],RZ(ADD(Extract(P[2],[f,-f]),[ NC[2][1],SIM(f1)])),Extract(P[3],[f,-f]),P[4],Z],[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]: fi: fi: return(NewP): end: #Differential-consistency subalgorithm............................................................................ CheckBranch:=proc(P,R) local i,sw,j,L,S,L1,Z,F: #option trace: S:={seq(type(P[1][i],'constant'),i=1..nops(P[1])),seq(type(P[2][i],'constant'),i=1..nops(P[2])),seq(type(P[3][i],'constant'),i=1..nops(P[3]))}: if member(true,S) then return(false): elif P[4]=[] and P[5]=[] then return(P): else Z:=P[4]: L:=[seq(DIVIDE(SIM(expand([DifferentialPrem(Z[i],P[1],reduction=full,R)][2])),[seq(Z[j],j=i+1..nops(Z))]),i=1..nops(P[4]))]: if member(0,L) then return(false): else L1:=L: for i from 1 to nops(L) do if type(L[i],'constant') then L1:=Extract(L1,[L[i]]):fi: od: F:=convert({op(P[4])} minus {op(L1)},list): return([P[1],P[2],P[3],L1,ADD(P[5],F)]):fi: fi: end: #Rosenfeld's lemma............................................................................. #it detects consistent regular systems by rosenfeld lemma with(PolynomialIdeals): chechregular := 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: #................................................................................................................. #Main algorithm MainProc:=proc(P,R) local inn,i,cn,cm,Br,NP,NewB,NC,chBr,q,Decom,q1,q2,chBr1,chBr2,T,count,Br4,cc,BBB,a,b,de,BB,dee,fe: #option trace: i:=1:count:=0: Decom:=[]: NP:=[]: Br:=[[[],P,[],[],[]]]:cc:=1: 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 q1:=DeltaPolynomial(op(chBr[3][1]),R):count:=count+1: q2:=expand([DifferentialPrem(q1,chBr[1],reduction=full,R)][2]):q:=SIM(q2): 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]],op(Br)]:else q:=DIVIDE(q,ADD(chBr[5],chBr[4])):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]]: Br4:=[Branch(chBr1,q,R)]: Br:=[op(Br4),op(Br)]:fi:fi: elif chBr[2]<>[] then q1:=chBr[2][1]:q2:=expand([DifferentialPrem(q1,chBr[1],reduction=full,R)][2]):q:=SIM(q2): if type(q,'constant')=true and q<>0 then NP:=[op(NP),chBr]:elif q=0 then Br:=[[chBr[1],Extract(chBr[2],[q1]),chBr [3],chBr[4],chBr[5]],op(Br)]:else q:=DIVIDE(q,ADD(chBr[5],chBr[4])):if type(q,'constant') then NP:=[op(NP),chBr]:else chBr1:=[chBr[1],Extract(chBr[2],[q1]),chBr[3],chBr[4],chBr[5]]: Br4:=[Branch(chBr1,q,R)]: Br:=[op(Br4),op(Br)]:fi:fi: else Decom:=[op(Decom),chBr]: fi: fi: end do: dee:=[seq([Decom[i][1],[op({seq(op([SIM(Initial(Decom[i][1][j],R)),SIM(Separant(Decom[i][1][j],R))]),j=1..nops(Decom[i][1]))} minus {1})]],i=1..nops(Decom))]: fe:=dee: for a in dee do for b in Extract(dee,[a]) do if {op(a[1])} intersect {op(b[1])}={op(a[1])} then fe:=Extract(fe,[b]): fi: end do: end do: de:={chechregular(fe)}: print(nops(NP)+nops(Decom)-nops(de),count): return(de): end: