Let $$$K$$$ be some field of characteristic $$$0$$$. Let $$$K[ [x_{1},...,x_{d}] ]$$$ be the ring of formal power series on it. Let $$$K[x_{1},...,x_{d}]$$$ be the polynomials of it. Let $$$K(x_{1},...,x_{d})=Frac(K[x_{1},...,x_{d}])$$$ be the field of rational functions over it.

Definition

A function $$$f \in K[ [x_{1},...,x_{d}] ]$$$ is called D-finite if for each $$$i$$$ there exists $$$n_{i}$$$ and polynomials (not all zeroes) $$$P_{i,0}(x_{1},...,x_{d}),...,P_{i,n_{i}}(x_{1},...,x_{d})$$$ such that $$$P_{i,0}f+P_{i,1}\frac{df}{dx_{i}}+...+P_{i,n_{i}}\frac{d^{n_{i}}(f)}{(dx_{i})^{n_{i}}}=0$$$.

Properties:

1 The definition is equivalent to the "linear space of derivatives of $$$f$$$ over $$$K(x_{1},...,x_{d})$$$ has finite dimension".

2 They forms a linear algebra over $$$K(x_{1},...,x_{d})$$$ (we can multiply them to scalar, add them and multiply them, and the result will be D-finite)

3 Every algebraic over $$$K(x_{1},...,x_{d})$$$ function is D-finite

4 If $$$u_{1},...,u_{d} \in K[ [u_{1},...,u_{d}] ]$$$ are algebraic, and $$$f$$$ is D-finite and composition $$$f(u_{1},...,u_{d})$$$ is well-defined in some sense, then it is D-finite.

