with(combinat); #-------- splitting_field: read(diffop); # Compute the splitting field of a local differential operator. # input: local differential operator f # output: [ ext, ram , ff] = [algebraic extension over Q , ramification , f # with the ramification applied on it ] # comment: it's not at all efficient, only to be used for operators of # small order spl_field:=proc(f,ext) local v,i,r,ff; v:=factor_op(f,`all alg factors`,ext); for i in v do if has(i,`alg over k((x))`) then if i[2]<>ext then RETURN(spl_field(f,i[2])) fi; # now i[2]=ext, i.e. no alg extension, hence a ramification: r:=i[3]; ff:=ramification_of(f,r,ext); v:=spl_field(ff,ext); RETURN([v[1],ramification_of(v[2],r,ext),ff]) fi od; # No algebraic extensions nor ramifications: [ext,x,f] end: # Input: f,ext such that there are no algebraic extensions needed to compute # a basis of solutions of f. # Output: a list containing: # variable names, each name corresponds to a solution # a list A such that subs(A,var) gives the exponential part # a list B such that subs(B,var) gives the rest of the value of that solution formal_sols:=proc(f,ext) options remember; local v,x,subep,subv,i,j,h; v:=formal_sol(f,ext); subep:=NULL; subv:=NULL; i:=1; for j in v do for h in j[1] do subep:=subep,x[i]=j[2]; subv:=subv,x[i]=h; i:=i+1; if nops(j)>2 then ERROR(`wrong input`) fi od od; # now the value for x[i] is x^subs(subep,x[i])*subs(subv,x[i]) res:=[[seq(x[i],i=1..degree(f,xDF))],[subep],[subv]]; userinfo(5,diffop,`formal solutions: `,res): res; end: # Input: an integer d # Output: a list of monomials of degree d which have an integer valuation # and their values. int_combs_sol:=proc(f,ext,d) local v,vs,res,i,ep,j,h; v:=formal_sols(f,ext); vs:=symmetric_sums(v[1],d); res:={}; for i in vs do ep:=subs(v[2],i); if type(ep,integer) then # convert the sum i to a product j: j:=convert([seq(h^coeff(i,h),h=indets(i))],`*`); res:=res union {j=eval_laurent(expand(x^ep*subs(v[3],j)),ext)} fi od; userinfo(3,diffop,`number constants to be computed: `, nops(res)); res end: # Not good with modm in the original code. lift_sol1:=proc(LL,aa,n,ext,acc) local L,a,t,i; L:=value_laurent[LL]; a:=nmterms_laurent(aa,acc-n); for t from accuracy_laurent[LL] to acc-1 do # divide the coefficient of x^t in a*L by n-t: L:=L+modm(g_expand(convert([seq(coeff(L,x,i)*coeff(a,x,t-i) ,i=ldegree(L,x)..degree(L,x))],`+`)/(n-t)*x^t,ext)) od; L end: sol_sympower:=proc(ff,d,N) local f,ext,b,i,v,r,s,l,j; global modulus,modulus2; # options trace; f:=convert(ff,diffop,x=0); ext:=g_ext(f); # v:=spl_field(f,ext); # For now I'll assume that no algebraic field extensions are needed b:=global_bounds(ff,d); if b[1]=0 then RETURN([0,0]) fi; s:=int_combs_sol(f,ext,d); if s={} then bug() fi; r:=convert([seq(v[i]*x^i,i=0..b[2])],`+`); i:=b[2]; l:=0; lv:=0; userinfo(3,diffop,`now computing the system for the coeffs of the invariant`); for j in s do i:=i+1; l:=l+v[i]*op(j)[1]; lv:=lv+v[i]*eval_laurent(expand(op(j)[2]/eval_laurent(b[1],ext)),ext) od; vars:={seq(v[j],j=0..i)}; # Now solve l=r in the variables vars # modulus2 is a prime and modulus is a power of modulus2 modulus2:=nextprime(1000000); zero:=compute_modm(2,nt,lv-r,N); # modular version of next line: # zero:=nt(lv-r,N); zero:={coeffs(zero,log(x))}; zero:={seq(coeffs(i,x),i=zero)}; userinfo(3,diffop,`computed the system for the coeffs of the invariant`); userinfo(3,diffop,`bow solving this system`); # solve needs modulus as well: modulus:=modulus2^2; result:=subs(g_solve(zero,vars),[l,r*b[1]]); modulus:=0; result end: #f:=1/9-94217/131072*x*DF+5482559/65536*x^2*DF^2+6384395/8192*x^3*DF^3 #+6135425/4096*x^4*DF^4+62285/64*x^5*DF^5+8161/32*x^6*DF^6 #+55/2*x^7*DF^7+x^8*DF^8; #f:=collect(f/lcoeff(f,DF),DF,factor); ## trace(bound_valuation_solsym,spl_field); ## global_bounds(f,2,[]); #------------------------------------------ sympowmatrix # From JA Weil: sympowmatrix:=proc(n,m,mat) local i,j,k,l,v,var,mini,maxi,N,SMU,w,M,liste; global _e,_y: # options remember; # That causes troubles, because if you change the output # matrix of this procedure, then call it again, you get the modified # matrix instead of the good one. userinfo(3,diffop,`started computing the sympower matrix`): if nargs=3 then for i from 1 to n do for j from 1 to n do _y[i][j]:=mat[i,j] od: od: else _y:=evaln(_y); fi; _e:=evaln(_e); for i from 1 to n do v[i]:=convert([seq(_y[i][j]*_e[j],j=1..n)],`+`); od; var:=[seq(_e[i],i=1..n)]; mini:=combinat[vectoint]([m,0$(n-1)]): maxi:=combinat[vectoint]([0$(n-1),m]): N:=maxi-mini+1: # = order of the m'th sym power of an n'th order operator SMU:=linalg[matrix](N,N): for i from mini to maxi do w:=combinat[inttovec](i,n); M:=collect(convert([seq(v[j]^w[j],j=1..n)],`*`),var); liste:=[multihom(M,var,m)]; for j from 1 to N do SMU[i+1-mini,j]:=liste[j] od; od; userinfo(3,diffop,`computed the sympower matrix`): evalm(SMU); end: multihombis:=proc(expr,var,m) local i,j,k,varbis,nieuw; if m=0 then RETURN(expr) fi; if nops(var)=1 then RETURN(subs(var[1]=1,expr)); fi; varbis:=[seq(var[i],i=2..nops(var))]; seq(multihombis(coeff(expr,var[1],m-i),varbis,i),i=0..m); end: multihom:=proc(expr,var,m) multihombis(collect(expr,var),var,m); end: #-------------------------------------------- # Assumptions for now: # the singularity we use is located at x=p, for now only p=0 # no ramifications at this singularity (but they are allowed at other singularities) # all generalized exponents at x=p are in Q(ext)[x^(-1)], # here ext gives the algebraic extension # eval_semireg_monomial: check if a monomial is semiregular (hence # in C((x))[log(x)] because we assumed no # ramifications), and if so, evaluate it mod x^n # derivative : Compute n'th derivative of an expression # [e,s]=Exp(e)*s # # Assumption: no ramifications # Input: a monomial M, # the output of formal_sols---NO, bit different # and an integer n # Output: value of the monomial mod x^n, or `non-trivial ep` if this value # is not in C((x))[log(x)] # sols[1] = list of vars # sols[2] = ep (exponential part) # sols[3] = semireg part of the solutions, must have valuation >=0 eval_semireg_monomial:=proc(M,sols,n) # The following isn't quite right yet, we must assume that either # all terms, or no terms have a trivial exp part if type(M,`+`) then RETURN(convert([seq(procname(i,args[2..nargs]),i=op(M))],`+`)) fi; ep:=convert([seq(degree(M,sols[1][i])*sols[2][i],i=1..nops(sols[1]))],`+`); if not type(ep,integer) then RETURN(`non-trivial ep`) fi; s3:=[seq(sols[3][i]$degree(M,sols[1][i]),i=1..nops(sols[1]))]; # Accuracy that is required from the semiregular parts of each # solution is n-ep, under the assumption that the valuations of these # semireg parts are >=0. s3:=nt(s3,n-ep); #lprint(hello_there); R:=expand(lcoeff(M)*convert(subs(hello=ln(x), taylor(subs(ln(x)=hello,convert(s3,`*`)),x=0,n-ep)),polynom)*x^ep); #lprint(bye_for_now); R end: # Assumption: no ramifications # input [e,s], meaning Exp(e)*s # optional n: take the n'th derivative # Output: the [e,s] of the derivative v(s)>=0. derivative:=proc(es,n) options remember; if nargs>1 then if n=0 then RETURN(es) elif n>0 then RETURN(procname(procname(es,n-1))) else ERROR(`wrong arguments`) fi fi; w:=es[2]; ep:=es[1]; w:=eval_laurent(expand(differentiate(w)+w*ep/x)); m:=min(seq(valuation_laurent(i,0),i=indets(w) intersect set_laurents)); if m<0 then ep:=ep+m; w:=eval_laurent(expand(w*x^(-m))) fi; [ep,w] end: # Compute the symmetric power matrix modulo x^ac eval_sympowmatrix:=proc(m,sols,ac) s1:=sols[1]; n:=nops(s1); s2:=subs(sols[2],s1); s3:=subs(sols[3],s1); s1:=[seq(seq(_y[i][j],i=1..n),j=1..n)]; s23:=[seq(seq(derivative([s2[i],s3[i]],j-1),j=1..n),i=1..n)]; s2:=[seq(i[1],i=s23)]; s3:=[seq(i[2],i=s23)]; M:=sympowmatrix(n,m); for j from 1 to linalg[coldim](M) do # lprint(j); userinfo(3,diffop,`checking exp part of collumn `,j): M[1,j]:=eval_semireg_monomial(expand(M[1,j]),[s1,s2,s3],ac) od; for j from 1 to linalg[coldim](M) do if not has(M[1,j],`non-trivial ep`) then # lprint(j); userinfo(3,diffop,`modular computing of column `,j): for i from 2 to linalg[coldim](M) do M[i,j]:=eval_semireg_monomial(expand(M[i,j]),[s1,s2,s3],ac) od fi od; evalm(M) end: if false then fl:=l_p(f,0); sols:=formal_sols(fl,[]); # trace(eval_sympowmatrix,eval_semireg_monomial,derivative); WW:=eval_sympowmatrix(6,sols,13): w:=NULL; for i from 1 to linalg[coldim](WW) do if not has(WW[1,i],`non-trivial ep`) then w:=w,WW[1,i] fi od; w:=[w]; Q:=convert([seq(a[i]*w[i],i=1..nops(w))],`+`); quo(Q,x^5,x); S:=solve({coeffs(expand("),x)}); subs(S,Q); fi; # printlevel:=5; if false then # From the file: hurwitz # G_168 deq:=x^2*(x-1)^2*diff(y(x),x$3)+(7*x-4)*x*(x-1)*diff(y(x),x$2)+(72/7*(x^2-x)-20/ 9*(x-1)+3/4*x)*diff(y(x),x)+(792/7^3*(x-1)+5/8+2/63)*y(x)=0; f:=convert(deq,diffop); f0:=l_p(f,0); sols:=formal_sols(f0,[]); M:=eval_sympowmatrix(6,sols,13): time(); # Verdraaid, dat kost 5142 seconden... # Nog eens, nu alleen eval_semireg_monomial aanroepen als # je geen `non-trivial ep` hebt. 4715 seconden, maakt dus niet veel uit. fi; if false then # From the file test_spl: F:=xDF^2*(xDF-1/2)*(xDF+1/(12))*(xDF-7/(12))*(xDF-1/(12))*(xDF-5/(12)) *(xDF+1/2)-x*(xDF+3/4)*(xDF+5/4); G:=collect(F,xDF); f:=convert(G,diffop,monic); f0:=l_p(f,0); # sols:=formal_sols(f0,[]); # M:=eval_sympowmatrix(2,sols,5): # time(); # 920 seconds # M:=eval_sympowmatrix(3,sols,5): # Feasible, but takes a long time. fi; #---------------------------------------- global_bounds_new # input: a global differential operator f # output: [r, bound] (with alg numbers in RootOf notation) # such that all rational solutions of symmetric_power(f,n) are of the form # r * polynomial of degree <= bound global_bounds:=proc(f,n,ext) local i,singularities,p,b,bound,d; if nargs=2 then i:=g_ext(f); RETURN(procname(subs(g_conversion1,f),n,i)) fi; # [[point,alg ext],..] singularities:=[seq([i[1],[op(i[3]),op(ext)]],i=[[infinity,[],[]], op(v_ext_m(g_factors(denom(g_normal(f/lcoeff(f,DF))),x,ext),x))])]; for i in singularities do userinfo(3,diffop,`computing bounds for singularity `,i): b:=bound_valuation_solsym(convert(f,diffop,x=i[1]),n,i[2]); if b=infinity then RETURN([0,0]) fi; if i[1]=infinity then # this one comes first. bound:=-b; if not type(n,integer) then # n is a list bound:=bound- convert([seq(n[p+1]*p*2,p=1..nops(n)-1)],`+`) fi; p:=1 elif i[2]=ext then p:=p*(x-i[1])^b; bound:=bound-b; else d:=subs(_Z=x,op(subs(g_conversion2,i[2][1]))); p:=p*d^b; bound:=bound-degree(d,x)*b fi od; userinfo(3,diffop,`ended computing bounds for singularities`): userinfo(3,diffop,`found `, [p,bound]): if bound<0 then RETURN([0,0]) fi; [p,bound] end: # a lower bound on the valuation of solutions in k((x)) of the nth symmetric # power of f. # If n is not an integer but a list n=[n_0,n_1,..n_k] then compute the bound # for symmetric_product(f_0$n_0,..,f_k$n_k) where f_i stands for the operator # for which the solutions are the i'th derivatives of the solutions of f. bound_valuation_solsym:=proc(f,n,ext) local v,e,i,m; if nargs<=1 then ERROR(`too few arguments`) elif nargs=2 then RETURN(procname(f,n,g_ext(f))) fi; v:=spl_field(f,ext); # v[3] = f in the splitting field # The minimal generalized exponents of v[3]: e:={seq(min(solve(nt(i[1],1)))+i[4],i=factor_op(v[3],`semireg`,v[1]))}; e:=symmetric_sums(e,n,degree(v[2],x)); m:=infinity; for i in e do if type(i,integer) then m:=min(m,i) fi od; if m=infinity then m else ceil(m/degree(v[2],x)) fi end: # used in bound_valuation_solsym(f,n,ext) # r = ramification index symmetric_sums:=proc(e,n,r) local e_prime,i; if n=0 then {0} elif n=1 then e elif type(n,integer) then symmetric_sum(e,procname(e,n-1)) elif nops(n)=1 then procname(e,op(n)) else e_prime:={seq(i+ldegree(i,x)-r,i=e)}; symmetric_sum(procname(e,n[1]), procname(e_prime,[op(2..nops(n),n)],r)) fi end: symmetric_sum:=proc(a,b) {seq(seq(i+j,i=a),j=b)} end: # Compute the symmetric power matrix modulo x^ac eval_sympowmatrix:=proc(m,sols,f) s1:=sols[1]; n:=nops(s1); s2:=subs(sols[2],s1); s3:=subs(sols[3],s1); s1:=[seq(seq(_y[i][j],i=1..n),j=1..n)]; s23:=[seq(seq(derivative([s2[i],s3[i]],j-1),j=1..n),i=1..n)]; s2:=[seq(i[1],i=s23)]; s3:=[seq(i[2],i=s23)]; t:=time(): M:=sympowmatrix(n,m); userinfo(3,diffop,`time taken for computing the sympower matrix`,time()-t); N:=linalg[coldim](M); # Now compute bounds for the accuracy for M[1..N,1]: bounds:=NULL; t:=time(): for i from 1 to N do bounds:=bounds, global_bounds(f,[seq(degree(M[i,1],_y[j][1]),j=1..n)], []) od; userinfo(3,diffop,`time taken for computing the bounds`,time()-t); b:=[seq(i[2]+ldegree(numer(i[1]),x)-ldegree(denom(i[1]),x) +2 -1 ,i=[bounds])]; # the +2 means that it computes 1 term too much for j from 1 to N do #lprint(j); userinfo(5,diffop,`check if non-trivial ep for col `,j): M[1,j]:=eval_semireg_monomial(expand(M[1,j]),[s1,s2,s3],b[1]) od; for j from 1 to N do if not has(M[1,j],`non-trivial ep`) then #lprint(j); userinfo(3,diffop,`evaluating column `,j): t:=time(): for i from 2 to N do if not has(M[i,1],_y[n][1]) then M[i,j]:=eval_semireg_monomial(expand(M[i,j]), [s1,s2,s3],b[i]) fi od ; userinfo(3,diffop,`time taken for column`,j,` is `,time()-t); fi; od; [evalm(M),[bounds]] end: #------------------------------------------ ../sympow_operator # Input f in C(x)[DF] # Output: the matrix of the symmetric power system #sympow_operator:=proc(f,n) # col1:=convert(linalg[col](sympowmatrix(degree(f,DF),n),1),list); # # # r:=subs(_y[1][1]=y(x),seq(_y[i][1]=diff(y(x),x$(i-1)),i=2..degree(f,DF)) # ,col1); # r:=[seq(diff(i,x),i=r)]; # DFn:=collect(DF^degree(f,DF)-f/lcoeff(f,DF),DF,normal); # r:=subs(diff(y(x),x$degree(f,DF))= # convert([seq(coeff(DFn,DF,i)*_y[i+1][1],i=0..degree(f,DF)-1)],`+`),r); # r:=expand(subs(seq(diff(y(x),x$(-i-1))=_y[-i][1],i=-degree(f,DF)..-2),y(x)=_y[1][1],r)); # linalg[matrix](nops(col1)$2,[seq(seq(coefff(i,j),j=col1),i=r)]) #end: #coefff:=proc(a,b) #if type(b,`*`) then # procname(coeff(a,op(1,b)),b/op(1,b)) #else # coeff(a,b) #fi #end: sympow_operator:=proc(f,m) local n,mi,ma,i,j,k,A,nu: with(combinat,vectoint):with(combinat,inttovec): with(linalg,matrix): n:=degree(f,DF); mi:=vectoint([m,0$(n-1)]); ma:=vectoint([0$(n-1),m]); nu:=ma-mi+1: A:=matrix(n,n,[0,1, seq(op([0$n,1]),j=1..n-2), seq(-normal(coeff(f,DF,i)/coeff(f,DF,n)),i=0..n-1)]); AA:=matrix(nu,nu,0): for i from mi to ma do vec:=inttovec(i,n): for j from 1 to n do if vec[j]<>0 then for k from 1 to n do if kj then wec:=[seq(vec[l],l=1..j-1),vec[j]-1, seq(vec[l],l=j+1..k-1),vec[k]+1,seq(vec[l],l=k+1..n)] fi: ind:=vectoint(wec)-mi+1: AA[i-mi+1,ind]:=vec[j]*A[j,k]: od: fi: od: od: evalm(AA); end: # f:=DF^3+DF+x; # n:=2; # trace(sympow_operator,coeff); #sympow_operator(f,n); apply_sympow_operator:=proc(f,n,sol) map(normal,evalm(sympow_operator(f,n) &* sol - [seq(diff(i,x),i=sol)])) end: #------------------------------------------------------------ # From the paper: if false then F:=xDF*(xDF-1/2)*(xDF-1/4)*(xDF+1/4)*(xDF-1/8)*(xDF-5/8) *(xDF+1/8)*(xDF+5/8)-x*(xDF+1/3)*(xDF-1/3); G:=collect(F,xDF); f:=convert(G,diffop,monic); fi: vi:=proc(vec) combinat[vectoint](vec) end: complete_vector:=proc(M,f,m,B) local F,n,N,i,j,k,mi,vec,aa; global a; n:=degree(f,DF): for i from 0 to n-1 do aa[i]:=normal(coeff(f,DF,i)/coeff(f,DF,n)) od: N:=linalg[coldim](M); mi:=vi([m,0$(n-1)]): ### THE ONLY THING THAT NEEDS TO BE CHANGED FOR HIGHER ORDER IS THESE NEXT BLOCK for j from 0 to m do vec:=[m-j,j,0]: index:=vi(vec)-mi+1; F[vec]:=convert([seq(`if`(has(M[1,i],`non-trivial ep`),0, a[i]*M[index,i]),i=1..N)],`+`): F[vec]:=normal(convert(series(normal(F[vec]/B[index][1]),x=0,B[index][2]+1) ,polynom)*B[index][1]) od: for j from m-1 to 0 by -1 do for i from 0 to m-j-1 do k:=m-i-j:vec:=[i,j,k]: F[vec]:= diff(F[[i,j+1,k-1]],x)-i*F[[i-1,j+2,k-1]] +(k-1)*(aa[2]*F[[i,j+1,k-1]]+aa[1]*F[[i,j+2,k-2]] +aa[0]*F[[i+1,j+1,k-2]]); F[vec]:=normal(1/(j+1)*F[vec]): od: od: ################### [seq(F[inttovec(i+mi,n)],i=0..N-1)]; end: #do_the_rest:=proc(f,NN) tester:=proc(f,NN) # option trace; n:=degree(f,DF):m:=NN: f0:=l_p(f,0); t:=time(): sols:=formal_sols(f0,[]); t1:=time(): userinfo(3,diffop,`time taken for formal sols`,t1-t); userinfo(3,diffop,`using formal solutions`,sols); M:=eval_sympowmatrix(NN,sols,f): t2:=time(); userinfo(3,diffop,`time for evaluating sympowmatrix`,t2-t1); if {op(M[2])}={[0,0]} then RETURN([0$linalg[coldim](M[1])]) fi; B:=subs([0,0]=[1,0],M[2]): # to avoid division by zero later on. userinfo(3,diffop,`bounds`,M[2]); M:=M[1]: # printlevel:=1; N:=linalg[coldim](M); eqns:={}; # The following got obsolete by using the `if` command in S:=... #### for i from 1 to N do if has(M[1,i],`non-trivial ep`) then #### eqns:=eqns union {a[i]} #### fi od; #### userinfo(3,diffop,`constants to be computed`,{seq(a[i],i=1..N)} minus eqns); ## CHANGE CUURENTLY VALID FOR ORDER 3 ONLY !! S:=complete_vector(M,f,NN,B); # S:=convert(evalm(convert([seq( # `if`(has(M[1,i],`non-trivial ep`),0,a[i]*linalg[col](M,i)) # ,i=1..N)],`+`)),list): userinfo(3,diffop,`constants to be computed`,indets(S) minus indets(f)); #### S:=subs(solve(eqns),S): S:=map(expand,S): S:=subs(log(x)=hello,S): So:=solve({seq(coeffs(coeff(i,hello,1),x),i=S)}); S:=map(normal,subs(So,S)): ## solve({seq(coeff(S[1],x,i),i=2..4)}); #solve({seq(coeff(i,x,4),i=S) , coeff(S[1],x,0)-1 }); #sol2:=map(factor,subs(",S)); # M3:=eval_sympowmatrix(3,sols,f): # Now compare S with the bounds B #SB:=[seq(convert(series(normal(S[i]/B[i][1]),x=0,B[i][2]+1),polynom),i=1..N)]: #SB:=subs(solve({seq(coeff(SB[i],x,B[i][2]+1),i=1..N)}),SB): # Si:=solve({seq(coeff(SB[i],x,B[i][2]+1),i=1..N)}): # userinfo(3,diffop,`relations on constants `, # {op(So),op(Si)}); # SB:=subs(Si,SB): # S:=[seq(normal(B[i][1]*SB[i]),i=1..N)]; # t3:=time(): # userinfo(3,diffop,`time taken for solving`,t3-t2); userinfo(3,diffop,`Now_testing.. and solving`); userinfo(3,diffop, convert(apply_sympow_operator(f,NN,S),list)); sol:=solve({seq(coeffs(expand(numer(i)),x) ,i=convert(apply_sympow_operator(f,NN,S),list))}); userinfo(3,diffop,`solutions :`,sol); userinfo(3,diffop,`total time`,time()-t); map(normal,subs(sol,S)) end: if false then deq:=x^2*(x-1)^2*diff(y(x),x$3)+(7*x-4)*x*(x-1)*diff(y(x),x$2) +(72/7*(x^2-x)-20/9*(x-1)+3/4*x)*diff(y(x),x) +(792/7^3*(x-1)+5/8+2/63)*y(x)=0; f:=convert(deq,diffop); fi: # From the file ../test_spl if false then f:=DF^3+5*(9*x^2-4*x+4)/48/x^2/(x-1)^2*DF -5/432/x^3/(x-1)^3*(81*x^3-58*x^2+102*x-44); f:=subs(x=x+1,f); f:=collect(f,DF,factor); t:=time(); result:=tester(f,6); t:=time()-t; fi: # Let r be an element of C(t), and f in C(x)[DF]. Now by writing x=r(t) # we get a differential equation in t instead of x, which will be the # output of this procedure (written in x instead of t notation). change_coordinates:=proc(f,r,t) # d/dx = d/dr = d/dt / dr/dt F:=substitute(op(subs(t=x,[DF/diff(r,t),subs(x=r,f)]))); collect(F/lcoeff(F,DF),DF,factor) end: # F_36: f:=DF^3+5*(9*x^2-4*x+4)/48/x^2/(x-1)^2*DF -5/432/x^3/(x-1)^3*(81*x^3-58*x^2+102*x-44); f2:=change_coordinates(f,t^2+1,t); do_it:=proc(NNN) options trace; for i from 0 to 1 do for j from 0 to 5 do # `diffop/back_up`(`clear all`); # For saving some memory fs:=substitute(DF+i/(2*NNN)/x+j/(6*NNN)*(2*x)/(x^2+1),f2); TT:=tester(fs,NNN); if nops(indets(TT[1]))>1 then lprint([i,j,nops(indets(TT[1]))-1]) fi od od end: # do_it(3); # result: the only one it gives is [1, 3, 1] # quit # So: that would be sqrt( x*(x^2+1) ); # Verdorie, the code doesn't yet work with algebraic numbers. # Dus zet ik er die LCLM in, maar dan gaat het me ook allemaal wat te lang # duren. do_it:=proc(NNN) options trace; for j in [3] do for i from 0 to 1 do # `diffop/back_up`(`clear all`); # For saving some memory fs:=substitute(DF+i/(2*NNN)/x+j/(6*NNN)/(x-RootOf(x^2+1)),f2); TT:=tester(LCLM(fs,`and conjugates`),NNN); if nops(indets(TT[1]))>1 then lprint([i,j,nops(indets(TT[1]))-1]) fi od od end: NNN:=3; fs:=substitute(DF+1/(2*NNN)/x+3/(6*NNN)*(2*x)/(x^2+1),f2); tester(fs,3);