Introduction↵
----------↵
↵
7 months ago [user:Elegia,2025-07-25] published a blog about D-finite functions ([part1](https://codeforces.me/blog/entry/138078), [part2](https://codeforces.me/blog/entry/138113)) in which he showed some examples of problems that can be solved by D-finite functions.↵
↵
Here I will explain a fast evaluation algorithm of the coefficients of a D-finite function.↵
↵
Let $K$ be some field of characteristic $0$. Let $K[ [x_{1},\ldots,x_{d}] ]$ be the ring of formal power series on it. Let $K[x_{1},\ldots,x_{d}]$ be the polynomials of it. Let $K(x_{1},\ldots,x_{d})=\mathop{\mathrm{Frac}}(K[x_{1},\ldots,x_{d}])$ be the field of rational functions over it. Let $ [ x_{1}^{i_{1}}\ldots x_{d}^{i_{d}} ] f$ be the coefficient in $f$ of the monomial $x_{1}^{i_{1}}\ldots x_{d}^{i_{d}}$.↵
↵
Definitions↵
----------↵
↵
A function $f \in K[ [x_{1},\ldots,x_{d}] ]$ is called _D-finite_ if for each $i$ there exists a positive integer $n_{i}$ and polynomials (not all zeroes) $P_{i,0}(x_{1}, \ldots,x_{d}),\ldots,P_{i,n_{i}}(x_{1},\ldots,x_{d})$ such that $P_{i,0}f+P_{i,1}\frac{df}{dx_{i}}+\ldots+P_{i,n_{i}}\frac{d^{n_{i}}(f)}{(dx_{i})^{n_{i}}}=0$.↵
↵
An infinite table of coefficients $a(i_{1},\ldots,i_{d})$ is called _P-recursive_ if there exists $k \in \mathbb{N}$ such that for each $j \in \\{ 1,\ldots,d \\} $ there exists a mapping $p: \\{ 0,\ldots,k \\} ^{d} \rightarrow K[x]$ such that at least one $p_{v} \neq 0$ and for all $i_{1},\ldots,i_{d} \geq k$ $\sum\limits_{v \in \\{0,\ldots,k\\}^{d}}^{} p_{v}(i_{j}) a(i_{1}-v_{1},\ldots,i_{d}-v_{d}) = 0$. ↵
↵
In particular, for $d=1$ it means that there exist such polynomials (not all zeroes) $P_{i}$: $ \forall n \geq k : \sum\limits_{0 \leq i \leq k}^{} P_{i}(n)a(n-i)=0$.↵
↵
Properties:↵
-----------↵
↵
**1** The definition is equivalent to the "linear space of partial derivatives of $f$ over $K(x_{1},\ldots,x_{d})$ has finite dimension". ↵
↵
In particular, every rational function is obviously D-finite↵
↵
**2** They form a linear algebra over $K(x_{1},\ldots,x_{d})$ (we can multiply them by a scalar, add them and multiply them, and the result will be D-finite)↵
↵
**3** Every algebraic over $K(x_{1},\ldots,x_{d})$ function is D-finite ↵
↵
**4** If $u_{1},\ldots,u_{d} \in K[ [u_{1},\ldots,u_{d}] ]$ are algebraic, and $f$ is D-finite and composition $f(u_{1},\ldots,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},\ldots,i_{d}}^{} x_{1}^{i_{12}}x_{3}^{i_{3}}\ldots x_{d}^{i_{d}} ([x_{1}^{i_{12}}x_{2}^{i_{12}}x_{3}^{i_{3}}\ldots 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 the diagonal) ↵
↵
**6(Lipshitz, 1989)** The function $f=\sum\limits a(i_{1},\ldots,i_{d}) x_{1}^{i_{1}}\ldots x_{d}^{i_{d}}$ is D-finite if and only if its coefficients are P-recursive.↵
↵
**7(Lipshitz, 1988)** If $f$ is D-finite, then for each $i_{1}$: $\sum\limits_{i_{2},i_{3},\ldots,i_{d}}^{} x_{2}^{i_{2}}x_{3}^{i_{3}}\ldots x_{d}^{i_{d}} ([x_{1}^{i_{1}}x_{2}^{i_{2}}x_{3}^{i_{3}}\ldots x_{d}^{i_{d}}] f)$ is D-finite. Moreover, the number of polynomials and their degrees in the D-finiteness relations are uniformly bounded. ↵
↵
Note: D-finite functions are not a field ↵
↵
Let's prove **7**↵
↵
<spoiler summary="Proof">↵
↵
Let's consider $g(x_{1},\ldots,x_{d})=x_{1}^{-i_{1}-1}f(x_{1},\ldots,x_{d})$. Then $f(x_{1},\ldots,x_{d})=x^{i_{1}+1}g(x_{1},\ldots,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},\ldots,x_{d})(d/dx_{1})^{k} g + P_{k-1}(x_{1},\ldots,x_{d})(d/dx_{1})^{k-1} g + \ldots +P_{0}(x_{1},\ldots,x_{d})g = 0$ and $Q_{l}(x_{1},\ldots,x_{d})(d/dx_{2})^{l} g + Q_{l-1}(x_{1},\ldots,x_{d})(d/dx_{2})^{l-1} g + \ldots +Q_{0}(x_{1},\ldots,x_{d})g = 0$. Let $L(x_{1},\ldots,x_{d})=P_{k}Q_{l}$. ↵
↵
Let $d=\max(\deg(P_{0},\ldots,P_{k}),Q_{0},\ldots,Q_{l})$. ↵
↵
If we write $(d/dx_{1})^{u}(d/dx_{2}^{v})(P_{k}(x_{1},\ldots,x_{d})(d/dx_{1})^{k} g + P_{k-1}(x_{1},\ldots,x_{d})(d/dx_{1})^{k-1} g + \ldots +P_{0}(x_{1},\ldots,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 + \ldots \text{(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'<s+t} p_{s',t'}(x_{1},\ldots,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}}\ldots x_{d}^{\alpha_{d}}(d/dx_{1})^{m_{1}}(d/dx_{2})^{m_{2}} g$, where $N$ is fixed, $\alpha_{2}+\alpha_{3}+\ldots+\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}}\ldots x_{d}^{\alpha_{d}}\sum\limits_{s,t}^{} P_{s,t}(x_{1},\ldots,x_{d})((d/dx_{i})^{s}(d/dx_{2})^{t}) g$ where $s<k,t<l$ and $P_{s,t}(x_{1},\ldots,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}}\ldots 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},\ldots,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}}\ldots 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}}\ldots 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}}\ldots x_{d}^{\alpha_{d}}(d/dx_{1})^{m_{1}}(d/dx_{2})^{m_{2}} g}=0$ (not all $c$ are zeroes), let $m$ be the minimal $m_{1}$ for which there is $c_{\ldots} \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}}\ldots 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 is 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}}\ldots 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},\ldots,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).↵
↵
</spoiler>↵
↵
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. ↵
↵
What functions are D-finite?↵
----------------------------↵
↵
Can you guess which of these functions are D-finite (that is, the sequence/table of coefficients is P-recursive)? ↵
↵
(1) $a_{n} = \[ x^{n+1} \] (1+x+x^{2})^{n}$ ($f(x)=\sum\limits_{}^{} a_{n}x^{n}$)↵
↵
(2) $a_{n}=n!$ ($f(x)=\sum\limits_{}^{} n!x^{n}$)↵
↵
(3) $a_{n,k}=S(n,k)$ [Stirling numbers of the second kind](https://en.wikipedia.org/wiki/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)}$ ↵
↵
(It is not so hard to see, that the D-finiteness of the OGF of some sequence/table is equivalent to the D-finiteness of EGF of this sequence/table)↵
↵
(4) [Eulerian numbers](https://en.wikipedia.org/wiki/Eulerian_number) $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!}$. ↵
↵
<spoiler summary="Answer">↵
↵
$(1),(2),(5)$ are D-finite, $(3),(4)$ I believe are not D-finite. ↵
↵
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, the reason it may not be D-finite is that we get the exponent of non-algebraic thing (and we can prove it is not D-finite). ↵
↵
For $(4)$ the problem, the reason it may not be D-finite is that when we take the last inverse of the whole thing, this thing is not algebraic (and we can prove it is not D-finite). ↵
↵
$(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)}$. ↵
↵
</spoiler>↵
↵
The D-finite functions with $d=1$ are just P-recursive sequences. And we can evaluate the 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))$. ↵
↵
<spoiler summary="Algorithm"> ↵
↵
The algorithm is pretty similar to the one for factorials. ↵
↵
Let $l$ be the length of the P-recursive relation ($P_{l}(n)a_{n-l}+\ldots+P_{0}(n)a_{n}=0$). And let $d=\max(\deg(P_{i}))$. We know $a_{0},\ldots,a_{l-1}$, we want to compute $a_{n}$. ($l$ and $d$ are assumed 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},\ldots,a_{kB+l-1}$, then we can easily compute $a_{n}$ in $O(\sqrt{n})$. ↵
↵
Suppose we have some $a_{x},\ldots,a_{x+l-1}$ and we want to compute $a_{x+\mathrm{diff}},\ldots,a_{x+\mathrm{diff}+l-1}$, then we will have↵
↵
↵
↵
$$\begin{pmatrix} a_{x+\mathrm{diff}}\\ a_{x+\mathrm{diff}+1}\\ . \\ . \\ . \\ a_{x+\mathrm{diff}+l-1} \end{pmatrix} = \frac{1}{S(x)}\begin{pmatrix} P_{00} (x) & P_{01}(x) \ldots & P_{0(l-1)}(x)\\ P_{10} (x) & P_{11}(x) \ldots & P_{1(l-1)}(x)\\ . \\ . \\ . \\ P_{(l-1)0} (x) & P_{(l-1)1}(x) \ldots & 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)\ldots P_{0}(x+\mathrm{diff}+l-1)$ and $\max(\deg(S),\deg(P_{i,j})) \leq \mathrm{diff} \cdot d$. ↵
↵
We will denote the polynomials for a given $\mathrm{diff}$ as $P_{\mathrm{diff},i,j}$ and $S_{\mathrm{diff}}$. ↵
↵
Let's store the polynomial of degree $d$ as an array of values in a points $0,K,2K,\ldots,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),\ldots,P_{2^s,i,j}(2^{2s} \cdot d)$. And same with $S$, that is $S_{2^s}(0),\ldots,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 are assuming $l$ and $d$ to be constants, then it is just $O(\sqrt{n}\log(n))$. ↵
↵
</spoiler>↵
↵
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 $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 **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),\ldots,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 say "Let's compute the coefficient of this P-recursive sequence", how do I know that the P-recursive relation holds? ↵
↵
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 first $k$ coefficients (for some $k$) 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 polynomials of $u$, and then we can try to evaluate them from some interpolation of a rational function, though I have not implemented it yet. ↵
↵
Example of this algorithm: ↵
↵
Assume 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}) \pmod{998\,244\,353}$. ↵
↵
This function $f$ is D-finite. ↵
↵
Let's compute its values for $0 \leq n <25$ and $0 \leq m < 25$. Then we can extend ↵
↵
$f(\cdot,0),f(\cdot,1),f(\cdot,2),\ldots,f(\cdot,24)$ to $f(10^{7},0),f(10^{7},1),\ldots,f(10^{7},24)$ by P-recursion. Then we can extend $f(10^{7},\cdot)$ from the 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. ↵
↵
Note: for $f(10^{8},10^{8})$ we can compute it from diagonal (property **5**). ↵
↵
Implementation:↵
---------------↵
↵
Here is a code↵
↵
<spoiler summary="Code">↵
```cpp↵
#include <bits/stdc++.h>↵
↵
using namespace std;↵
//#define int long long↵
///You can add it if you want↵
const int md = 998244353;↵
mt19937 rnd;↵
#define app push_back↵
#define all(x) (x).begin(),(x).end()↵
#ifdef LOCAL↵
#define debug(...) [](auto...a){ ((cout << a << ' '), ...) << endl;}(#__VA_ARGS__, ":", __VA_ARGS__)↵
#define debugv(v) do {cout<< #v <<" : {"; for(int izxc=0;izxc<v.size();++izxc) {cout << v[izxc];if(izxc+1!=v.size()) cout << ","; }cout <<"}"<< endl;} while(0)↵
#else↵
#define debug(...)↵
#define debugv(v)↵
#endif↵
#define lob(a,x) lower_bound(all(a),x)↵
#define upb(a,x) upper_bound(all(a),x)↵
↵
template<int M, int K, int G>↵
struct Fft↵
{↵
// 1, 1/4, 1/8, 3/8, 1/16, 5/16, 3/16, 7/16, ...↵
int g[1 << (K - 1)];↵
↵
Fft() : g()↵
{ //if tl constexpr...↵
// static_assert(K >= 2, "Fft: K >= 2 must hold");↵
g[0] = 1;↵
g[1 << (K - 2)] = G;↵
for (int l = 1 << (K - 2); l >= 2; l >>= 1)↵
{↵
g[l >> 1] = (g[l] * 1LL * g[l]) % M;↵
}↵
assert((g[1]*1LL * g[1]) % M == M - 1);↵
for (int l = 2; l <= 1 << (K - 2); l <<= 1)↵
{↵
for (int i = 1; i < l; ++i)↵
{↵
g[l + i] = (g[l] * 1LL * g[i]) % M;↵
}↵
}↵
}↵
↵
void fft(vector<int> &x) const↵
{↵
const int n = x.size();↵
assert(n <= 1 << K);↵
for (int h = __builtin_ctz(n); h--;)↵
{↵
const int l = (1 << h);↵
for (int i = 0; i < n >> (h + 1); ++i)↵
{↵
for (int j = i << (h + 1); j < (((i << 1) + 1) << h); ++j)↵
{↵
const int t = (g[i] * 1LL * x[j | l]) % M;↵
x[j | l] = x[j] - t;↵
if (x[j | l] < 0)↵
x[j | l] += M;↵
x[j] += t;↵
if (x[j] >= M)↵
x[j] -= M;↵
}↵
}↵
}↵
for (int i = 0, j = 0; i < n; ++i)↵
{↵
if (i < j)↵
std::swap(x[i], x[j]);↵
for (int l = n; (l >>= 1) && !((j ^= l) & l);)↵
{↵
}↵
}↵
}↵
↵
vector<int> convolution(vector<int> a, vector<int> b) const↵
{↵
if (a.empty() || b.empty())↵
return {};↵
const int p =mdM;↵
for (int &x: a)↵
{↵
x %= p;↵
if (x >= p)↵
x -= p;↵
if (x < 0)↵
x += p;↵
}↵
for (int &x: b)↵
{↵
x %= p;↵
if (x >= p)↵
x -= p;↵
if (x < 0)↵
x += p;↵
}↵
const int na = a.size(), nb = b.size();↵
int n, invN = 1;↵
for (n = 1; n < na + nb - 1; n <<= 1)↵
invN = ((invN & 1) ? (invN + M) : invN) >> 1;↵
vector<int> x(n, 0), y(n, 0);↵
std::copy(a.begin(), a.end(), x.begin());↵
std::copy(b.begin(), b.end(), y.begin());↵
fft(x);↵
fft(y);↵
for (int i = 0; i < n; ++i)↵
x[i] = (((static_cast<long long>(x[i]) * y[i]) % M) * invN) % M;↵
std::reverse(x.begin() + 1, x.end());↵
fft(x);↵
x.resize(na + nb - 1);↵
return x;↵
}↵
};↵
↵
Fft<998244353, 21, 31 * 31 * 31 * 31> muls;↵
↵
template<int32_t MOD>↵
struct ModInt↵
{↵
int32_t value;↵
↵
ModInt() : value(0)↵
{↵
}↵
↵
ModInt(long long v) : value(v % MOD)↵
{↵
if (value < 0)↵
value += MOD;↵
}↵
↵
ModInt(int32_t v): value(v % MOD)↵
{↵
if (value < 0)↵
value += MOD;↵
}↵
↵
ModInt operator+=(ModInt m)↵
{↵
value += m.value;↵
if (value >= MOD)↵
value -= MOD;↵
return value;↵
}↵
↵
ModInt operator-=(ModInt m)↵
{↵
value -= m.value;↵
if (value < 0)↵
value += MOD;↵
return value;↵
}↵
↵
ModInt operator*=(ModInt m)↵
{↵
value = (value * 1LL * m.value) % MOD;↵
return value;↵
}↵
↵
ModInt power(long long exp) const↵
{↵
if (exp == 0)↵
return 1;↵
ModInt res = (exp & 1 ? value : 1);↵
ModInt half = power(exp >> 1);↵
return res * half * half;↵
}↵
↵
ModInt operator/=(ModInt m) { return *this *= m.power(MOD - 2); }↵
↵
friend std::istream &operator>>(std::istream &is, ModInt &m)↵
{↵
is >> m.value;↵
return is;↵
}↵
↵
friend std::ostream &operator<<(std::ostream &os, const ModInt &m)↵
{↵
os << m.value;↵
return os;↵
}↵
↵
explicit operator int32_t() const { return value; }↵
↵
explicit operator long long() const { return value; }↵
↵
static int32_t mod() { return MOD; }↵
};↵
↵
template<int32_t MOD>↵
ModInt<MOD> operator+(ModInt<MOD> a, ModInt<MOD> b) { return a += b; }↵
↵
template<int32_t MOD, typename L>↵
ModInt<MOD> operator+(L a, ModInt<MOD> b) { return ModInt<MOD>(a) += b; }↵
↵
template<int32_t MOD, typename R>↵
ModInt<MOD> operator+(ModInt<MOD> a, R b) { return a += b; }↵
↵
template<int32_t MOD>↵
ModInt<MOD> operator-(ModInt<MOD> a, ModInt<MOD> b) { return a -= b; }↵
↵
template<int32_t MOD, typename L>↵
ModInt<MOD> operator-(L a, ModInt<MOD> b) { return ModInt<MOD>(a) -= b; }↵
↵
template<int32_t MOD, typename R>↵
ModInt<MOD> operator-(ModInt<MOD> a, R b) { return a -= b; }↵
↵
template<int32_t MOD>↵
ModInt<MOD> operator*(ModInt<MOD> a, ModInt<MOD> b) { return a *= b; }↵
↵
template<int32_t MOD, typename L>↵
ModInt<MOD> operator*(L a, ModInt<MOD> b) { return ModInt<MOD>(a) *= b; }↵
↵
template<int32_t MOD, typename R>↵
ModInt<MOD> operator*(ModInt<MOD> a, R b) { return a *= b; }↵
↵
template<int32_t MOD>↵
ModInt<MOD> operator/(ModInt<MOD> a, ModInt<MOD> b) { return a /= b; }↵
↵
template<int32_t MOD, typename L>↵
ModInt<MOD> operator/(L a, ModInt<MOD> b) { return ModInt<MOD>(a) /= b; }↵
↵
template<int32_t MOD, typename R>↵
ModInt<MOD> operator/(ModInt<MOD> a, R b) { return a /= b; }↵
↵
template<int32_t MOD>↵
bool operator==(ModInt<MOD> a, ModInt<MOD> b) { return a.value == b.value; }↵
↵
template<int32_t MOD, typename L>↵
bool operator==(L a, ModInt<MOD> b) { return a == b.value; }↵
↵
template<int32_t MOD, typename R>↵
bool operator==(ModInt<MOD> a, R b) { return a.value == b; }↵
↵
template<int32_t MOD>↵
bool operator!=(ModInt<MOD> a, ModInt<MOD> b) { return a.value != b.value; }↵
↵
template<int32_t MOD, typename L>↵
bool operator!=(L a, ModInt<MOD> b) { return a != b.value; }↵
↵
template<int32_t MOD, typename R>↵
bool operator!=(ModInt<MOD> a, R b) { return a.value != b; }↵
↵
using mint = ModInt<md>;↵
mint inv(mint x) { return 1 / x; }↵
↵
__int128 gcd(__int128 a, __int128 b, __int128 &x, __int128 &y)↵
{↵
if (b == 0)↵
{↵
x = 1;↵
y = 0;↵
return a;↵
}↵
__int128 d = gcd(b, a % b, y, x);↵
y -= a / b * x;↵
return d;↵
}↵
↵
__int128 inv(__int128 r, __int128 m)↵
{↵
__int128 x, y;↵
gcd(r, m, x, y);↵
return (x + m) % m;↵
}↵
↵
__int128 crt(__int128 r, __int128 n, __int128 c, __int128 m) { return r + ((c - r) % m + m) * inv(n, m) % m * n; }↵
const int m2 = 167772161, m3 = 469762049;↵
Fft<m2, 21, 147771621> muls2;↵
Fft<m3, 21, 297449090> muls3;↵
↵
vector<mint> operator*(vector<mint> a, vector<mint> b)↵
{ ///modulo-dependent convolution↵
if (a.empty() || b.empty())↵
return {};↵
if (md == 998244353)↵
{↵
vector<int> a1(a.size());↵
for (int i = 0; i < a.size(); ++i)↵
a1[i] = a[i].value;↵
vector<int> b1(b.size());↵
for (int i = 0; i < b.size(); ++i)↵
b1[i] = b[i].value;↵
vector<int> c1 = muls.convolution(a1, b1);↵
vector<mint> c;↵
for (int x: c1)↵
c.app(x);↵
return c;↵
}↵
else↵
{↵
vector<int> a1(a.size());↵
for (int i = 0; i < a.size(); ++i)↵
a1[i] = a[i].value;↵
vector<int> b1(b.size());↵
for (int i = 0; i < b.size(); ++i)↵
b1[i] = b[i].value;↵
vector<int> c1 = muls.convolution(a1, b1);↵
vector<int> c2 = muls2.convolution(a1, b1);↵
vector<int> c3 = muls3.convolution(a1, b1);↵
assert(c1.size()==c2.size() && c2.size()==c3.size());↵
vector<int> c4(c1.size());↵
for (int i = 0; i < c1.size(); ++i)↵
{↵
__int128 ost1 = c1[i];↵
__int128 m1 = 998244353;↵
__int128 ost2 = c2[i];↵
__int128 ost3 = c3[i];↵
__int128 ost = crt(crt(ost1, m1, ost2, m2), m1 * 1LL * m2, ost3, m3);↵
c4[i] = (ost % md);↵
}↵
vector<mint> c;↵
for (int x: c4)↵
c.app(x);↵
return c;↵
}↵
}↵
↵
vector<vector<mint> > gaussbasis(vector<vector<mint> > A) ///returns basis of Av=0↵
{↵
int n = A.size();↵
int m = A[0].size();↵
int bi = 0;↵
for (int i = 0; i < n; ++i)↵
{↵
if (bi == m)↵
break;↵
for (int j = i; j < n; ++j)↵
{↵
if (A[j][bi] != 0)↵
{↵
if (j != i) { swap(A[i], A[j]); }↵
break;↵
}↵
}↵
if (A[i][bi] != 0)↵
{↵
mint o = 1 / A[i][bi];↵
for (int j = i + 1; j < n; ++j)↵
{↵
mint we = (A[j][bi] * o);↵
for (int k = bi; k < m; ++k)↵
{↵
A[j][k] -= we * A[i][k];↵
}↵
}↵
}↵
else↵
{↵
++bi;↵
--i;↵
continue;↵
}↵
}↵
vector<int> indices(m);↵
iota(all(indices), 0);↵
for (int i = n - 1; i >= 0; --i)↵
{↵
int bi = 0;↵
while (bi < m && A[i][bi] == 0) { ++bi; }↵
if (bi < m)↵
{↵
indices.erase(find(all(indices), bi));↵
}↵
}↵
vector<vector<mint> > v(indices.size(), vector<mint>(m, 0));↵
for (int i = 0; i < indices.size(); ++i)↵
{↵
v[i][indices[i]] = 1;↵
}↵
for (int i = n - 1; i >= 0; --i)↵
{↵
int bi = 0;↵
while (bi < m && A[i][bi] == 0) { ++bi; }↵
if (bi == m)↵
continue;↵
for (int k = 0; k < indices.size(); ++k)↵
{↵
mint cur = 0;↵
for (int j = bi + 1; j < m; ++j)↵
{↵
cur -= A[i][j] * v[k][j];↵
}↵
v[k][bi] = cur / A[i][bi];↵
}↵
}↵
return v;↵
}↵
↵
optional<vector<mint> > gauss(vector<vector<mint> > A, vector<mint> b) ///returns v such that Av=b↵
{↵
int n = A.size();↵
assert(b.size()==n);↵
int m = A[0].size();↵
int bi = 0;↵
for (int i = 0; i < n; ++i)↵
{↵
if (bi == m)↵
break;↵
for (int j = i; j < n; ++j)↵
{↵
if (A[j][bi] != 0)↵
{↵
if (j != i)↵
{↵
swap(A[i], A[j]);↵
swap(b[i], b[j]);↵
}↵
break;↵
}↵
}↵
if (A[i][bi] != 0)↵
{↵
mint o = inv(A[i][bi]);↵
for (int j = i + 1; j < n; ++j)↵
{↵
mint we = (A[j][bi] * o);↵
b[j] -= we * b[i];↵
for (int k = bi; k < m; ++k)↵
{↵
A[j][k] -= we * A[i][k];↵
}↵
}↵
}↵
else↵
{↵
++bi;↵
--i;↵
continue;↵
}↵
}↵
vector<mint> v(m);↵
for (int i = n - 1; i >= 0; --i)↵
{↵
int bi = 0;↵
while (bi < m && A[i][bi] == 0) { ++bi; }↵
if (bi == m)↵
{↵
if (b[i] != 0) { return nullopt; }↵
else { continue; }↵
} {↵
mint cur = b[i];↵
for (int j = bi + 1; j < m; ++j)↵
{↵
cur -= A[i][j] * v[j];↵
}↵
v[bi] = cur * inv(A[i][bi]);↵
}↵
}↵
return v;↵
}↵
↵
optional<vector<vector<mint> > > findPrecursion(vector<mint> a)↵
{ ///finds P-recursion of a given sequence A by gauss↵
for (int snd = 0; snd <= 20; ++snd)↵
{↵
for (int n = 1; n <= snd - 1; ++n)↵
{↵
vector<vector<mint> > A;↵
int d = snd - n;↵
int eq = ((int) (a.size())) - (n - 1);↵
if (eq < n * d) { continue; }↵
for (int i = n - 1; i < a.size(); ++i)↵
{↵
vector<mint> u;↵
for (int j = 0; j < n; ++j)↵
{↵
mint de = 1;↵
for (int k = 0; k < d; ++k)↵
{↵
u.app(a[i - j] * de);↵
de *= i;↵
}↵
}↵
A.app(u);↵
}↵
vector<vector<mint> > zx = gaussbasis(A);↵
if (zx.empty())↵
continue;↵
//debug(n, d);↵
vector<mint> ans = zx[0];↵
vector<vector<mint> > res;↵
for (int j = 0; j < n; ++j)↵
{↵
res.app({});↵
for (int k = 0; k < d; ++k)↵
{↵
res[j].app(ans[j * d + k]);↵
}↵
}↵
return res;↵
}↵
}↵
return nullopt;↵
}↵
↵
optional<vector<mint> > evaluatePrecursion(vector<mint> a, vector<vector<mint> > rec, int sz)↵
{ ///a(0),...,a(a.size()-1) -> (by P-recursion rec) a(0),...,a(sz-1)↵
int n = rec.size();↵
int d = rec[0].size();↵
int given = a.size();↵
if (given >= sz)↵
{↵
a.resize(sz);↵
return a;↵
}↵
if (a.size() < n) { return nullopt; }↵
vector<mint> tore;↵
for (int i = given; i < sz; ++i)↵
{↵
mint de = 1;↵
mint s = 0;↵
for (int k = 0; k < d; ++k)↵
{↵
s += de * rec[0][k];↵
de *= i;↵
}↵
if (s == 0) { return nullopt; }↵
tore.app(s);↵
}↵
vector<mint> pref(tore.size() + 1);↵
pref[0] = 1;↵
for (int i = 0; i < tore.size(); ++i) { pref[i + 1] = pref[i] * tore[i]; }↵
mint pro = pref[tore.size()];↵
mint invpro = 1 / pro;↵
mint cur = invpro;↵
vector<mint> invtore(tore.size());↵
for (int i = tore.size() - 1; i >= 0; --i)↵
{↵
invtore[i] = cur * pref[i];↵
cur *= tore[i];↵
}↵
for (int i = given; i < sz; ++i)↵
{↵
mint chi = 0;↵
for (int j = 1; j < n; ++j)↵
{↵
mint de = 1;↵
for (int k = 0; k < d; ++k)↵
{↵
chi += a[i - j] * de * rec[j][k];↵
de *= i;↵
}↵
}↵
a.app(((mint) (0)) - chi * invtore[i - given]);↵
}↵
return a;↵
}↵
↵
mint value(vector<mint> a, mint x)↵
{ ///A(x)↵
mint de = 1;↵
mint ans = 0;↵
for (int i = 0; i < a.size(); ++i)↵
{↵
ans += a[i] * de;↵
de *= x;↵
}↵
return ans;↵
}↵
↵
vector<mint> shiftofsamplingpoints(vector<mint> a)↵
{ ///P(0),...,P(t) we want to compute P(0),...,P(4t+1)↵
int t = a.size() - 1;↵
vector<mint> fact(4 * t + 2);↵
fact[0] = 1;↵
for (int i = 1; i < 4 * t + 2; ++i)↵
fact[i] = fact[i - 1] * i;↵
vector<mint> invf(4 * t + 2);↵
invf[4 * t + 1] = 1 / fact[4 * t + 1];↵
for (int i = 4 * t; i >= 0; --i) { invf[i] = (invf[i + 1] * (i + 1)); }↵
assert(invf[0]==1);↵
vector<mint> invm(4 * t + 2, 0);↵
for (int i = 1; i < 4 * t + 2; ++i) { invm[i] = fact[i - 1] * invf[i]; }↵
vector<mint> values(t + 1, 0);↵
for (int k = 0; k <= t; ++k)↵
{↵
mint o = 1;↵
if ((t - k) % 2 == 1) { o = (((mint) (0)) - 1); }↵
values[k] = (a[k] * invf[k] * invf[t - k] * o);↵
}↵
vector<mint> h = invm * values;↵
vector<mint> res;↵
for (int i = 0; i <= t; ++i)↵
{↵
res.app(a[i]);↵
}↵
for (int x = t + 1; x <= 4 * t + 1; ++x)↵
{↵
mint ans = fact[x];↵
ans *= invf[x - t - 1];↵
ans *= h[x];↵
res.app(ans);↵
}↵
return res;↵
}↵
↵
optional<mint> evaluatePrecursionfast(vector<mint> a, vector<vector<mint> > rec, int id)↵
{ ///a(0),...,a(a.size()-1) -> (by P-recursion rec) a(id), O(sqrt(id)*log(id))↵
if (id < a.size())↵
return a[id];↵
int n = rec.size();↵
int d = 1;↵
for (auto &v: rec) { d = max(d, ((int) (v.size() - 1))); }↵
if (a.size() < n - 1) { return nullopt; }↵
if (n == 1) { return 0; }↵
int l = n - 1;↵
int shift = 0;↵
while (a.size() > l)↵
{↵
a.erase(a.begin());↵
++shift;↵
--id;↵
}↵
int u = 0;↵
while ((1LL << u) * (1LL << u) <= id) { ++u; }↵
int B = (1 << u);↵
vector<mint> S;↵
vector<vector<vector<mint> > > A(l, vector<vector<mint> >(l));↵
int sz = d;↵
for (int k = 0; k <= d; ++k)↵
{↵
S.app(value(rec[0], k + l + shift));↵
}↵
for (int i = 0; i < l - 1; ++i)↵
{↵
for (int j = 0; j < l; ++j)↵
{↵
if (j == i + 1)↵
{↵
for (int k = 0; k <= d; ++k)↵
{↵
A[i][j].app(value(rec[0], k + l + shift));↵
}↵
}↵
else↵
{↵
for (int k = 0; k <= d; ++k)↵
{↵
A[i][j].app(0);↵
}↵
}↵
}↵
}↵
for (int j = 0; j < l; ++j)↵
{↵
for (int k = 0; k <= d; ++k)↵
{↵
A[l - 1][j].app(((mint) (0)) - value(rec[l - j], k + l + shift));↵
}↵
}↵
for (int s = 0; s < u; ++s)↵
{↵
S = shiftofsamplingpoints(S);↵
assert(S.size()==4*sz+2);↵
for (int i = 0; i < l; ++i)↵
{↵
for (int j = 0; j < l; ++j)↵
{↵
A[i][j] = shiftofsamplingpoints(A[i][j]);↵
assert(A[i][j].size()==4*sz+2);↵
}↵
}↵
vector<vector<vector<mint> > > newA(l, vector<vector<mint> >(l, vector<mint>(2 * sz + 1, 0)));↵
for (int k = 0; k <= 2 * sz; ++k)↵
{↵
for (int ii = 0; ii < l; ++ii)↵
{↵
for (int jj = 0; jj < l; ++jj)↵
{↵
for (int kk = 0; kk < l; ++kk)↵
{↵
newA[ii][kk][k] += A[ii][jj][2 * k + 1] * A[jj][kk][2 * k];↵
}↵
}↵
}↵
}↵
vector<mint> newS(2 * sz + 1, 0);↵
for (int k = 0; k <= 2 * sz; ++k) { newS[k] = S[2 * k] * S[2 * k + 1]; }↵
sz *= 2;↵
S = newS;↵
A = newA;↵
}↵
int k = (id - l) / B;↵
assert(k>=0); ///id>=l+1 here↵
mint pro = 1;↵
vector<mint> v;↵
for (int i = 0; i < l; ++i)↵
v.app(a[i]);↵
for (int i = 0; i < k; ++i)↵
{↵
vector<mint> newv(l, 0);↵
vector<vector<mint> > M(l, vector<mint>(l, 0));↵
for (int ii = 0; ii < l; ++ii)↵
{↵
for (int jj = 0; jj < l; ++jj)↵
{↵
assert(i<A[ii][jj].size());↵
M[ii][jj] = A[ii][jj][i];↵
}↵
}↵
pro *= S[i];↵
for (int ii = 0; ii < l; ++ii)↵
{↵
for (int jj = 0; jj < l; ++jj)↵
{↵
newv[ii] += M[ii][jj] * v[jj];↵
}↵
}↵
v = newv;↵
}↵
int cur = k * B + l - 1;↵
assert(cur<id);↵
while (cur < id)↵
{↵
mint newval = 0;↵
for (int i = 0; i < l; ++i)↵
{↵
newval -= v[i] * value(rec[l - i], cur + 1 + shift);↵
}↵
mint de = value(rec[0], cur + 1 + shift);↵
pro *= de;↵
for (mint &x: v) { x *= de; }↵
v.erase(v.begin());↵
v.app(newval);↵
++cur;↵
}↵
if (pro == 0)↵
return nullopt;↵
return v.back() / pro;↵
}↵
↵
optional<vector<mint> > sequenceextender(vector<mint> a, int sz)↵
{ ///finds P-recursion, and if was found, calculates a(0),...,a(sz-1)↵
if (a.size() >= sz)↵
{↵
a.resize(sz);↵
return a;↵
}↵
auto uu = findPrecursion(a);↵
if (!uu)↵
return nullopt;↵
auto rec = (*uu);↵
auto ans = evaluatePrecursion(a, rec, sz);↵
if (!ans)↵
return nullopt;↵
return (*ans);↵
}↵
↵
optional<mint> fastgetvaluebyid(vector<mint> a, int id)↵
{ ///finds P-recursion, and if was found, calculates a(id) in O(sqrt(id)*log(id))↵
if (a.size() > id) { return a[id]; }↵
auto uu = findPrecursion(a);↵
if (!uu)↵
return nullopt;↵
auto rec = (*uu);↵
auto ans = evaluatePrecursionfast(a, rec, id);↵
if (!ans)↵
return nullopt;↵
return (*ans);↵
}↵
↵
optional<mint> optimalgetvaluebyid(vector<mint> a, int id)↵
{ ///finds P-recursion, and if was found, calculates a(id) by choosing optimal of O(id) method and O(sqrt(id)*log(id)) method↵
if (a.size() > id) { return a[id]; }↵
auto uu = findPrecursion(a);↵
if (!uu)↵
return nullopt;↵
auto rec = (*uu);↵
int n = rec.size();↵
int d = rec[0].size();↵
double C = 1;↵
if (md != 998244353)↵
C = 3;↵
double val1 = sqrt(id) * log(id) * C * n * n * d + sqrt(id) * n * n * n * d;↵
double val2 = n * 1.0 * d * 1.0 * id;↵
//debug(val1, val2);↵
if (val1 < val2)↵
{↵
//debug("fastgetvalue");↵
auto ans = evaluatePrecursionfast(a, rec, id);↵
if (!ans)↵
return nullopt;↵
return (*ans);↵
}↵
else↵
{↵
//debug("sequenceextender");↵
auto ans = evaluatePrecursion(a, rec, id + 1);↵
if (!ans)↵
return nullopt;↵
return (*ans)[id];↵
}↵
}↵
↵
vector<vector<mint> > transpose(vector<vector<mint> > a)↵
{ ///transposes the table A↵
if (a.empty())↵
return a;↵
int n = a.size();↵
int m = a[0].size();↵
vector<vector<mint> > b(m, vector<mint>(n, 0));↵
for (int i = 0; i < n; ++i)↵
for (int j = 0; j < m; ++j)↵
b[j][i] = a[i][j];↵
return b;↵
}↵
↵
/// If the size of a is big TL (too many Gauss), if the size of a is small WA (not enough for finding the P-recursion), you should keep balance↵
optional<vector<mint> > extendtableslow(vector<vector<mint> > a, vector<pair<int, int> > que)↵
{ /// extends table A, finding A(que[i].first,que[i].second), in a O(sum(que[i]))↵
if (que.empty())↵
{↵
vector<mint> ret = {};↵
return ret;↵
}↵
int ma = 0;↵
for (auto [i,j]: que) { ma = max(ma, j); }↵
int n = a.size();↵
int m = a[0].size();↵
vector<vector<mint> > ex(n);↵
for (int i = 0; i < n; ++i)↵
{↵
auto h = sequenceextender(a[i], ma + 1);↵
if (!h) { return nullopt; }↵
ex[i] = (*h);↵
}↵
vector<mint> res;↵
for (auto [i,j]: que)↵
{↵
vector<mint> e;↵
for (int k = 0; k < n; ++k)↵
{↵
e.app(ex[k][j]);↵
}↵
auto h = sequenceextender(e, i + 1);↵
if (!h) { return nullopt; }↵
res.app((*h)[i]);↵
}↵
return res;↵
}↵
↵
optional<vector<mint> > extendtablefast(vector<vector<mint> > a, vector<pair<int, int> > que)↵
{ /// extends table A, finding A(que[i].first,que[i].second)↵
if (que.empty())↵
{↵
vector<mint> ret = {};↵
return ret;↵
}↵
int n = a.size();↵
int m = a[0].size();↵
vector<vector<mint> > ex(n);↵
vector<vector<mint> > rec[n];↵
for (int k = 0; k < n; ++k)↵
{↵
auto uu = findPrecursion(a[k]);↵
if (!uu) { return nullopt; }↵
rec[k] = (*uu);↵
}↵
vector<mint> res;↵
for (auto [i,j]: que)↵
{↵
vector<mint> e;↵
for (int k = 0; k < n; ++k)↵
{↵
auto uu = evaluatePrecursionfast(a[k], rec[k], j);↵
if (!uu) { return nullopt; }↵
e.app(*uu);↵
}↵
auto h = optimalgetvaluebyid(e, i);↵
if (!h) { return nullopt; }↵
res.app(*h);↵
}↵
return res;↵
}↵
↵
optional<vector<mint> > extendtable(vector<vector<mint> > a, vector<pair<int, int> > que)↵
{ /// extends table A, finding A(que[i].first,que[i].second)↵
double C = 1;↵
if (md != 998244353) { C = 3; }↵
double op1 = 0;↵
double op2 = 0;↵
for (auto [i,j]: que)↵
{↵
op1 += C * 1.0 * ((int) (a.size())) * 1.0 * sqrt(i + 1) * 1.0 * log(i + 2) * 5 * 5 * 5;↵
op1 += C * sqrt(j + 1) * log(j + 2) * 5 * 5 * 5;↵
}↵
for (auto [i,j]: que)↵
{↵
op2 += C * 1.0 * ((int) (a.size())) * 1.0 * (i + 1) * 1.0 * 5 * 5;↵
op2 += C * 1.0 * (j + 1) * 1.0 * 5 * 5;↵
}↵
//debug(op1, op2);↵
if (op1 < op2) { return extendtablefast(a, que); }↵
else { return extendtableslow(a, que); }↵
}↵
↵
optional<vector<mint> > getcolumnoftable(vector<vector<mint> > a, int col, int size)↵
{ /// get A(0,col),A(1,col),...,A(size-1,col)↵
int n = a.size();↵
int m = a[0].size();↵
vector<mint> e;↵
for (int i = 0; i < n; ++i)↵
{↵
auto h = sequenceextender(a[i], col + 1);↵
if (!h) { return nullopt; }↵
e.app((*h)[col]);↵
}↵
auto h = sequenceextender(e, size);↵
if (!h)↵
return nullopt;↵
return *h;↵
}↵
↵
optional<vector<mint> > getrowoftable(vector<vector<mint> > a, int row, int size)↵
{ /// get A(row,0),A(row,1),...,A(row,size-1)↵
return getcolumnoftable(transpose(a), row, size);↵
}↵
↵
void test()↵
{↵
vector<mint> v1 = {1, 1, 3, 7, 19, 51, 141, 393, 1107, 3139}; ///(1+x+x^2)^n [x^n]↵
auto uu1 = sequenceextender(v1, 30);↵
if (uu1)↵
{↵
auto u1 = (*uu1);↵
debugv(u1);↵
}↵
vector<mint> v2 = {1, 30, 465, 4930, 40020, 264306, 1474795, 7133130, 30462615, 116470380}; ///(1+x+x^2)^30 [x^n]↵
auto uu2 = sequenceextender(v2, 70);↵
if (uu2)↵
{↵
auto u2 = (*uu2);↵
debugv(u2);↵
}↵
vector<mint> v3 = {0, 567646151, 513265721, 604121291, 715018514, 398975714, 610803800, 499563577, 491416403,↵
913506524↵
}; ///s(30,n), not D-finite↵
auto uu3 = sequenceextender(v3, 35);↵
if (uu3)↵
{↵
auto u3 = (*uu3);↵
debugv(u3);↵
}↵
vector<mint> v4 = {0, 1, 2, 9, 44, 265, 1854, 14833, 133496, 1334961, 14684570};↵
///n!*x*e^(-x) [x^n] (number of permutations of size n without stable points)↵
auto uu4 = sequenceextender(v4, 20);↵
if (uu4)↵
{↵
auto u4 = (*uu4);↵
debugv(u4);↵
}↵
auto uu5 = fastgetvaluebyid({1, 1, 2, 6, 24, 120, 720}, 11); ///factorials ,by id↵
if (uu5)↵
{↵
auto u5 = (*uu5);↵
debug(u5);↵
}↵
auto uu6 = fastgetvaluebyid(v1, 29); ///(1+x+x^2)^n [x^n], by id↵
if (uu6)↵
{↵
auto u6 = (*uu6);↵
debug(u6);↵
}↵
auto uu7 = fastgetvaluebyid(v1, 500000000); ///(1+x+x^2)^n [x^n], by id↵
if (uu7)↵
{↵
auto u7 = (*uu7);↵
debug(u7);↵
}↵
auto uu8 = optimalgetvaluebyid(v1, 500000000); ///(1+x+x^2)^n [x^n], by id↵
if (uu8)↵
{↵
auto u8 = (*uu8);↵
debug(u8);↵
}↵
auto uu9 = optimalgetvaluebyid({1, 1, 2, 6, 24, 120, 720}, 11); ///factorials ,by id↵
if (uu9)↵
{↵
auto u9 = (*uu9);↵
debug(u9);↵
}↵
auto uu10 = optimalgetvaluebyid({1, 1, 2, 6, 24, 120, 720}, 998244352); ///factorials ,by id↵
if (uu10)↵
{↵
auto u10 = (*uu10);↵
debug(u10);↵
}↵
vector<vector<mint> > table1 =↵
{{1, 2, 3, 4, 5}, {2, 3, 4, 5, 6}, {3, 4, 5, 6, 7}, {4, 5, 6, 7, 8}, {5, 6, 7, 8, 9}}; ///f(i,j)=i+j+1↵
vector<pair<int, int> > que1 = {{1, 1}, {0, 0}, {5, 7}, {11342, 1333}};↵
auto uu11 = extendtablefast(table1, que1);↵
auto uu12 = extendtableslow(table1, que1);↵
auto uu13 = extendtable(table1, que1);↵
debug((bool) (uu11));↵
debug((bool) (uu12));↵
debug((bool) (uu13));↵
if (uu11 && uu12 && uu13)↵
{↵
auto u11 = (*uu11);↵
auto u12 = (*uu12);↵
auto u13 = (*uu13);↵
debugv(u11);↵
debugv(u12);↵
debugv(u13);↵
}↵
auto uu14 = getcolumnoftable(table1, 10, 100);↵
if (uu14)↵
{↵
auto u14 = (*uu14);↵
debugv(u14);↵
}↵
auto uu15 = getrowoftable(table1, 50, 70);↵
if (uu15)↵
{↵
auto u15 = (*uu15);↵
debugv(u15);↵
}↵
///exp((x+y)/((1-x)(1-y)), I removed it from examples↵
}↵
↵
int32_t main()↵
{↵
ios_base::sync_with_stdio(false);↵
cin.tie(0);↵
test();↵
return 0;↵
}↵
```↵
</spoiler>↵
↵
Code works for all prime moduli. Explanation:↵
↵
* `findPrecursion` finds P-recursion ↵
↵
* `evaluatePrecursion` evaluates P-recursion from the initial values if possible ↵
↵
* `evaluatePrecursionfast` evaluates the element of a 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` and `getrowoftable` get column or row of an extended table ↵
↵
↵
A "small" problem↵
-----------------↵
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)+\ldots+a_{n-l}P_{l}(n)=0$.↵
↵
You know $a(0),\ldots,a(l-1)$ and you want to extend this sequence to $a(0),\ldots,a(N)$. But what if $P_{0}(k)=0$ for some $l \leq k \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 functions when we can use it, it can't be zero. ↵
↵
Example where it matters↵
------------------------↵
↵
Let's try to compute $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. ↵
↵
Examples of the problems↵
------------------------↵
↵
**1)** [Atcoder Regular 202 problem D](https://atcoder.jp/contests/arc202/tasks/arc202_d)↵
↵
<spoiler>↵
Here if you are solving it in $O((H+w+T)^{1.5})$ using the reflection principle you have a subproblem to calculate $(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,\ldots,k+C$, and then extend it to the whole sequence in linear time.↵
</spoiler>↵
↵
**2)** We can also solve [problem: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: [submission:330486092]).↵
↵
I think there are some other examples of problems that can be solved by this method. ↵
↵
Conclusion (Warning)↵
====================↵
↵
When I got it at first, it seemed to me that "every" function in combinatorics is D-finite, but it is actually wrong (see before). I think that one should use this algorithm only if one has the proof that some function is D-finite, otherwise it can be a waste of time. So I don't recommend using this algorithm as magic.↵
↵
↵
↵
↵
Thanks to [user:orz,2025-07-25] for editing this blog!
----------↵
↵
7 months ago [user:Elegia,2025-07-25] published a blog about D-finite functions ([part1](https://codeforces.me/blog/entry/138078), [part2](https://codeforces.me/blog/entry/138113)) in which he showed some examples of problems that can be solved by D-finite functions.↵
↵
Here I will explain a fast evaluation algorithm of the coefficients of a D-finite function.↵
↵
Let $K$ be some field of characteristic $0$. Let $K[ [x_{1},\ldots,x_{d}] ]$ be the ring of formal power series on it. Let $K[x_{1},\ldots,x_{d}]$ be the polynomials of it. Let $K(x_{1},\ldots,x_{d})=\mathop{\mathrm{Frac}}(K[x_{1},\ldots,x_{d}])$ be the field of rational functions over it. Let $ [ x_{1}^{i_{1}}\ldots x_{d}^{i_{d}} ] f$ be the coefficient in $f$ of the monomial $x_{1}^{i_{1}}\ldots x_{d}^{i_{d}}$.↵
↵
Definitions↵
----------↵
↵
A function $f \in K[ [x_{1},\ldots,x_{d}] ]$ is called _D-finite_ if for each $i$ there exists a positive integer $n_{i}$ and polynomials (not all zeroes) $P_{i,0}(x_{1}, \ldots,x_{d}),\ldots,P_{i,n_{i}}(x_{1},\ldots,x_{d})$ such that $P_{i,0}f+P_{i,1}\frac{df}{dx_{i}}+\ldots+P_{i,n_{i}}\frac{d^{n_{i}}(f)}{(dx_{i})^{n_{i}}}=0$.↵
↵
An infinite table of coefficients $a(i_{1},\ldots,i_{d})$ is called _P-recursive_ if there exists $k \in \mathbb{N}$ such that for each $j \in \\{ 1,\ldots,d \\} $ there exists a mapping $p: \\{ 0,\ldots,k \\} ^{d} \rightarrow K[x]$ such that at least one $p_{v} \neq 0$ and for all $i_{1},\ldots,i_{d} \geq k$ $\sum\limits_{v \in \\{0,\ldots,k\\}^{d}}^{} p_{v}(i_{j}) a(i_{1}-v_{1},\ldots,i_{d}-v_{d}) = 0$. ↵
↵
In particular, for $d=1$ it means that there exist such polynomials (not all zeroes) $P_{i}$: $ \forall n \geq k : \sum\limits_{0 \leq i \leq k}^{} P_{i}(n)a(n-i)=0$.↵
↵
Properties:↵
-----------↵
↵
**1** The definition is equivalent to the "linear space of partial derivatives of $f$ over $K(x_{1},\ldots,x_{d})$ has finite dimension". ↵
↵
In particular, every rational function is obviously D-finite↵
↵
**2** They form a linear algebra over $K(x_{1},\ldots,x_{d})$ (we can multiply them by a scalar, add them and multiply them, and the result will be D-finite)↵
↵
**3** Every algebraic over $K(x_{1},\ldots,x_{d})$ function is D-finite ↵
↵
**4** If $u_{1},\ldots,u_{d} \in K[ [u_{1},\ldots,u_{d}] ]$ are algebraic, and $f$ is D-finite and composition $f(u_{1},\ldots,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},\ldots,i_{d}}^{} x_{1}^{i_{12}}x_{3}^{i_{3}}\ldots x_{d}^{i_{d}} ([x_{1}^{i_{12}}x_{2}^{i_{12}}x_{3}^{i_{3}}\ldots 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 the diagonal) ↵
↵
**6(Lipshitz, 1989)** The function $f=\sum\limits a(i_{1},\ldots,i_{d}) x_{1}^{i_{1}}\ldots x_{d}^{i_{d}}$ is D-finite if and only if its coefficients are P-recursive.↵
↵
**7(Lipshitz, 1988)** If $f$ is D-finite, then for each $i_{1}$: $\sum\limits_{i_{2},i_{3},\ldots,i_{d}}^{} x_{2}^{i_{2}}x_{3}^{i_{3}}\ldots x_{d}^{i_{d}} ([x_{1}^{i_{1}}x_{2}^{i_{2}}x_{3}^{i_{3}}\ldots x_{d}^{i_{d}}] f)$ is D-finite. Moreover, the number of polynomials and their degrees in the D-finiteness relations are uniformly bounded. ↵
↵
Note: D-finite functions are not a field ↵
↵
Let's prove **7**↵
↵
<spoiler summary="Proof">↵
↵
Let's consider $g(x_{1},\ldots,x_{d})=x_{1}^{-i_{1}-1}f(x_{1},\ldots,x_{d})$. Then $f(x_{1},\ldots,x_{d})=x^{i_{1}+1}g(x_{1},\ldots,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},\ldots,x_{d})(d/dx_{1})^{k} g + P_{k-1}(x_{1},\ldots,x_{d})(d/dx_{1})^{k-1} g + \ldots +P_{0}(x_{1},\ldots,x_{d})g = 0$ and $Q_{l}(x_{1},\ldots,x_{d})(d/dx_{2})^{l} g + Q_{l-1}(x_{1},\ldots,x_{d})(d/dx_{2})^{l-1} g + \ldots +Q_{0}(x_{1},\ldots,x_{d})g = 0$. Let $L(x_{1},\ldots,x_{d})=P_{k}Q_{l}$. ↵
↵
Let $d=\max(\deg(P_{0},\ldots,P_{k}),Q_{0},\ldots,Q_{l})$. ↵
↵
If we write $(d/dx_{1})^{u}(d/dx_{2}^{v})(P_{k}(x_{1},\ldots,x_{d})(d/dx_{1})^{k} g + P_{k-1}(x_{1},\ldots,x_{d})(d/dx_{1})^{k-1} g + \ldots +P_{0}(x_{1},\ldots,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 + \ldots \text{(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'<s+t} p_{s',t'}(x_{1},\ldots,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}}\ldots x_{d}^{\alpha_{d}}(d/dx_{1})^{m_{1}}(d/dx_{2})^{m_{2}} g$, where $N$ is fixed, $\alpha_{2}+\alpha_{3}+\ldots+\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}}\ldots x_{d}^{\alpha_{d}}\sum\limits_{s,t}^{} P_{s,t}(x_{1},\ldots,x_{d})((d/dx_{i})^{s}(d/dx_{2})^{t}) g$ where $s<k,t<l$ and $P_{s,t}(x_{1},\ldots,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}}\ldots 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},\ldots,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}}\ldots 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}}\ldots 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}}\ldots x_{d}^{\alpha_{d}}(d/dx_{1})^{m_{1}}(d/dx_{2})^{m_{2}} g}=0$ (not all $c$ are zeroes), let $m$ be the minimal $m_{1}$ for which there is $c_{\ldots} \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}}\ldots 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 is 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}}\ldots 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},\ldots,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).↵
↵
</spoiler>↵
↵
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. ↵
↵
What functions are D-finite?↵
----------------------------↵
↵
Can you guess which of these functions are D-finite (that is, the sequence/table of coefficients is P-recursive)? ↵
↵
(1) $a_{n} = \[ x^{n+1} \] (1+x+x^{2})^{n}$ ($f(x)=\sum\limits_{}^{} a_{n}x^{n}$)↵
↵
(2) $a_{n}=n!$ ($f(x)=\sum\limits_{}^{} n!x^{n}$)↵
↵
(3) $a_{n,k}=S(n,k)$ [Stirling numbers of the second kind](https://en.wikipedia.org/wiki/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)}$ ↵
↵
(It is not so hard to see, that the D-finiteness of the OGF of some sequence/table is equivalent to the D-finiteness of EGF of this sequence/table)↵
↵
(4) [Eulerian numbers](https://en.wikipedia.org/wiki/Eulerian_number) $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!}$. ↵
↵
<spoiler summary="Answer">↵
↵
$(1),(2),(5)$ are D-finite, $(3),(4)$ I believe are not D-finite. ↵
↵
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, the reason it may not be D-finite is that we get the exponent of non-algebraic thing (and we can prove it is not D-finite). ↵
↵
For $(4)$ the problem, the reason it may not be D-finite is that when we take the last inverse of the whole thing, this thing is not algebraic (and we can prove it is not D-finite). ↵
↵
$(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)}$. ↵
↵
</spoiler>↵
↵
The D-finite functions with $d=1$ are just P-recursive sequences. And we can evaluate the 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))$. ↵
↵
<spoiler summary="Algorithm"> ↵
↵
The algorithm is pretty similar to the one for factorials. ↵
↵
Let $l$ be the length of the P-recursive relation ($P_{l}(n)a_{n-l}+\ldots+P_{0}(n)a_{n}=0$). And let $d=\max(\deg(P_{i}))$. We know $a_{0},\ldots,a_{l-1}$, we want to compute $a_{n}$. ($l$ and $d$ are assumed 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},\ldots,a_{kB+l-1}$, then we can easily compute $a_{n}$ in $O(\sqrt{n})$. ↵
↵
Suppose we have some $a_{x},\ldots,a_{x+l-1}$ and we want to compute $a_{x+\mathrm{diff}},\ldots,a_{x+\mathrm{diff}+l-1}$, then we will have↵
↵
↵
↵
$$\begin{pmatrix} a_{x+\mathrm{diff}}\\ a_{x+\mathrm{diff}+1}\\ . \\ . \\ . \\ a_{x+\mathrm{diff}+l-1} \end{pmatrix} = \frac{1}{S(x)}\begin{pmatrix} P_{00} (x) & P_{01}(x) \ldots & P_{0(l-1)}(x)\\ P_{10} (x) & P_{11}(x) \ldots & P_{1(l-1)}(x)\\ . \\ . \\ . \\ P_{(l-1)0} (x) & P_{(l-1)1}(x) \ldots & 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)\ldots P_{0}(x+\mathrm{diff}+l-1)$ and $\max(\deg(S),\deg(P_{i,j})) \leq \mathrm{diff} \cdot d$. ↵
↵
We will denote the polynomials for a given $\mathrm{diff}$ as $P_{\mathrm{diff},i,j}$ and $S_{\mathrm{diff}}$. ↵
↵
Let's store the polynomial of degree $d$ as an array of values in a points $0,K,2K,\ldots,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),\ldots,P_{2^s,i,j}(2^{2s} \cdot d)$. And same with $S$, that is $S_{2^s}(0),\ldots,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 are assuming $l$ and $d$ to be constants, then it is just $O(\sqrt{n}\log(n))$. ↵
↵
</spoiler>↵
↵
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 $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 **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),\ldots,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 say "Let's compute the coefficient of this P-recursive sequence", how do I know that the P-recursive relation holds? ↵
↵
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 first $k$ coefficients (for some $k$) 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 polynomials of $u$, and then we can try to evaluate them from some interpolation of a rational function, though I have not implemented it yet. ↵
↵
Example of this algorithm: ↵
↵
Assume 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}) \pmod{998\,244\,353}$. ↵
↵
This function $f$ is D-finite. ↵
↵
Let's compute its values for $0 \leq n <25$ and $0 \leq m < 25$. Then we can extend ↵
↵
$f(\cdot,0),f(\cdot,1),f(\cdot,2),\ldots,f(\cdot,24)$ to $f(10^{7},0),f(10^{7},1),\ldots,f(10^{7},24)$ by P-recursion. Then we can extend $f(10^{7},\cdot)$ from the 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. ↵
↵
Note: for $f(10^{8},10^{8})$ we can compute it from diagonal (property **5**). ↵
↵
Implementation:↵
---------------↵
↵
Here is a code↵
↵
<spoiler summary="Code">↵
```cpp↵
#include <bits/stdc++.h>↵
↵
using namespace std;↵
//#define int long long↵
///You can add it if you want↵
const int md = 998244353;↵
mt19937 rnd;↵
#define app push_back↵
#define all(x) (x).begin(),(x).end()↵
#ifdef LOCAL↵
#define debug(...) [](auto...a){ ((cout << a << ' '), ...) << endl;}(#__VA_ARGS__, ":", __VA_ARGS__)↵
#define debugv(v) do {cout<< #v <<" : {"; for(int izxc=0;izxc<v.size();++izxc) {cout << v[izxc];if(izxc+1!=v.size()) cout << ","; }cout <<"}"<< endl;} while(0)↵
#else↵
#define debug(...)↵
#define debugv(v)↵
#endif↵
#define lob(a,x) lower_bound(all(a),x)↵
#define upb(a,x) upper_bound(all(a),x)↵
↵
template<int M, int K, int G>↵
struct Fft↵
{↵
// 1, 1/4, 1/8, 3/8, 1/16, 5/16, 3/16, 7/16, ...↵
int g[1 << (K - 1)];↵
↵
Fft() : g()↵
{ //if tl constexpr...↵
// static_assert(K >= 2, "Fft: K >= 2 must hold");↵
g[0] = 1;↵
g[1 << (K - 2)] = G;↵
for (int l = 1 << (K - 2); l >= 2; l >>= 1)↵
{↵
g[l >> 1] = (g[l] * 1LL * g[l]) % M;↵
}↵
assert((g[1]*1LL * g[1]) % M == M - 1);↵
for (int l = 2; l <= 1 << (K - 2); l <<= 1)↵
{↵
for (int i = 1; i < l; ++i)↵
{↵
g[l + i] = (g[l] * 1LL * g[i]) % M;↵
}↵
}↵
}↵
↵
void fft(vector<int> &x) const↵
{↵
const int n = x.size();↵
assert(n <= 1 << K);↵
for (int h = __builtin_ctz(n); h--;)↵
{↵
const int l = (1 << h);↵
for (int i = 0; i < n >> (h + 1); ++i)↵
{↵
for (int j = i << (h + 1); j < (((i << 1) + 1) << h); ++j)↵
{↵
const int t = (g[i] * 1LL * x[j | l]) % M;↵
x[j | l] = x[j] - t;↵
if (x[j | l] < 0)↵
x[j | l] += M;↵
x[j] += t;↵
if (x[j] >= M)↵
x[j] -= M;↵
}↵
}↵
}↵
for (int i = 0, j = 0; i < n; ++i)↵
{↵
if (i < j)↵
std::swap(x[i], x[j]);↵
for (int l = n; (l >>= 1) && !((j ^= l) & l);)↵
{↵
}↵
}↵
}↵
↵
vector<int> convolution(vector<int> a, vector<int> b) const↵
{↵
if (a.empty() || b.empty())↵
return {};↵
const int p =
for (int &x: a)↵
{↵
x %= p;↵
if (x >= p)↵
x -= p;↵
if (x < 0)↵
x += p;↵
}↵
for (int &x: b)↵
{↵
x %= p;↵
if (x >= p)↵
x -= p;↵
if (x < 0)↵
x += p;↵
}↵
const int na = a.size(), nb = b.size();↵
int n, invN = 1;↵
for (n = 1; n < na + nb - 1; n <<= 1)↵
invN = ((invN & 1) ? (invN + M) : invN) >> 1;↵
vector<int> x(n, 0), y(n, 0);↵
std::copy(a.begin(), a.end(), x.begin());↵
std::copy(b.begin(), b.end(), y.begin());↵
fft(x);↵
fft(y);↵
for (int i = 0; i < n; ++i)↵
x[i] = (((static_cast<long long>(x[i]) * y[i]) % M) * invN) % M;↵
std::reverse(x.begin() + 1, x.end());↵
fft(x);↵
x.resize(na + nb - 1);↵
return x;↵
}↵
};↵
↵
Fft<998244353, 21, 31 * 31 * 31 * 31> muls;↵
↵
template<int32_t MOD>↵
struct ModInt↵
{↵
int32_t value;↵
↵
ModInt() : value(0)↵
{↵
}↵
↵
ModInt(long long v) : value(v % MOD)↵
{↵
if (value < 0)↵
value += MOD;↵
}↵
↵
ModInt(int32_t v): value(v % MOD)↵
{↵
if (value < 0)↵
value += MOD;↵
}↵
↵
ModInt operator+=(ModInt m)↵
{↵
value += m.value;↵
if (value >= MOD)↵
value -= MOD;↵
return value;↵
}↵
↵
ModInt operator-=(ModInt m)↵
{↵
value -= m.value;↵
if (value < 0)↵
value += MOD;↵
return value;↵
}↵
↵
ModInt operator*=(ModInt m)↵
{↵
value = (value * 1LL * m.value) % MOD;↵
return value;↵
}↵
↵
ModInt power(long long exp) const↵
{↵
if (exp == 0)↵
return 1;↵
ModInt res = (exp & 1 ? value : 1);↵
ModInt half = power(exp >> 1);↵
return res * half * half;↵
}↵
↵
ModInt operator/=(ModInt m) { return *this *= m.power(MOD - 2); }↵
↵
friend std::istream &operator>>(std::istream &is, ModInt &m)↵
{↵
is >> m.value;↵
return is;↵
}↵
↵
friend std::ostream &operator<<(std::ostream &os, const ModInt &m)↵
{↵
os << m.value;↵
return os;↵
}↵
↵
explicit operator int32_t() const { return value; }↵
↵
explicit operator long long() const { return value; }↵
↵
static int32_t mod() { return MOD; }↵
};↵
↵
template<int32_t MOD>↵
ModInt<MOD> operator+(ModInt<MOD> a, ModInt<MOD> b) { return a += b; }↵
↵
template<int32_t MOD, typename L>↵
ModInt<MOD> operator+(L a, ModInt<MOD> b) { return ModInt<MOD>(a) += b; }↵
↵
template<int32_t MOD, typename R>↵
ModInt<MOD> operator+(ModInt<MOD> a, R b) { return a += b; }↵
↵
template<int32_t MOD>↵
ModInt<MOD> operator-(ModInt<MOD> a, ModInt<MOD> b) { return a -= b; }↵
↵
template<int32_t MOD, typename L>↵
ModInt<MOD> operator-(L a, ModInt<MOD> b) { return ModInt<MOD>(a) -= b; }↵
↵
template<int32_t MOD, typename R>↵
ModInt<MOD> operator-(ModInt<MOD> a, R b) { return a -= b; }↵
↵
template<int32_t MOD>↵
ModInt<MOD> operator*(ModInt<MOD> a, ModInt<MOD> b) { return a *= b; }↵
↵
template<int32_t MOD, typename L>↵
ModInt<MOD> operator*(L a, ModInt<MOD> b) { return ModInt<MOD>(a) *= b; }↵
↵
template<int32_t MOD, typename R>↵
ModInt<MOD> operator*(ModInt<MOD> a, R b) { return a *= b; }↵
↵
template<int32_t MOD>↵
ModInt<MOD> operator/(ModInt<MOD> a, ModInt<MOD> b) { return a /= b; }↵
↵
template<int32_t MOD, typename L>↵
ModInt<MOD> operator/(L a, ModInt<MOD> b) { return ModInt<MOD>(a) /= b; }↵
↵
template<int32_t MOD, typename R>↵
ModInt<MOD> operator/(ModInt<MOD> a, R b) { return a /= b; }↵
↵
template<int32_t MOD>↵
bool operator==(ModInt<MOD> a, ModInt<MOD> b) { return a.value == b.value; }↵
↵
template<int32_t MOD, typename L>↵
bool operator==(L a, ModInt<MOD> b) { return a == b.value; }↵
↵
template<int32_t MOD, typename R>↵
bool operator==(ModInt<MOD> a, R b) { return a.value == b; }↵
↵
template<int32_t MOD>↵
bool operator!=(ModInt<MOD> a, ModInt<MOD> b) { return a.value != b.value; }↵
↵
template<int32_t MOD, typename L>↵
bool operator!=(L a, ModInt<MOD> b) { return a != b.value; }↵
↵
template<int32_t MOD, typename R>↵
bool operator!=(ModInt<MOD> a, R b) { return a.value != b; }↵
↵
using mint = ModInt<md>;↵
mint inv(mint x) { return 1 / x; }↵
↵
__int128 gcd(__int128 a, __int128 b, __int128 &x, __int128 &y)↵
{↵
if (b == 0)↵
{↵
x = 1;↵
y = 0;↵
return a;↵
}↵
__int128 d = gcd(b, a % b, y, x);↵
y -= a / b * x;↵
return d;↵
}↵
↵
__int128 inv(__int128 r, __int128 m)↵
{↵
__int128 x, y;↵
gcd(r, m, x, y);↵
return (x + m) % m;↵
}↵
↵
__int128 crt(__int128 r, __int128 n, __int128 c, __int128 m) { return r + ((c - r) % m + m) * inv(n, m) % m * n; }↵
const int m2 = 167772161, m3 = 469762049;↵
Fft<m2, 21, 147771621> muls2;↵
Fft<m3, 21, 297449090> muls3;↵
↵
vector<mint> operator*(vector<mint> a, vector<mint> b)↵
{ ///modulo-dependent convolution↵
if (a.empty() || b.empty())↵
return {};↵
if (md == 998244353)↵
{↵
vector<int> a1(a.size());↵
for (int i = 0; i < a.size(); ++i)↵
a1[i] = a[i].value;↵
vector<int> b1(b.size());↵
for (int i = 0; i < b.size(); ++i)↵
b1[i] = b[i].value;↵
vector<int> c1 = muls.convolution(a1, b1);↵
vector<mint> c;↵
for (int x: c1)↵
c.app(x);↵
return c;↵
}↵
else↵
{↵
vector<int> a1(a.size());↵
for (int i = 0; i < a.size(); ++i)↵
a1[i] = a[i].value;↵
vector<int> b1(b.size());↵
for (int i = 0; i < b.size(); ++i)↵
b1[i] = b[i].value;↵
vector<int> c1 = muls.convolution(a1, b1);↵
vector<int> c2 = muls2.convolution(a1, b1);↵
vector<int> c3 = muls3.convolution(a1, b1);↵
assert(c1.size()==c2.size() && c2.size()==c3.size());↵
vector<int> c4(c1.size());↵
for (int i = 0; i < c1.size(); ++i)↵
{↵
__int128 ost1 = c1[i];↵
__int128 m1 = 998244353;↵
__int128 ost2 = c2[i];↵
__int128 ost3 = c3[i];↵
__int128 ost = crt(crt(ost1, m1, ost2, m2), m1 * 1LL * m2, ost3, m3);↵
c4[i] = (ost % md);↵
}↵
vector<mint> c;↵
for (int x: c4)↵
c.app(x);↵
return c;↵
}↵
}↵
↵
vector<vector<mint> > gaussbasis(vector<vector<mint> > A) ///returns basis of Av=0↵
{↵
int n = A.size();↵
int m = A[0].size();↵
int bi = 0;↵
for (int i = 0; i < n; ++i)↵
{↵
if (bi == m)↵
break;↵
for (int j = i; j < n; ++j)↵
{↵
if (A[j][bi] != 0)↵
{↵
if (j != i) { swap(A[i], A[j]); }↵
break;↵
}↵
}↵
if (A[i][bi] != 0)↵
{↵
mint o = 1 / A[i][bi];↵
for (int j = i + 1; j < n; ++j)↵
{↵
mint we = (A[j][bi] * o);↵
for (int k = bi; k < m; ++k)↵
{↵
A[j][k] -= we * A[i][k];↵
}↵
}↵
}↵
else↵
{↵
++bi;↵
--i;↵
continue;↵
}↵
}↵
vector<int> indices(m);↵
iota(all(indices), 0);↵
for (int i = n - 1; i >= 0; --i)↵
{↵
int bi = 0;↵
while (bi < m && A[i][bi] == 0) { ++bi; }↵
if (bi < m)↵
{↵
indices.erase(find(all(indices), bi));↵
}↵
}↵
vector<vector<mint> > v(indices.size(), vector<mint>(m, 0));↵
for (int i = 0; i < indices.size(); ++i)↵
{↵
v[i][indices[i]] = 1;↵
}↵
for (int i = n - 1; i >= 0; --i)↵
{↵
int bi = 0;↵
while (bi < m && A[i][bi] == 0) { ++bi; }↵
if (bi == m)↵
continue;↵
for (int k = 0; k < indices.size(); ++k)↵
{↵
mint cur = 0;↵
for (int j = bi + 1; j < m; ++j)↵
{↵
cur -= A[i][j] * v[k][j];↵
}↵
v[k][bi] = cur / A[i][bi];↵
}↵
}↵
return v;↵
}↵
↵
optional<vector<mint> > gauss(vector<vector<mint> > A, vector<mint> b) ///returns v such that Av=b↵
{↵
int n = A.size();↵
assert(b.size()==n);↵
int m = A[0].size();↵
int bi = 0;↵
for (int i = 0; i < n; ++i)↵
{↵
if (bi == m)↵
break;↵
for (int j = i; j < n; ++j)↵
{↵
if (A[j][bi] != 0)↵
{↵
if (j != i)↵
{↵
swap(A[i], A[j]);↵
swap(b[i], b[j]);↵
}↵
break;↵
}↵
}↵
if (A[i][bi] != 0)↵
{↵
mint o = inv(A[i][bi]);↵
for (int j = i + 1; j < n; ++j)↵
{↵
mint we = (A[j][bi] * o);↵
b[j] -= we * b[i];↵
for (int k = bi; k < m; ++k)↵
{↵
A[j][k] -= we * A[i][k];↵
}↵
}↵
}↵
else↵
{↵
++bi;↵
--i;↵
continue;↵
}↵
}↵
vector<mint> v(m);↵
for (int i = n - 1; i >= 0; --i)↵
{↵
int bi = 0;↵
while (bi < m && A[i][bi] == 0) { ++bi; }↵
if (bi == m)↵
{↵
if (b[i] != 0) { return nullopt; }↵
else { continue; }↵
} {↵
mint cur = b[i];↵
for (int j = bi + 1; j < m; ++j)↵
{↵
cur -= A[i][j] * v[j];↵
}↵
v[bi] = cur * inv(A[i][bi]);↵
}↵
}↵
return v;↵
}↵
↵
optional<vector<vector<mint> > > findPrecursion(vector<mint> a)↵
{ ///finds P-recursion of a given sequence A by gauss↵
for (int snd = 0; snd <= 20; ++snd)↵
{↵
for (int n = 1; n <= snd - 1; ++n)↵
{↵
vector<vector<mint> > A;↵
int d = snd - n;↵
int eq = ((int) (a.size())) - (n - 1);↵
if (eq < n * d) { continue; }↵
for (int i = n - 1; i < a.size(); ++i)↵
{↵
vector<mint> u;↵
for (int j = 0; j < n; ++j)↵
{↵
mint de = 1;↵
for (int k = 0; k < d; ++k)↵
{↵
u.app(a[i - j] * de);↵
de *= i;↵
}↵
}↵
A.app(u);↵
}↵
vector<vector<mint> > zx = gaussbasis(A);↵
if (zx.empty())↵
continue;↵
//debug(n, d);↵
vector<mint> ans = zx[0];↵
vector<vector<mint> > res;↵
for (int j = 0; j < n; ++j)↵
{↵
res.app({});↵
for (int k = 0; k < d; ++k)↵
{↵
res[j].app(ans[j * d + k]);↵
}↵
}↵
return res;↵
}↵
}↵
return nullopt;↵
}↵
↵
optional<vector<mint> > evaluatePrecursion(vector<mint> a, vector<vector<mint> > rec, int sz)↵
{ ///a(0),...,a(a.size()-1) -> (by P-recursion rec) a(0),...,a(sz-1)↵
int n = rec.size();↵
int d = rec[0].size();↵
int given = a.size();↵
if (given >= sz)↵
{↵
a.resize(sz);↵
return a;↵
}↵
if (a.size() < n) { return nullopt; }↵
vector<mint> tore;↵
for (int i = given; i < sz; ++i)↵
{↵
mint de = 1;↵
mint s = 0;↵
for (int k = 0; k < d; ++k)↵
{↵
s += de * rec[0][k];↵
de *= i;↵
}↵
if (s == 0) { return nullopt; }↵
tore.app(s);↵
}↵
vector<mint> pref(tore.size() + 1);↵
pref[0] = 1;↵
for (int i = 0; i < tore.size(); ++i) { pref[i + 1] = pref[i] * tore[i]; }↵
mint pro = pref[tore.size()];↵
mint invpro = 1 / pro;↵
mint cur = invpro;↵
vector<mint> invtore(tore.size());↵
for (int i = tore.size() - 1; i >= 0; --i)↵
{↵
invtore[i] = cur * pref[i];↵
cur *= tore[i];↵
}↵
for (int i = given; i < sz; ++i)↵
{↵
mint chi = 0;↵
for (int j = 1; j < n; ++j)↵
{↵
mint de = 1;↵
for (int k = 0; k < d; ++k)↵
{↵
chi += a[i - j] * de * rec[j][k];↵
de *= i;↵
}↵
}↵
a.app(((mint) (0)) - chi * invtore[i - given]);↵
}↵
return a;↵
}↵
↵
mint value(vector<mint> a, mint x)↵
{ ///A(x)↵
mint de = 1;↵
mint ans = 0;↵
for (int i = 0; i < a.size(); ++i)↵
{↵
ans += a[i] * de;↵
de *= x;↵
}↵
return ans;↵
}↵
↵
vector<mint> shiftofsamplingpoints(vector<mint> a)↵
{ ///P(0),...,P(t) we want to compute P(0),...,P(4t+1)↵
int t = a.size() - 1;↵
vector<mint> fact(4 * t + 2);↵
fact[0] = 1;↵
for (int i = 1; i < 4 * t + 2; ++i)↵
fact[i] = fact[i - 1] * i;↵
vector<mint> invf(4 * t + 2);↵
invf[4 * t + 1] = 1 / fact[4 * t + 1];↵
for (int i = 4 * t; i >= 0; --i) { invf[i] = (invf[i + 1] * (i + 1)); }↵
assert(invf[0]==1);↵
vector<mint> invm(4 * t + 2, 0);↵
for (int i = 1; i < 4 * t + 2; ++i) { invm[i] = fact[i - 1] * invf[i]; }↵
vector<mint> values(t + 1, 0);↵
for (int k = 0; k <= t; ++k)↵
{↵
mint o = 1;↵
if ((t - k) % 2 == 1) { o = (((mint) (0)) - 1); }↵
values[k] = (a[k] * invf[k] * invf[t - k] * o);↵
}↵
vector<mint> h = invm * values;↵
vector<mint> res;↵
for (int i = 0; i <= t; ++i)↵
{↵
res.app(a[i]);↵
}↵
for (int x = t + 1; x <= 4 * t + 1; ++x)↵
{↵
mint ans = fact[x];↵
ans *= invf[x - t - 1];↵
ans *= h[x];↵
res.app(ans);↵
}↵
return res;↵
}↵
↵
optional<mint> evaluatePrecursionfast(vector<mint> a, vector<vector<mint> > rec, int id)↵
{ ///a(0),...,a(a.size()-1) -> (by P-recursion rec) a(id), O(sqrt(id)*log(id))↵
if (id < a.size())↵
return a[id];↵
int n = rec.size();↵
int d = 1;↵
for (auto &v: rec) { d = max(d, ((int) (v.size() - 1))); }↵
if (a.size() < n - 1) { return nullopt; }↵
if (n == 1) { return 0; }↵
int l = n - 1;↵
int shift = 0;↵
while (a.size() > l)↵
{↵
a.erase(a.begin());↵
++shift;↵
--id;↵
}↵
int u = 0;↵
while ((1LL << u) * (1LL << u) <= id) { ++u; }↵
int B = (1 << u);↵
vector<mint> S;↵
vector<vector<vector<mint> > > A(l, vector<vector<mint> >(l));↵
int sz = d;↵
for (int k = 0; k <= d; ++k)↵
{↵
S.app(value(rec[0], k + l + shift));↵
}↵
for (int i = 0; i < l - 1; ++i)↵
{↵
for (int j = 0; j < l; ++j)↵
{↵
if (j == i + 1)↵
{↵
for (int k = 0; k <= d; ++k)↵
{↵
A[i][j].app(value(rec[0], k + l + shift));↵
}↵
}↵
else↵
{↵
for (int k = 0; k <= d; ++k)↵
{↵
A[i][j].app(0);↵
}↵
}↵
}↵
}↵
for (int j = 0; j < l; ++j)↵
{↵
for (int k = 0; k <= d; ++k)↵
{↵
A[l - 1][j].app(((mint) (0)) - value(rec[l - j], k + l + shift));↵
}↵
}↵
for (int s = 0; s < u; ++s)↵
{↵
S = shiftofsamplingpoints(S);↵
assert(S.size()==4*sz+2);↵
for (int i = 0; i < l; ++i)↵
{↵
for (int j = 0; j < l; ++j)↵
{↵
A[i][j] = shiftofsamplingpoints(A[i][j]);↵
assert(A[i][j].size()==4*sz+2);↵
}↵
}↵
vector<vector<vector<mint> > > newA(l, vector<vector<mint> >(l, vector<mint>(2 * sz + 1, 0)));↵
for (int k = 0; k <= 2 * sz; ++k)↵
{↵
for (int ii = 0; ii < l; ++ii)↵
{↵
for (int jj = 0; jj < l; ++jj)↵
{↵
for (int kk = 0; kk < l; ++kk)↵
{↵
newA[ii][kk][k] += A[ii][jj][2 * k + 1] * A[jj][kk][2 * k];↵
}↵
}↵
}↵
}↵
vector<mint> newS(2 * sz + 1, 0);↵
for (int k = 0; k <= 2 * sz; ++k) { newS[k] = S[2 * k] * S[2 * k + 1]; }↵
sz *= 2;↵
S = newS;↵
A = newA;↵
}↵
int k = (id - l) / B;↵
assert(k>=0); ///id>=l+1 here↵
mint pro = 1;↵
vector<mint> v;↵
for (int i = 0; i < l; ++i)↵
v.app(a[i]);↵
for (int i = 0; i < k; ++i)↵
{↵
vector<mint> newv(l, 0);↵
vector<vector<mint> > M(l, vector<mint>(l, 0));↵
for (int ii = 0; ii < l; ++ii)↵
{↵
for (int jj = 0; jj < l; ++jj)↵
{↵
assert(i<A[ii][jj].size());↵
M[ii][jj] = A[ii][jj][i];↵
}↵
}↵
pro *= S[i];↵
for (int ii = 0; ii < l; ++ii)↵
{↵
for (int jj = 0; jj < l; ++jj)↵
{↵
newv[ii] += M[ii][jj] * v[jj];↵
}↵
}↵
v = newv;↵
}↵
int cur = k * B + l - 1;↵
assert(cur<id);↵
while (cur < id)↵
{↵
mint newval = 0;↵
for (int i = 0; i < l; ++i)↵
{↵
newval -= v[i] * value(rec[l - i], cur + 1 + shift);↵
}↵
mint de = value(rec[0], cur + 1 + shift);↵
pro *= de;↵
for (mint &x: v) { x *= de; }↵
v.erase(v.begin());↵
v.app(newval);↵
++cur;↵
}↵
if (pro == 0)↵
return nullopt;↵
return v.back() / pro;↵
}↵
↵
optional<vector<mint> > sequenceextender(vector<mint> a, int sz)↵
{ ///finds P-recursion, and if was found, calculates a(0),...,a(sz-1)↵
if (a.size() >= sz)↵
{↵
a.resize(sz);↵
return a;↵
}↵
auto uu = findPrecursion(a);↵
if (!uu)↵
return nullopt;↵
auto rec = (*uu);↵
auto ans = evaluatePrecursion(a, rec, sz);↵
if (!ans)↵
return nullopt;↵
return (*ans);↵
}↵
↵
optional<mint> fastgetvaluebyid(vector<mint> a, int id)↵
{ ///finds P-recursion, and if was found, calculates a(id) in O(sqrt(id)*log(id))↵
if (a.size() > id) { return a[id]; }↵
auto uu = findPrecursion(a);↵
if (!uu)↵
return nullopt;↵
auto rec = (*uu);↵
auto ans = evaluatePrecursionfast(a, rec, id);↵
if (!ans)↵
return nullopt;↵
return (*ans);↵
}↵
↵
optional<mint> optimalgetvaluebyid(vector<mint> a, int id)↵
{ ///finds P-recursion, and if was found, calculates a(id) by choosing optimal of O(id) method and O(sqrt(id)*log(id)) method↵
if (a.size() > id) { return a[id]; }↵
auto uu = findPrecursion(a);↵
if (!uu)↵
return nullopt;↵
auto rec = (*uu);↵
int n = rec.size();↵
int d = rec[0].size();↵
double C = 1;↵
if (md != 998244353)↵
C = 3;↵
double val1 = sqrt(id) * log(id) * C * n * n * d + sqrt(id) * n * n * n * d;↵
double val2 = n * 1.0 * d * 1.0 * id;↵
//debug(val1, val2);↵
if (val1 < val2)↵
{↵
//debug("fastgetvalue");↵
auto ans = evaluatePrecursionfast(a, rec, id);↵
if (!ans)↵
return nullopt;↵
return (*ans);↵
}↵
else↵
{↵
//debug("sequenceextender");↵
auto ans = evaluatePrecursion(a, rec, id + 1);↵
if (!ans)↵
return nullopt;↵
return (*ans)[id];↵
}↵
}↵
↵
vector<vector<mint> > transpose(vector<vector<mint> > a)↵
{ ///transposes the table A↵
if (a.empty())↵
return a;↵
int n = a.size();↵
int m = a[0].size();↵
vector<vector<mint> > b(m, vector<mint>(n, 0));↵
for (int i = 0; i < n; ++i)↵
for (int j = 0; j < m; ++j)↵
b[j][i] = a[i][j];↵
return b;↵
}↵
↵
/// If the size of a is big TL (too many Gauss), if the size of a is small WA (not enough for finding the P-recursion), you should keep balance↵
optional<vector<mint> > extendtableslow(vector<vector<mint> > a, vector<pair<int, int> > que)↵
{ /// extends table A, finding A(que[i].first,que[i].second), in a O(sum(que[i]))↵
if (que.empty())↵
{↵
vector<mint> ret = {};↵
return ret;↵
}↵
int ma = 0;↵
for (auto [i,j]: que) { ma = max(ma, j); }↵
int n = a.size();↵
int m = a[0].size();↵
vector<vector<mint> > ex(n);↵
for (int i = 0; i < n; ++i)↵
{↵
auto h = sequenceextender(a[i], ma + 1);↵
if (!h) { return nullopt; }↵
ex[i] = (*h);↵
}↵
vector<mint> res;↵
for (auto [i,j]: que)↵
{↵
vector<mint> e;↵
for (int k = 0; k < n; ++k)↵
{↵
e.app(ex[k][j]);↵
}↵
auto h = sequenceextender(e, i + 1);↵
if (!h) { return nullopt; }↵
res.app((*h)[i]);↵
}↵
return res;↵
}↵
↵
optional<vector<mint> > extendtablefast(vector<vector<mint> > a, vector<pair<int, int> > que)↵
{ /// extends table A, finding A(que[i].first,que[i].second)↵
if (que.empty())↵
{↵
vector<mint> ret = {};↵
return ret;↵
}↵
int n = a.size();↵
int m = a[0].size();↵
vector<vector<mint> > ex(n);↵
vector<vector<mint> > rec[n];↵
for (int k = 0; k < n; ++k)↵
{↵
auto uu = findPrecursion(a[k]);↵
if (!uu) { return nullopt; }↵
rec[k] = (*uu);↵
}↵
vector<mint> res;↵
for (auto [i,j]: que)↵
{↵
vector<mint> e;↵
for (int k = 0; k < n; ++k)↵
{↵
auto uu = evaluatePrecursionfast(a[k], rec[k], j);↵
if (!uu) { return nullopt; }↵
e.app(*uu);↵
}↵
auto h = optimalgetvaluebyid(e, i);↵
if (!h) { return nullopt; }↵
res.app(*h);↵
}↵
return res;↵
}↵
↵
optional<vector<mint> > extendtable(vector<vector<mint> > a, vector<pair<int, int> > que)↵
{ /// extends table A, finding A(que[i].first,que[i].second)↵
double C = 1;↵
if (md != 998244353) { C = 3; }↵
double op1 = 0;↵
double op2 = 0;↵
for (auto [i,j]: que)↵
{↵
op1 += C * 1.0 * ((int) (a.size())) * 1.0 * sqrt(i + 1) * 1.0 * log(i + 2) * 5 * 5 * 5;↵
op1 += C * sqrt(j + 1) * log(j + 2) * 5 * 5 * 5;↵
}↵
for (auto [i,j]: que)↵
{↵
op2 += C * 1.0 * ((int) (a.size())) * 1.0 * (i + 1) * 1.0 * 5 * 5;↵
op2 += C * 1.0 * (j + 1) * 1.0 * 5 * 5;↵
}↵
//debug(op1, op2);↵
if (op1 < op2) { return extendtablefast(a, que); }↵
else { return extendtableslow(a, que); }↵
}↵
↵
optional<vector<mint> > getcolumnoftable(vector<vector<mint> > a, int col, int size)↵
{ /// get A(0,col),A(1,col),...,A(size-1,col)↵
int n = a.size();↵
int m = a[0].size();↵
vector<mint> e;↵
for (int i = 0; i < n; ++i)↵
{↵
auto h = sequenceextender(a[i], col + 1);↵
if (!h) { return nullopt; }↵
e.app((*h)[col]);↵
}↵
auto h = sequenceextender(e, size);↵
if (!h)↵
return nullopt;↵
return *h;↵
}↵
↵
optional<vector<mint> > getrowoftable(vector<vector<mint> > a, int row, int size)↵
{ /// get A(row,0),A(row,1),...,A(row,size-1)↵
return getcolumnoftable(transpose(a), row, size);↵
}↵
↵
void test()↵
{↵
vector<mint> v1 = {1, 1, 3, 7, 19, 51, 141, 393, 1107, 3139}; ///(1+x+x^2)^n [x^n]↵
auto uu1 = sequenceextender(v1, 30);↵
if (uu1)↵
{↵
auto u1 = (*uu1);↵
debugv(u1);↵
}↵
vector<mint> v2 = {1, 30, 465, 4930, 40020, 264306, 1474795, 7133130, 30462615, 116470380}; ///(1+x+x^2)^30 [x^n]↵
auto uu2 = sequenceextender(v2, 70);↵
if (uu2)↵
{↵
auto u2 = (*uu2);↵
debugv(u2);↵
}↵
vector<mint> v3 = {0, 567646151, 513265721, 604121291, 715018514, 398975714, 610803800, 499563577, 491416403,↵
913506524↵
}; ///s(30,n), not D-finite↵
auto uu3 = sequenceextender(v3, 35);↵
if (uu3)↵
{↵
auto u3 = (*uu3);↵
debugv(u3);↵
}↵
vector<mint> v4 = {0, 1, 2, 9, 44, 265, 1854, 14833, 133496, 1334961, 14684570};↵
///n!*x*e^(-x) [x^n] (number of permutations of size n without stable points)↵
auto uu4 = sequenceextender(v4, 20);↵
if (uu4)↵
{↵
auto u4 = (*uu4);↵
debugv(u4);↵
}↵
auto uu5 = fastgetvaluebyid({1, 1, 2, 6, 24, 120, 720}, 11); ///factorials ,by id↵
if (uu5)↵
{↵
auto u5 = (*uu5);↵
debug(u5);↵
}↵
auto uu6 = fastgetvaluebyid(v1, 29); ///(1+x+x^2)^n [x^n], by id↵
if (uu6)↵
{↵
auto u6 = (*uu6);↵
debug(u6);↵
}↵
auto uu7 = fastgetvaluebyid(v1, 500000000); ///(1+x+x^2)^n [x^n], by id↵
if (uu7)↵
{↵
auto u7 = (*uu7);↵
debug(u7);↵
}↵
auto uu8 = optimalgetvaluebyid(v1, 500000000); ///(1+x+x^2)^n [x^n], by id↵
if (uu8)↵
{↵
auto u8 = (*uu8);↵
debug(u8);↵
}↵
auto uu9 = optimalgetvaluebyid({1, 1, 2, 6, 24, 120, 720}, 11); ///factorials ,by id↵
if (uu9)↵
{↵
auto u9 = (*uu9);↵
debug(u9);↵
}↵
auto uu10 = optimalgetvaluebyid({1, 1, 2, 6, 24, 120, 720}, 998244352); ///factorials ,by id↵
if (uu10)↵
{↵
auto u10 = (*uu10);↵
debug(u10);↵
}↵
vector<vector<mint> > table1 =↵
{{1, 2, 3, 4, 5}, {2, 3, 4, 5, 6}, {3, 4, 5, 6, 7}, {4, 5, 6, 7, 8}, {5, 6, 7, 8, 9}}; ///f(i,j)=i+j+1↵
vector<pair<int, int> > que1 = {{1, 1}, {0, 0}, {5, 7}, {11342, 1333}};↵
auto uu11 = extendtablefast(table1, que1);↵
auto uu12 = extendtableslow(table1, que1);↵
auto uu13 = extendtable(table1, que1);↵
debug((bool) (uu11));↵
debug((bool) (uu12));↵
debug((bool) (uu13));↵
if (uu11 && uu12 && uu13)↵
{↵
auto u11 = (*uu11);↵
auto u12 = (*uu12);↵
auto u13 = (*uu13);↵
debugv(u11);↵
debugv(u12);↵
debugv(u13);↵
}↵
auto uu14 = getcolumnoftable(table1, 10, 100);↵
if (uu14)↵
{↵
auto u14 = (*uu14);↵
debugv(u14);↵
}↵
auto uu15 = getrowoftable(table1, 50, 70);↵
if (uu15)↵
{↵
auto u15 = (*uu15);↵
debugv(u15);↵
}↵
///exp((x+y)/((1-x)(1-y)), I removed it from examples↵
}↵
↵
int32_t main()↵
{↵
ios_base::sync_with_stdio(false);↵
cin.tie(0);↵
test();↵
return 0;↵
}↵
```↵
</spoiler>↵
↵
Code works for all prime moduli. Explanation:↵
↵
* `findPrecursion` finds P-recursion ↵
↵
* `evaluatePrecursion` evaluates P-recursion from the initial values if possible ↵
↵
* `evaluatePrecursionfast` evaluates the element of a 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` and `getrowoftable` get column or row of an extended table ↵
↵
↵
A "small" problem↵
-----------------↵
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)+\ldots+a_{n-l}P_{l}(n)=0$.↵
↵
You know $a(0),\ldots,a(l-1)$ and you want to extend this sequence to $a(0),\ldots,a(N)$. But what if $P_{0}(k)=0$ for some $l \leq k \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 functions when we can use it, it can't be zero. ↵
↵
Example where it matters↵
------------------------↵
↵
Let's try to compute $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. ↵
↵
Examples of the problems↵
------------------------↵
↵
**1)** [Atcoder Regular 202 problem D](https://atcoder.jp/contests/arc202/tasks/arc202_d)↵
↵
<spoiler>↵
Here if you are solving it in $O((H+w+T)^{1.5})$ using the reflection principle you have a subproblem to calculate $(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,\ldots,k+C$, and then extend it to the whole sequence in linear time.↵
</spoiler>↵
↵
**2)** We can also solve [problem: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: [submission:330486092]).↵
↵
I think there are some other examples of problems that can be solved by this method. ↵
↵
Conclusion (Warning)↵
====================↵
↵
When I got it at first, it seemed to me that "every" function in combinatorics is D-finite, but it is actually wrong (see before). I think that one should use this algorithm only if one has the proof that some function is D-finite, otherwise it can be a waste of time. So I don't recommend using this algorithm as magic.↵
↵
↵
↵
↵
Thanks to [user:orz,2025-07-25] for editing this blog!
