题目简述

给定 $n, p, k$，其中 $p$ 是素数，求在 $\mathbb{F}_p$ 上的 $n \times n$ 矩阵中，秩恰好等于 $r$ 的矩阵个数，对 $998244353$ 取模，输出 $r=0,1,\dots,k$ 的答案。
约束：$1 \le n \le 10^{18}$，$2 \le p < 998244353$，$0 \le k \le 5000$。

---

问题转化

设 $dp_{m,k}$ 表示 $\mathbb{F}_p$ 上 $m \times n$ 矩阵中秩为 $k$ 的矩阵个数。
考虑从 $m-1$ 行矩阵 $B$ 扩展到最后一行 $C$，有：

• 若 $\operatorname{rank}(B) = \operatorname{rank}(A) = k$，则 $C$ 必须落在 $B$ 的行空间内，该空间大小为 $p^k$。
• 若 $\operatorname{rank}(B) = k-1$，$\operatorname{rank}(A)=k$，则 $C$ 不能落在 $B$ 的行空间内，但必须与 $B$ 的行空间一起张成 $k$ 维空间。这样的 $C$ 有 $p^n - p^{k-1}$ 个。

于是递推式：

$$dp_{m,k} = dp_{m-1,k}\cdot p^k + dp_{m-1,k-1}\cdot (p^n - p^{k-1}). $$

该递推中 $p^n - p^{k-1}$ 项与 $m$ 无关但依赖于 $n$，处理不便。定义

$$g_{m,k} = \frac{dp_{m,k}}{\prod_{i=0}^{k-1} (p^n - p^i)}. $$

代入可得 $g$ 满足：

$$g_{m,k} = g_{m-1,k}\cdot p^k + g_{m-1,k-1},$$

初始 $g_{0,0}=1$，且当 $k<0$ 或 $k>m$ 时 $g_{m,k}=0$。

于是最终答案：

$$\text{ans}_r = g_{n,r} \cdot \prod_{i=0}^{r-1} (p^n - p^i) \pmod{998244353}. $$

因此问题归结为：给定 $n$（极大），计算 $g_{n,0}, g_{n,1}, \dots, g_{n,k}$，其中 $g$ 满足上述递推。

---

计算 $g_{n,k}$

1. 递推的矩阵形式与快速幂

将 $g_m = (g_{m,0}, g_{m,1}, \dots, g_{m,k})$ 视为行向量，递推可写作：

$$g_{m+1} = g_m \cdot T,$$

其中 $T$ 是上三角矩阵：

$$T_{i,i} = p^i,\quad T_{i,i+1} = 1,\quad \text{其余为 }0. $$

故 $g_n = g_0 \cdot T^n$。
直接计算 $T^n$ 需要 $O(k^3 \log n)$，不可行。但注意到 $T$ 的特殊结构，我们可以利用倍增法和卷积优化。

2. 从 $g_m$ 计算 $g_{2m}$

观察递推的组合意义：$g_{m,k}$ 是从 $(0,0)$ 走到 $(m,k)$ 的路径数，每一步 $(i,j) \to (i+1, j)$ 权值为 $p^j$，$(i,j) \to (i+1, j+1)$ 权值为 $1$。
因此 $g_{m,k}$ 的生成函数为：

$$\sum_{k} g_{m,k} x^k = \prod_{i=0}^{m-1} (1 + p^i x). $$

于是：

$$g_{2m,k} = [x^k] \prod_{i=0}^{2m-1} (1 + p^i x) = [x^k] \left( \prod_{i=0}^{m-1} (1 + p^i x) \right) \left( \prod_{i=m}^{2m-1} (1 + p^i x) \right). $$

第二个乘积中提取 $p^{m}$：

$$\prod_{i=m}^{2m-1} (1 + p^i x) = \prod_{j=0}^{m-1} (1 + p^{m+j} x) = \prod_{j=0}^{m-1} (1 + p^m \cdot p^j x). $$

记 $q = p^m$，则第二个乘积为 $\prod_{j=0}^{m-1} (1 + q \cdot p^j x) = A(qx)$，其中 $A(x)=\prod_{i=0}^{m-1}(1+p^i x)$。
所以

$$\sum_{k} g_{2m,k} x^k = A(x) \cdot A(qx) = \left( \sum_{a} g_{m,a} x^a \right) \left( \sum_{b} g_{m,b} q^b x^b \right) = \sum_{k} \left( \sum_{a+b=k} g_{m,a} g_{m,b} q^b \right) x^k. $$

