############################################################################## ### INPUT: a trivariate rational function f in s,t,x. ### OUTPUT: a linear differential operator L(x,Dx) annihilating Diag(f). ### NOTES on the algorithm: it is based on Remark 3 in [Lipshitz89]: let ### F(s,t,x) = 1/s/t * f(s,t/s,x/t). If P(x;Dx,Ds,Dt) annihilates F, then ### P(x;Dx,0,0) kills Diag(f). ### The idea is to determine such a P by undetermined coefficients: what makes ### things work is the observation that P(F)=0 amounts to a linear homogeneous ### system in the undetermined coefficients of P, having nb_eqs = O(N^3) ### equations and nb_unk = O(N^4) unknowns. ### In our approach (which is more in Fasenmyer's style), the notable ### difference is that the equation P(F)=0 is also equivalent to a ### homogeneous linear system in the (undetermined) coefficients c_{ijk}(x) of ### P = sum(c_{ijk}*Dx^i*Ds^j*Dt^k). The system has nb_eqs = O(N^2) equations ### and nb_unk = O(N^3) unknowns, these unknowns being polynomials in x, ### instead of scalars. ### The algorithm works by increasing the value of N until nb_unk > nb_eqs. ### This value of N is typically smaller than the one in Lipshitz's approach. ### On our particular example of 3D Rooks, we have: ### nb_eqs:=(13*N+14)*(N+1)/2 ### nb_unk:=binomial(N+3,3); ############################################################################## Fasenmyer3rat:=proc(f, delay := 0, \$) local F, Q, L, N, NB_eqs, NB_unk, del, NumDer, j, i, k, eqs, sys, unk, nb_eqs, nb_unk, x0, sol0, sol, Eq, Lstu, pol_eqs, pol_unk, tmp, Qx, Qs, Qt, zob; F:=normal(1/s/t*subs({x=s,s=t/s,t=x/t},f)); Q:=denom(F); Qx := diff(Q, x) ; Qs := diff(Q, s) ; Qt := diff(Q, t) ; L:=0; NB_eqs:=[]: NB_unk:=NB_eqs: del := delay ; NumDer[0, 0, 0] := numer(F) ; # Invariant: NumDer[i, j, k] / Q^(N+1) = diff(F, [x\$i, s\$j, t\$k]). for N while L = 0 or del > 0 do if L <> 0 then del := del - 1 end if ; print(N); print("setting up the linear system"); print("derivatives"); tmp := NumDer[0, N-1, 0] ; NumDer[0, N, 0] := expand( diff(tmp, s) * Q - N * tmp * Qs ) ; tmp := NumDer[0, 0, N-1] ; NumDer[0, 0, N] := expand( diff(tmp, t) * Q - N * tmp * Qt ) ; for i to N do tmp := NumDer[0, i-1, N-i] ; NumDer[0, i, N-i]:=expand( diff(tmp, s) * Q - N * tmp * Qs ) od ; for i from 1 to N do for j from 0 to N-i do tmp := NumDer[i-1, j, N-i-j] ; NumDer[i, j, N-i-j] := expand( diff(tmp, x) * Q - N * tmp * Qx ) od od ; NumDer[0, 0, 0] := normal(NumDer[0, 0, 0] * Q) ; for j from 0 to N-2 do NumDer[0, j+1, 0] := normal(NumDer[0, j+1, 0] * Q) ; NumDer[0, 0, j+1] := normal(NumDer[0, 0, j+1] * Q) ; od ; for i to N do for j to N-i-1 do NumDer[0, i, j] := normal(NumDer[0, i, j] * Q) od od ; for i from 1 to N do for j from 0 to N-i do for k from 0 to N-i-j-1 do NumDer[i, j, k] := normal(NumDer[i, j, k] * Q) od od od ; if L <> 0 and del > 0 then next end if ; print("add equations"); eqs:=normal(add(add(add(c[i, j, k]*NumDer[i, j, k], k=0..N-j-i),j=0..N-i),i=0..N)): print("extract equations"); eqs := expand(eqs) ; sys := [coeffs(eqs, {s, t})] ; unk := subs(x = NULL, indets(sys)) ; nb_eqs:=nops(sys); NB_eqs:=[op(NB_eqs),nb_eqs]: nb_unk:=nops(unk); NB_unk:=[op(NB_unk),nb_unk]: print("number of equations and unkowns", nb_eqs, nb_unk); print("degree in x", max(seq(degree(eq,x), eq=eqs))); print("solving the linear system"); x0:=nextprime(rand(2^10)(1)); print("at a random point x =", x0); sol0 := SolveTools:-Linear(eval(sys, x = x0), unk) ; if map2(op, 2, sol0) <> {0} and L * del = 0 then print("it seems that there is a solution"); sol := SolveTools:-Linear(sys, unk) ; print("computing annihilating operator"); Eq:=0; for i from 0 to N do for j from 0 to N-i do for k from 0 to N-i-j do Eq:= normal(Eq+c[i, j, k]*Dx^i*Ds^j*Dt^k); od od od; Lstu:=subs(sol,numer(Eq)); if Lstu<>0 then L:=subs({Ds=0,Dt=0},Lstu); end if; end if; if N = 6 then pol_eqs:=factor(CurveFitting[PolynomialInterpolation]([\$1..nops(NB_eqs)], NB_eqs,Z)); pol_unk:=factor(CurveFitting[PolynomialInterpolation]([\$1..nops(NB_unk)], NB_unk,Z)); print("number of equations",pol_eqs); print("number of unknowns", pol_unk); end if od ; return collect(numer(L),Dx), numer(Lstu), pol_eqs, pol_unk; end proc; #### MinimalFasenmyer3rat := proc(f, delay := 0, \$) local tt, FR, L, ind, L1, L2, G; tt:=time(): FR:=Fasenmyer3rat(f, delay); print("time Fasenmyer3rat", time()-tt); L:=FR[1]; ind := indets(L) minus {x, Dx} ; L1:=subs([seq(s=rand(100)(1),s=ind)],L): L2:=subs([seq(s=rand(100)(1),s=ind)],L): tt:=time(): G:=DEtools[GCRD](L1,L2,[Dx,x]); print("time GCRD", time()-tt); return G; end proc; #### simplest example #f:=1/(1-x-s-t): #MinimalFasenmyer3rat(f); #### 3D rook f:=1/(1 - s/(1-s) - t/(1-t) - x/(1-x) ); MinimalFasenmyer3rat(f);