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' \lt 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 \lt k,t \lt 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).

$$$(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)}$$$.

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

.

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))$$$.

#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 = M;
        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;
}

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.