因此：

$$g_{2m,k} = \sum_{a+b=k} g_{m,a} g_{m,b} \cdot (p^m)^b. $$

3. 利用二次型变换为标准卷积

上式中的 $p^{m b}$ 项阻碍了直接卷积。利用恒等式：

$$m b = \binom{m+b}{2} - \binom{m}{2} - \binom{b}{2} \quad \text{和} \quad a(m-b) = am + \binom{a}{2} + \binom{b}{2} - \binom{a+b}{2}. $$

题解方法一给出的正确公式为：

$$g_{2m,k} = \sum_{a+b=k} g_{m,a} g_{m,b} \cdot p^{a(m-b)}. $$

（可以从对称性推导，但直接使用原递推更容易验证。）
利用 $a(m-b) = am + \binom{a}{2} + \binom{b}{2} - \binom{k}{2}$，得：

$$g_{2m,k} = p^{-\binom{k}{2}} \sum_{a+b=k} \left( g_{m,a} p^{am + \binom{a}{2}} \right) \left( g_{m,b} p^{\binom{b}{2}} \right). $$

因此定义：

$$A_i = g_{m,i} \cdot p^{m i} \cdot p^{\binom{i}{2}}, \quad B_i = g_{m,i} \cdot p^{\binom{i}{2}}, $$

计算卷积 $C = A * B$，则：

$$g_{2m,k} = p^{-\binom{k}{2}} \cdot C_k.$$

4. 倍增法计算 $g_n$

利用上述平方操作和“加一行”操作（$g_{m+1,k} = p^k g_{m,k} + g_{m,k-1}$），我们可以用类似快速幂的方法，将 $n$ 二进制分解，从 $g_0$ 开始，依次应用平方和加一操作，得到 $g_n$。

• 平方：$g \gets \text{square}(g)$，同时 $m \gets 2m$。
• 加一：g \gets \text{add_one}(g)，同时 $m \gets m+1$。

时间复杂度：每次平方需要一次卷积，长度为 $O(k)$，用 NTT 为 $O(k \log k)$。总共 $O(\log n)$ 次平方和 $O(\log n)$ 次加一，故总复杂度 $O(k \log k \log n)$，在 $k\le 5000$ 时可行。

5. 模运算细节

• 模数 $998244353$ 是 NTT 友好模数，原根为 $3$。
• 需要计算 $p^e \bmod MOD$，其中指数 $e$ 可能很大（如 $n$ 或 $m$）。利用费马小定理：$p^{MOD-1} \equiv 1$，所以指数可以对 $MOD-1$ 取模。注意 $p$ 与 $MOD$ 互质（$p<MOD$ 且 $MOD$ 是素数），因此可行。
• 在平方操作中需要 $p^{m i}$，其中 $m$ 是当前行数，可能很大（但 $m \le n \le 10^{18}$）。我们维护 cur_m_mod = m % (MOD-1)，这样 base = p^{cur_m_mod}，然后 base_pow[i] = base^i。
• 预处理 $p^{\binom{i}{2}}$ 及其逆元，以及 $p^i$。

---

最终答案计算

得到 $g_{n,r}$ 后，答案：

$$\text{ans}_r = g_{n,r} \cdot \prod_{i=0}^{r-1} (p^n - p^i) \pmod{998244353}. $$

其中 $p^n$ 也通过快速幂计算，指数 $n$ 对 $MOD-1$ 取模。

---

C++ 代码实现

#include <bits/stdc++.h>
using namespace std;

const int MOD = 998244353;
const int G = 3;  // primitive root

long long mod_pow(long long a, long long e) {
    long long res = 1;
    a %= MOD;
    while (e) {
        if (e & 1) res = res * a % MOD;
        a = a * a % MOD;
        e >>= 1;
    }
    return res;
}

