题目理解

我们要求的是：

• n 位十进制数（允许前导 0）
• 是 p 的倍数
• 各位数字之和 ≤ m_i，其中 m_i = 0, 1, 2, ..., m
• 对每个 m_i 输出答案，模 998244353

---

思路分析

1. 状态定义

对于这种“n 位十进制数”且 n 很大（最大到 (10^9)）的情况，直接枚举每一位是不可能的。

我们考虑 DP：

设
[ dp[len][r][s] ] 表示长度为 len，模 p 的余数为 r，数位和为 s 的数字个数。

最终对于给定的 m_i，答案是： [ \sum_{r=0}^{p-1} \sum_{s=0}^{m_i} dp[n][r][s] \quad \text{但只取 r=0 的部分} ] 因为要求是 p 的倍数，所以余数必须为 0。

---

2. 转移方程

从 len 到 len+1，我们枚举新的一位数字 d（0 到 9）：

• 新的余数：( r' = (r \times 10 + d) \bmod p )
• 新的数位和：( s' = s + d )

所以： [ dp[len+1][r'][s'] += dp[len][r][s] ]

初始状态：dp[0][0][0] = 1（0 位数，余数 0，和 0）

---

3. 问题规模

• n 最大 (10^9)，不能直接循环长度。
• p 最大 50，m 最大 1000（但 m 可能很大，不过我们只关心 ≤ m 的和）。
• 直接二维状态 (r, s) 大小是 p * (m+1)，最大 50 * 1001 ≈ 50050。

我们需要用 矩阵快速幂 来加速递推。

---

4. 矩阵构造

我们把状态 (r, s) 展平成一维索引：
idx = r * (m+1) + s，总状态数 P = p * (m+1)。

那么递推关系可以写成： [ \mathbf{v}{len+1} = M \cdot \mathbf{v}{len} ] 其中 (M) 是 P × P 的转移矩阵，M[idx2][idx1] 表示从状态 (r,s) 通过一个数字 d 转移到 (r',s') 的转移系数（这里为 1 如果 s' ≤ m，否则 0）。

---

5. 最终答案

初始向量 v0 只有 (r=0, s=0) 位置为 1，其余为 0。

计算： [ v_n = M^n \cdot v_0 ] 然后对于每个 m_i，答案是： [ \text{ans}[m_i] = \sum_{s=0}^{m_i} v_n[(r=0, s)] ] 即余数为 0，数位和 ≤ m_i 的所有方案数。

---

6. 优化

• 矩阵大小 P = p*(m+1) 最大可能 50*1001 = 50050，但 50000 维的矩阵快速幂不可行（O(P^3 log n) 太大）。
• 必须注意：m 最大 1000 时，p 最大 16（见数据范围），所以 P ≤ 16*1001 = 16016，仍然太大，无法直接矩阵快速幂。

因此需要优化状态：
我们其实不需要同时存储 (r,s) 的所有组合，因为转移时 s 只增加 0~9，所以我们可以用 滚动数组 按 s 从小到大处理，但 n 很大时仍然需要矩阵快速幂。

---

7. 更优方法：按余数矩阵，数位和作为生成函数

我们可以把数位和视为生成函数的指数：

设
[ F_len(x) = [f_{len,0}(x), f_{len,1}(x), \dots, f_{len,p-1}(x)] ] 其中 (f_{len,r}(x)) 是一个多项式，(x^s) 的系数表示长度 len、余数 r、数位和 s 的方案数。

那么转移是： [ f_{len+1, r'}(x) += x^d \cdot f_{len, r}(x) ] 其中 (r' = (r*10 + d) \bmod p)。

这等价于： [ F_{len+1} = F_{len} \cdot A(x) ] 其中 (A(x)) 是一个 p × p 的矩阵，每个元素是一个多项式（次数 ≤ 9），表示从余数 i 到余数 j 的转移多项式。

---

8. 多项式乘法与快速幂

这样矩阵大小只有 p × p，每个元素是次数 ≤ 9 的多项式。
多项式乘法可以用 FFT/NTT 加速（模 998244353 正好是 NTT 模数）。

那么计算 (A(x)^n) 再乘以初始向量 [1, 0, 0, ...]（多项式形式），就得到 (F_n(x))。

最终答案： [ \text{ans}[m_i] = \text{多项式 } f_{n,0}(x) \text{ 的 } x^0 \text{ 到 } x^{m_i} \text{ 的系数和} ] 即前缀和。

---

9. 复杂度

• 矩阵大小 p ≤ 50，多项式次数 ≤ 9，乘法用 NTT，每次矩阵乘法 O(p^3 * m log m) 中的 m 是多项式次数，其实固定为 9，所以是 O(p^3 * K log K) 其中 K=16（向上取2的幂）。
• 快速幂 O(log n) 次乘法，完全可行。

---

10. 实现细节

• 用 NTT 做多项式乘法，模 998244353。
• 矩阵快速幂时，矩阵元素是多项式。
• 初始向量是 (1, 0, ..., 0) 对应的多项式是 1（常数多项式）。
• 最后提取 f_{0}(x) 的前缀和。

---

代码框架（C++）

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

const int MOD = 998244353;

// 多项式乘法 (NTT)
vector<int> multiply(const vector<int>& a, const vector<int>& b) {
    // 这里实现 NTT 卷积，长度 = len(a)+len(b)-1
    // 简化为暴力，因为次数最多9+9=18
    int n = a.size(), m = b.size();
    vector<int> res(n + m - 1, 0);
    for (int i = 0; i < n; i++)
        for (int j = 0; j < m; j++)
            res[i+j] = (res[i+j] + 1LL * a[i] * b[j]) % MOD;
    return res;
}

// 多项式加法
vector<int> add(const vector<int>& a, const vector<int>& b) {
    int n = max(a.size(), b.size());
    vector<int> res(n, 0);
    for (int i = 0; i < a.size(); i++) res[i] = (res[i] + a[i]) % MOD;
    for (int i = 0; i < b.size(); i++) res[i] = (res[i] + b[i]) % MOD;
    return res;
}

struct PolyMatrix {
    int p;
    int deg;
    vector<vector<vector<int>>> M; // M[i][j] 是一个多项式

    PolyMatrix(int p, int deg) : p(p), deg(deg) {
        M.assign(p, vector<vector<int>>(p, vector<int>(deg+1, 0)));
    }

    PolyMatrix operator*(const PolyMatrix& other) const {
        PolyMatrix res(p, deg);
        for (int i = 0; i < p; i++)
            for (int k = 0; k < p; k++)
                for (int j = 0; j < p; j++) {
                    auto prod = multiply(M[i][k], other.M[k][j]);
                    // 截断到 deg
                    if (prod.size() > deg+1) prod.resize(deg+1);
                    res.M[i][j] = add(res.M[i][j], prod);
                }
        return res;
    }
};

PolyMatrix mat_pow(PolyMatrix A, long long n) {
    int p = A.p;
    int deg = A.deg;
    PolyMatrix res(p, deg);
    // 单位矩阵
    for (int i = 0; i < p; i++) res.M[i][i][0] = 1;
    while (n) {
        if (n & 1) res = res * A;
        A = A * A;
        n >>= 1;
    }
    return res;
}

int main() {
    long long n;
    int p, m;
    cin >> n >> p >> m;

    // 构造转移矩阵 A
    PolyMatrix A(p, m); // 注意这里度数是 m，不是 9，因为最终需要 ≤ m 的和
    // 但这样多项式乘法代价大，需要优化：实际上我们只需要度数 ≤ m，且初始度数=9
    // 优化：用真正的小度数乘法，最后再截断

    for (int r = 0; r < p; r++) {
        for (int d = 0; d <= 9; d++) {
            int r2 = (r * 10 + d) % p;
            if (d <= m) A.M[r][r2][d] = (A.M[r][r2][d] + 1) % MOD;
        }
    }

    PolyMatrix T = mat_pow(A, n);

    // 初始向量 [1,0,0,...]
    vector<int> init(p, 0);
    init[0] = 1;

    // 结果向量 = init * T
    vector<int> res_poly(m+1, 0);
    for (int j = 0; j < p; j++) {
        for (int k = 0; k <= m; k++) {
            res_poly[k] = (res_poly[k] + 1LL * init[j] * T.M[j][0][k]) % MOD;
        }
    }

    // 前缀和
    vector<int> ans(m+1);
    ans[0] = res_poly[0];
    for (int i = 1; i <= m; i++) {
        ans[i] = (ans[i-1] + res_poly[i]) % MOD;
    }

    for (int i = 0; i <= m; i++) {
        cout << ans[i] << " \n"[i == m];
    }

    return 0;
}

---

总结

这道题的核心是：

1. 将“数位和”和“模 p 余数”结合成二维状态。
2. 用多项式表示数位和的分布，将 DP 转移转化为矩阵形式。
3. 利用矩阵快速幂 + 多项式乘法（NTT）处理大 n。
4. 最终取余数 0 对应的多项式的前缀和作为答案。

这样就可以在 (O(p^3 \cdot K \log K \cdot \log n)) 的复杂度内解决，其中 (K) 是多项式最大次数（这里 ≤ m，但可通过优化限制到 9，再最终扩展到 m）。