5(Lipshitz, 1988) If $$$f$$$ is D-finite, then $$$\sum\limits_{i_{12},i_{3},...,i_{d}}^{} x_{1}^{i_{12}}x_{3}^{i_{3}}...x_{d}^{i_{d}} ([x_{1}^{i_{12}}x_{2}^{i_{12}}x_{3}^{i_{3}}...x_{d}^{i_{d}}] f)$$$ is D-finite. That is you can contract the function $$$f$$$ to the diagonal coefficients. (You can do a series of contractions to select only coefficients on diagonal) \ \textbf{6} \textbf{(Lipshitz, 1989)} \ The function $$$f=\sum\limits a(i_{1},...,i_{d}) x_{1}^{i_{1}}...x_{d}^{i_{d}}$$$ is D-finite if and only if its coefficients are P-recursive, that is for each $$$j$$$ there exists some recurence $$$\sum\limits_{v}^{} p_{v}(i_{j})a(i_{1}-v_{1},...,i_{d}-v_{d})=0$$$ where $$$v$$$ are some fixed arrays, $$$p_{v}$$$ are some fixed polynomials for each $$$v$$$. \ \texbf{7} \textbf{(Lipshitz, 1988)} \ If $$$f$$$ is D-finite, then for each $$$i_{1}$$$: $$$\sum\limits_{i_{2},i_{3},...,i_{d}}^{} x_{2}^{i_{2}}x_{3}^{i_{3}}...x_{d}^{i_{d}} ([x_{1}^{i_{1}}x_{2}^{i_{2}}x_{3}^{i_{3}}...x_{d}^{i_{d}}] f)$$$ is D-finite. Moreover, the number of polynomials and their degrees in the D-finiteness relations are uniformly bounded. \ Let's prove $$$\textbf{7}$$$ \ \begin{proof} Let's consider $$$g(x_{1},...,x_{d})=x_{1}^{-i_{1}-1}f(x_{1},...,x_{d})$$$. Then $$$f(x_{1},...,x_{d})=x^{i_{1}+1}g(x_{1},...,x_{d})$$$. So $$$g$$$ is also D-finite (Though it is not a formal power series, because it has the negative degrees of $$$x_{1}$$$, but we mean the same). \ Let's consider relations $$$P_{k}(x_{1},...,x_{d})(d/dx_{1})^{k} g + P_{k-1}(x_{1},...,x_{d})(d/dx_{1})^{k-1} g + ... +P_{0}(x_{1},...,x_{d})g = 0$$$ and $$$Q_{l}(x_{1},...,x_{d})(d/dx_{2})^{l} g + Q_{l-1}(x_{1},...,x_{d})(d/dx_{2})^{l-1} g + ... +Q_{0}(x_{1},...,x_{d})g = 0$$$. Let $$$L(x_{1},...,x_{d})=P_{k}Q_{l}$$$. \ Let $$$d=max(deg(P_{0},...,P_{k}),Q_{0},...,Q_{l})$$$. \ If we write $$$(d/dx_{1})^{u}(d/dx_{2}^{v})(P_{k}(x_{1},...,x_{d})(d/dx_{1})^{k} g + P_{k-1}(x_{1},...,x_{d})(d/dx_{1})^{k-1} g + ... +P_{0}(x_{1},...,x_{d})g)=0$$$ (we are differentiating first relation many times) then we get $$$ P_{k} (d/dx_{1})^{u+k}(d/dx_{2})^{v} g + ... $$$(lower differentials)$$$=0$$$ with the degree of all polynomials here being at most $$$d$$$. Same with second relation. (Here we use the fact that after differentiating the degree of polynomial decreases) \ So if we have $$$L(d/dx_{1})^{s}(d/dx_{2})^{t} g$$$ we can simplify it to the lower differentials getting $$$L(d/dx_{1})^{s}(d/dx_{2})^{t} g = \sum_{s'+t' \lt s+t} p_{s',t'}(x_{1},...,x_{d})((d/dx_{1})^{s'}(d/dx_{2})^{t'} g)$$$ where $$$deg(p_{s',t'}) \leq d$$$. \ Then let's write $$$L^{N}x_{2}^{\alpha_{2}}x_{3}^{\alpha_{3}}...x_{d}^{\alpha_{d}}(d/dx_{1})^{m_{1}}(d/dx_{2})^{m_{2}} g$$$, where $$$N$$$ is fixed, $$$\alpha_{2}+\alpha_{3}+...+\alpha_{d}+m_{1}+m_{2} \leq N$$$. Using the above things we can simplify such things to a form of $$$L^{N-m_{1}-m_{2}}x_{2}^{\alpha_{2}}x_{3}^{\alpha_{3}}...x_{d}^{\alpha_{d}}\sum\limits_{s,t}^{} P_{s,t}(x_{1},...,x_{d})((d/dx_{i})^{s}(d/dx_{2})^{t}) g$$$ where $$$s \lt k,t \lt l$$$ and $$$P_{s,t}(x_{1},...,x_{d})$$$ are polynomials of degree at most $$$2(m_{1}+m_{2})d$$$ (Each time we are simplifying the highest by $$$s+t$$$ differentials using one $$$L$$$). \ So $$$L^{N}x_{2}^{\alpha_{2}}x_{3}^{\alpha_{3}}...x_{d}^{\alpha_{d}}(d/dx_{1})^{m_{1}}(d/dx_{2})^{m_{2}} g$$$ are the linear combination with coefficients in $$$K$$$ of $$$sum\limits_{s,t}^{} p_{s,t}(x_{1},...,x_{d})((d/dx_{i})^{s}(d/dx_{2})^{t}) g$$$ where $$$deg(p_{s,t}) \leq 2Nd$$$. \ So the dimension of the linear space over $$$K$$$ generated by them is at most $$$kl\binom{2Nd+d}{d}$$$. \ But there are $$$\binom{N+d+1}{d+1}$$$ of them, for big $$$N$$$: $$$\binom{N+d+1}{d+1} \geq kl\binom{2Nd+d}{d}$$$, so they are linear dependent! \ So for this big $$$N$$$ there is a relation $$$\sum c_{\alpha_{i},m_{1},m_{2}}L^{N}x_{2}^{\alpha_{2}}x_{3}^{\alpha_{3}}...x_{d}^{\alpha_{d}}(d/dx_{1})^{m_{1}}(d/dx_{2})^{m_{2}} g=0$$$. Let's divide it to $$$L^{N}$$$ and write $$$\sum c_{\alpha_{i},m_{1},m_{2}}x_{2}^{\alpha_{2}}x_{3}^{\alpha_{3}}...x_{d}^{\alpha_{d}}(d/dx_{1})^{m_{1}}(d/dx_{2})^{m_{2}} g=0$$$.\ Then let's explain why we actually divide $$$f$$$ to $$$x_{1}^{i_{1}+1}$$$ in the first step. \ We have $$$\sum {c_{\alpha_{i},m_{1},m_{2}}x_{2}^{\alpha_{2}}x_{3}^{\alpha_{3}}...x_{d}^{\alpha_{d}}(d/dx_{1})^{m_{1}}(d/dx_{2})^{m_{2}} g}=0$$$ (not all $$$c$$$ are zeroes), let's $$$m$$$ be the minimal $$$m_{1}$$$ for which there is $$$c_{...} \neq 0$$$ (not all coefficients with this $$$m_{1}$$$ are zeroes). Then let's consider the terms in the $$$\sum c_{\alpha_{i},m_{1},m_{2}}x_{2}^{\alpha_{2}}x_{3}^{\alpha_{3}}...x_{d}^{\alpha_{d}}(d/dx_{1})^{m_{1}}(d/dx_{2})^{m_{2}} g=0$$$ with $$$x_{1}^{-m-1}$$$. And they are coming to us only from differentiating coefficients with $$$x_{1}^{-1}$$$ $$$m$$$ times, because all coefficients with non-negative degree of $$$x_{1}$$$ can't differentiate to negative degrees (when being $$$x_{1}^{0}$$$ it differentiates (by $$$x_{1}$$$) to $$$0$$$). And that's why we divided everything to $$$x_{1}^{i_{1}+1}$$$ in the first step. \ So consider the coefficient with $$$x_{1}^{-m-1}$$$, then we would have here (after dividing to $$$m_{1}!(-1)^{m_{1}}$$$) $$$\sum {c_{\alpha_{i},m,m_{2}}x_{2}^{\alpha_{2}}x_{3}^{\alpha_{3}}...x_{d}^{\alpha_{d}}(d/dx_{2})^{m_{2}} ([x_{1}^{-1}]g)}=0$$$. So we got a recurence relation for $$$(d/dx_{2})^{k} ([x_{1}^{-1}]g)$$$. Doing this thing with not only $$$x_{2}$$$, but $$$x_{3},...,x_{d}$$$ too we get that $$$([x_{1}^{-1}]g)$$$ is D-finite. So $$$([x_{1}^{i_{1}}] f)$$$ is D-finite. (with the uniformly bounded length of relations and the degrees of polynomials). \ \end{proof} Note: Here in the proof we actually have shown that the non-negative (by degree of $$$x_{1}$$$) part of a D-finite function is D-finite. \ \ \ \ \textbf{What functions are D-finite?} \ Can you guess, what of this functions are D-finite? \ (1) $$$a_{n}=(1+x+x^{2})^{n}[x^{n+1}]$$$ \ (2) $$$a_{n}=n!$$$ \ (3) $$$a_{n,k}=S(n,k)$$$ — Stirling numbers of the second kind, number of ways to partition a set of $$$n$$$ objects into $$$k$$$ non-empty subsets. \ $$$\sum\limits_{n=0}^{\infty} \sum\limits_{k=0}^{n} S(n,k) \frac{x^n}{n!} y^k = e^{y(e^x-1)}$$$ \ (4) Eulerian numbers $$$a_{n,k}=E(n,k)$$$ is a number of $$$n$$$-element permutations with $$$k$$$ ascents. \ $$$\sum\limits_{n=0}^{\infty} \sum\limits_{k=0}^{n} E(n,k)t^{k} \frac{x^n}{n!} = \frac{t-1}{t - e^{(t-1)\,x}} = \left(1-\frac{e^{(t-1)x}-1}{t-1}\right)^{-1}$$$ \ (5) $a_{u,v,w}=\sum\limits_{0 \leq x \leq u, 0\leq y \leq v, 0\leq z \leq w} \frac{(x+y)!(y+z)!(z+x)!}{x!y!z!}$. \ Answer: $$$(1),(2),(5)$$$ are D-finite, $$$(3),(4)$$$ I suppose to be not D-finite (though I don't have proof). \ For $$$(1)$$$ it is a shift by $$$x^{-1}$$$ + contraction to the diagonal of the D-finite function $$$\frac{1}{1-y(1+x+x^2)}$$$ \ For $$$(2)$$$ it is clear \ For $$$(3)$$$ the problem, why it can be not D-finite is that we get the exponent of non-algebraic thing. \ For $$$(4)$$$ the problem, why it can be not D-finite is that when we take the last inverse of the whole thing, this thing is not algebraic. \ $$$(5)$$$ is D-finite because it is contraction to the diagonal of some product of factorials and inverse factorials, then multiplied by $$$\frac{1}{(1-x)(1-y)(1-z)}$$$. \ \ \ \ The D-finite functions with $$$d=1$$$ are just P-recursive sequences. And we can evaluate first n terms of a P-recursive sequence in a stupid way in $$$O(n)$$$, or there is an algorithm to evaluate the $$$n$$$-th term (only one) in $$$O(\sqrt{n}log(n))$$$. \ Algorithm: \ The algorithm is pretty similar to the one for factorials. \ Let $$$l$$$ is the length of the P-recursive relation ($$$P_{l}(n)a_{n-l}+...+P_{0}(n)a_{n}=0$$$). And let $$$d=max(deg(P_{i}))$$$. We know $$$a_{0},...,a_{l-1}$$$, we want to compute $$$a_{n}$$$. ($$$l$$$ and $$$d$$$ are assuming to be constants) \ Let $$$B=2^{u} \geq \sqrt{n}$$$, $$$B$$$ is a minimal such. Let's divide everything to a blocks of size $$$B$$$. Then let $$$k=\lfloor \frac{n}{B} \rfloor$$$, let's compute $$$a_{kB},...,a_{kB+l-1}$$$, then we can easily compute $$$a_{n}$$$ in $$$O(\sqrt{n})$$$. \ Let we have some $$$a_{x},...,a_{x+l-1}$$$ and we want to compute $$$a_{x+diff},...,a_{x+diff+l-1}$$$, then we will have \ $\begin{pmatrix} a_{x+diff}\ a_{x+diff+1}\ . \ . \ . \ a_{x+diff+l-1} \end{pmatrix} = \frac{1}{S(x)}\begin{pmatrix} P_{00} (x) & P_{01}(x) ... & P_{0(l-1)}(x)\ P_{10} (x) & P_{11}(x) ... & P_{1(l-1)}(x)\ . \ . \ . \ P_{(l-1)0} (x) & P_{(l-1)1}(x) ... & P_{(l-1)(l-1)}(x) \end{pmatrix} \begin{pmatrix} a_{x}\ a_{x+1}\ . \ . \ . \ a_{x+l-1} \end{pmatrix}$. \ \ \ Where $$$S(x)=P_{0}(x+l)...P_{0}(x+diff+l-1)$$$ and $$$max(deg(S),deg(P_{i,j})) \leq diff \cdot d$$$. \ We will denote the polynomials for a given $$$diff$$$ as $$$P_{diff,i,j}$$$ and $$$S_{diff}$$$. \ Let's store the polynomial of degree $$$d$$$ as an array of values in a points $$$0,K,2K,...,K(d-1),Kd$$$ for some $$$K$$$. For a given $$$K$$$ we can multiply polynomials in such form (to get such form of a product) in $$$O((d_{1}+d_{2})log(d_{1}+d_{2}))$$$ by doing Taylor shifts. \ Let $$$P_{d,i,j}$$$ and $$$S_{d}$$$ will be such polynomials for given $$$d$$$. \ Then for some $$$s$$$ let's store $$$P_{2^s,i,j}(0),P_{2^s,i,j}(2^s),...,P_{2^s,i,j}(2^{2s} \cdot d)$$$. And same with $$$S$$$, that is $$$S_{2^s}(0),...,S_{2^s,i,j}(2^{2s} \cdot d)$$$. \ Then from $$$s=0$$$ to $$$s=u$$$ we can compute all of this by shifts of sampling points (this is the operation that requires FFT) and matrix multiplications. With $$$s=u$$$ we are done. \ The final complexity is $$$O(\sqrt{n}log(n)l^{2}d+\sqrt{n}l^{3}d)$$$, but as we assuming $$$l$$$ and $$$d$$$ constants, then it is just $$$O(\sqrt{n}log(n))$$$. \ \ \ \

But how to solve problems using this? \ Let us assume that some function $$$f$$$ is D-finite, for example $$$f=\frac{P(x,y)}{Q(x,y)}$$$ where we know $$$P(x,y)$$$ and $$$Q(x,y)$$$ are fixed polynomials. And we want to compute $$$[x^{n}y^{m}] f$$$ (number of variables actually does not matter). Let's notice that $$$h(u)=[x^{n}y^{u}] f$$$ is D-finite (this is the property \textbf{7}) that is equivalent to P-recursive, and we can compute $$$h(m)$$$ in $$$O(\sqrt{m}log(m))$$$ from some first values of $$$h$$$: $$$h(0),...,h(k)$$$ where $$$k$$$ is a constant. \ To compute $$$h(l)$$$ we will do the same, we will notice that $$$g(u)=[x^{u}y^{l}] f$$$ is P-recursive, so we can compute $$$g(l)$$$ in $$$O(\sqrt{n}log(n))$$$. \ So we can compute $$$[x^{n}y^{m}] f$$$ in $$$O(\sqrt{n+m}log(n+m))$$$. \ But when I saying let's compute the coefficient of this P-recursive sequence, how do I know the P-recursive relation? \ One can compute it from the proof of property $$$\textbf{7}$$$, but I think this is a lot of pain. So to find P-recursive relation let's compute the some first $$$k$$$ coefficients and run a Gauss. \ Note: We can also notice from the proof of property $$$\textbf{7}$$$, that the coefficients of P-recursion of $$$h(u)$$$ can be expressed as a polynomials of $$$u$$$, and then we can try to evaluate them from interpolation of a projective function, though I have not implemented it yet. \ \ \ \ Example of this algorithm: \ Let we have a function $$$f(n,m)=[x^ny^m] e^{\frac{x+y}{(1-x)(1-y)}}$$$. And we want to compute $$$f(10^{7},10^{8})$$$ (mod 998244353) \ This function $$$f$$$ is D-finite. \ Let's compute it values for $$$0 \leq n \lt 25$$$ and $$$0 \leq m \lt 25$$$. Then we can extend $$$f(\cdot,0),f(\cdot,1),f(\cdot,2),...,f(\cdot,24)$$$ to $$$f(10^{7},0),f(10^{7},1),...,f(10^{7},24)$$$ by P-recursion. Then we can extend $$$f(10^{7},\cdot)$$$ from first $$$25$$$ values to $$$f(10^{7},10^{8})$$$ by P-recursion too. It works in $$$O(25\sqrt{10^{7}}log(10^{7})+\sqrt{10^{8}}log(10^{8}))$$$ which takes ~3.1 sec. \ \ \ \ Implementation: \ Here is a code (paste!) \ Code works for all prime modulos \ findPreccursion finds P-recursion \ evaluatePrecursion evaluates P-recursion from the initial values if possible \ evaluatePrecursionfast evaluates the element of sequence from P-recursion and the initial values \ sequenceextender is trying to extend the sequence \ optimalgetvaluebyid is trying to get the element of an extended sequence by its position \ extendtable extends table A, finding A(que[i].first,que[i].second) \ getcolumnoftable/getrowoftable get column and row of an extended table \
\ \ \ But there is one caveat in this algorithm. Assume that there is a sequence $$$a_{i}$$$ with P-recursion with relation $$$a_{n}P_{0}(n)+...+a_{n-l}P_{l}(n)=0$$$. You know $$$a(0),...,a(l-1)$$$ and you want to extend this sequence to $$$a(0),...,a(N)$$$. But what if $$$P_{0}(k)=0$$$ for some $$$l \leq n \leq N$$$. Then we can't extend this sequence. And I don't know what to do with that, but luckily, I suppose that for almost all function when we can use it, it can't be zero. \ \ \ Example where it matters: \ Let's try to $$$f(10^{6},10^{6})$$$ with $$$f(n,m)=[x^{n}y^{m}] \frac{1}{1-xy}$$$. We will try to calculate some first columns, and then a row. But even if we managed to somehow guess the first columns, we would be unable to get a row, because our initial values in this row will all be zeroes, and we want to get $$$1$$$ somehow, but we can't. \ \ \ Example of the problems: \
\hyperlink{https://atcoder.jp/contests/arc202/tasks/arc202_d}{Atcoder Regular 202 problem D} \ Here if you are solving it in $$$O((H+w+T)^{1.5}$$$ using the reflection principle you have a subproblem calculate the $$$(1+x+x^{2})^{a}[x^{a+k}]$$$ for all $$$0 \leq a \leq T$$$. This can be solved by FFT (and even in such a way, that it will not add the logarithm to the asymptotic of the solution), but this also can be solved using the fact that it is P-recursive (because it is а column shift of a diagonal) and find coefficients with $$$a=k,...,k+C$$$, and then extend it to the whole sequence in a linear time. \ \ We can also solve \hyperlink{https://codeforces.me/contest/1747/problem/E}{1747E} in $$$O(\sqrt{n}log(n))$$$ per test (compute the coefficient of $$$(1-x)^{2}(1-y)^{2}(1-2x-2y+2xy)^{-2}$$$), but actually, because there is a multitest and the constraints are big this solution is hard to pass (I have managed, but after a lot of pain and with some strange FFT \hyperlink{https://codeforces.me/contest/1747/submission/330486092}{submission}) \ \ I think there are some other examples of problems that can be solved by this method. \ \end{document}