void ntt(vector<int>& a, bool invert) {
    int n = a.size();
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1) j ^= bit;
        j ^= bit;
        if (i < j) swap(a[i], a[j]);
    }
    for (int len = 2; len <= n; len <<= 1) {
        int wlen = mod_pow(G, (MOD - 1) / len);
        if (invert) wlen = mod_pow(wlen, MOD - 2);
        for (int i = 0; i < n; i += len) {
            int w = 1;
            for (int j = 0; j < len / 2; j++) {
                int u = a[i + j];
                int v = (long long)a[i + j + len / 2] * w % MOD;
                a[i + j] = (u + v) % MOD;
                a[i + j + len / 2] = (u - v + MOD) % MOD;
                w = (long long)w * wlen % MOD;
            }
        }
    }
    if (invert) {
        int inv_n = mod_pow(n, MOD - 2);
        for (int& x : a) x = (long long)x * inv_n % MOD;
    }
}

vector<int> convolution(const vector<int>& a, const vector<int>& b) {
    int n = 1;
    int sz = a.size() + b.size() - 1;
    while (n < sz) n <<= 1;
    vector<int> fa(n), fb(n);
    copy(a.begin(), a.end(), fa.begin());
    copy(b.begin(), b.end(), fb.begin());
    ntt(fa, false);
    ntt(fb, false);
    for (int i = 0; i < n; i++) fa[i] = (long long)fa[i] * fb[i] % MOD;
    ntt(fa, true);
    fa.resize(sz);
    return fa;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    long long n;
    int p, k;
    cin >> n >> p >> k;

    // precompute p^i and p^{C(i,2)}
    vector<int> p_pow(k + 1), p_c2(k + 1);
    p_pow[0] = 1;
    for (int i = 1; i <= k; i++) p_pow[i] = (long long)p_pow[i - 1] * p % MOD;
    p_c2[0] = 1;
    for (int i = 1; i <= k; i++) p_c2[i] = (long long)p_c2[i - 1] * p_pow[i - 1] % MOD;

    // precompute inverses of p_c2
    vector<int> p_c2_inv(k + 1);
    for (int i = 0; i <= k; i++) p_c2_inv[i] = mod_pow(p_c2[i], MOD - 2);

    // p^n mod MOD
    int pn = mod_pow(p, n % (MOD - 1));

    // dp state: g_m for current m
    vector<int> g(k + 1, 0);
    g[0] = 1;
    long long cur_m_mod = 0;  // m modulo (MOD-1)

    // function to add one row: g_{m+1} from g_m
    auto add_one = [&](const vector<int>& g) {
        vector<int> ng(k + 1, 0);
        ng[0] = g[0];
        for (int i = 1; i <= k; i++) {
            ng[i] = ((long long)g[i] * p_pow[i] + g[i - 1]) % MOD;
        }
        return ng;
    };

    // function to double rows: g_{2m} from g_m, given m_mod = m mod (MOD-1)
    auto square = [&](const vector<int>& g, long long m_mod) {
        // base = p^{m_mod}
        int base = mod_pow(p, m_mod);
        vector<int> base_pow(k + 1);
        base_pow[0] = 1;
        for (int i = 1; i <= k; i++) base_pow[i] = (long long)base_pow[i - 1] * base % MOD;

        vector<int> A(k + 1), B(k + 1);
        for (int i = 0; i <= k; i++) {
            A[i] = (long long)g[i] * base_pow[i] % MOD * p_c2[i] % MOD;
            B[i] = (long long)g[i] * p_c2[i] % MOD;
        }
        vector<int> C = convolution(A, B);
        vector<int> res(k + 1, 0);
        for (int i = 0; i <= k; i++) {
            if (i >= (int)C.size()) break;
            int factor = p_c2_inv[i];
            res[i] = (long long)C[i] * factor % MOD;
        }
        return res;
    };

    // binary exponentiation: process bits from high to low
    for (int bit = 60; bit >= 0; bit--) {
        // square
        g = square(g, cur_m_mod);
        cur_m_mod = (cur_m_mod * 2) % (MOD - 1);
        // if current bit is 1, add one
        if ((n >> bit) & 1) {
            g = add_one(g);
            cur_m_mod = (cur_m_mod + 1) % (MOD - 1);
        }
    }

    // precompute prefix products of (p^n - p^i)
    vector<int> prod(k + 1, 1);
    for (int i = 1; i <= k; i++) {
        int diff = (pn - p_pow[i - 1] + MOD) % MOD;
        prod[i] = (long long)prod[i - 1] * diff % MOD;
    }

    // output answers
    for (int r = 0; r <= k; r++) {
        long long ans = (long long)g[r] * prod[r] % MOD;
        cout << ans << " \n"[r == k];
    }

    return 0;
}