发现这个 \(m\) 很大，且这个式子是 \(\times\)。
一个想法是拆成 \(m = \prod {p_i}^{e_i}(p_i\in \mathbb{P})\) 然后对于 \(M = p_i^{e_i}\) 依次考虑 \(b_i = a_i\bmod M\) 和 \(N = n\bmod M\)。
根据 \(\text{CRT}\)，对于任意一个 \(M\) 得到的不同的 \(b_i\) 对于最后的 \(a_i\) 都是不同的，于是方案数可以直接相乘。
接下来考虑求解 \(M, N\)。
因为最后的 \(b_i\) 是 \(\bmod M\) 意义下的，可以把值域当作 \([1, M]\)。
记 \(f(c, i)\) 为 \(c\) 个 \([1, M]\) 的数乘积 \(p\) 的指数为 \(i\) 的方案数。
首先对于 \(N = 0\)，那么就是要求 \(p\) 的指数需要 \(\ge e\)，可以考虑容斥，总方案数即为 \(M^k\)，不合法的方案数为 \(\sum\limits_{i = 0}^{e - 1} f(k, i)\)。
然后对于 \(N \not = 0\)，考虑令 \(N = r\times p^z(p\nmid r)\)。
先得到 \(f(k - 1, i)\)，然后用 \(b_k\) 去凑出 \(r\times p^z\)。
那么对于 \(u\times p^x(p\nmid u)\)，就需要 \(b_k = v\times p^{z - x}(p\nmid v)\) 满足 \(uv\times p^z\equiv r\times p^z(\bmod M)\)。
一种直观的想法是令 \(v = u^{-1}\times r\)，因为 \(p\nmid u\) 肯定存在 \(u^{-1}\)。
进一步，发现 \(v = u^{-1}\times r + kp^{e - x}\) 也是合法的，这是因为后面那部分乘后的结果为 \(uk\times p^e\bmod M = 0\)。
那么就会每 \(p^{e - x}\) 个值就会存在一个合法的 \(v\)，就会有 \(p^x\) 个合法的 \(v\)。
于是这部分的方案数为 \(\sum\limits_{i = 0}^e f(k - 1, i)\times p^i\)。
接下来就只需要考虑求 \(f\) 了。
考虑到令 \(g(i)\) 为 \([1, M]\) 范围内的数 \(p\) 的指数为 \(i\) 的方案数，显然有 \(g(i) = \varphi(\frac{M}{p^i})(0\le i\le e)\)。
那么就有转移式 \(f(k, i) = \sum\limits_{j = 0}^i f(k - 1, j)g(i - j)\)。
进一步的，能发现这是卷积形式。
就可以写成 \(G(x) = \sum\limits_{i = 0}^e g(i)x^i\)，有 \(F(k, x) = \sum\limits_{i = 0}^e f(k, i)x^i = G^k(x)\)。
那么对多项式做快速幂即可，因为项数是 \(O(\log m)\) 的可以直接暴力乘。
时间复杂度 \(O(\sqrt{m} + \log^2 m\log k)\)。
代码
#include<bits/stdc++.h>
using i64 = long long;
const i64 mod = 998244353;
inline i64 qpow(i64 a, i64 b) {i64 v = 1; a %= mod; while (b) b & 1 && ((v *= a) %= mod), (a *= a) %= mod, b >>= 1; return v;}
struct node_t {
std::vector<i64> a; int n;
inline node_t(int n_ = 0) {n = n_, a.resize(n + 1);}
inline i64 operator [] (int x) const {return a[x];}
inline i64 &operator [] (int x) {return a[x];}
inline node_t operator * (const node_t &o) const {
node_t v(n);
for (int i = 0; i <= n; i++) for (int j = 0; i + j <= n; j++) (v[i + j] += a[i] * o[j]) %= mod;
return v;
}
inline node_t &operator *= (const node_t &o) {return (*this) = (*this) * o;}
};
node_t qpow(node_t a, i64 b) {
node_t v(a.n); v[0] = 1;
while (b) b & 1 && (v *= a, 1), a *= a, b >>= 1;
return v;
}
int main() {
i64 k, n, m; scanf("%lld%lld%lld", &k, &n, &m);
std::map<i64, int> mp;
for (i64 x = 2; x * x <= m; x++) if (m % x == 0) {
while (m % x == 0) m /= x, mp[x]++;
}
if (m > 1) mp[m] = 1;
i64 ans = 1;
for (auto _ : mp) {
i64 p = _.first; int e = _.second;
i64 M = 1;
for (int j = 1; j <= e; j++) M *= p;
node_t B(e); i64 M_ = M;
for (int i = 0; i < e; i++) B[i] = (M_ - M_ / p) % mod, M_ /= p;
B[e] = 1;
i64 N = n % M, v;
if (! N) {
node_t qp = qpow(B, k);
v = qpow(M, k);
for (int i = 0; i < e; i++) (v += mod - qp[i]) %= mod;
} else {
node_t qp = qpow(B, k - 1);
int c = 0;
while (N % p == 0) N /= p, c++;
v = 0;
for (int i = 0; i <= c; i++) (v += qp[i] * qpow(p, i)) %= mod;
}
(ans *= v) %= mod;
}
printf("%lld\n", ans);
return 0;
}