看到次幂就得老老实实的把答案算出来。
不难发现对于一个排列的答案，它的贡献是 \(i\to p_i\) 后所有环大小的 \(\text{lcm}\)，于是设 \(f(i,j)\) 表示 \(i\) 个点组成的若干个环，环大小的 \(\text{lcm}\) 是 \(j\) 的方案数。
每次加入一个新的大小为 \(k\) 的环的时候就乘上选出 \(k\) 个点的方案数以及它们构成环的方案数。显然它们构成环的方案数是 \((k-1)!\)
但是这样会算重，因为一个环加入的顺序没有关系，那么强制选剩下的点中编号最小的点就行了。
于是转移为
\[f(i+k,\text{lcm}(j,k))\gets f(i+k,\text{lcm}(j,k))+f(i,j)\times \binom{n-i-1}{k-1}\times (k-1)! \]第二维特别大，于是用 map 保存，第二维的非空状态在 \(n=50\) 时为 \(1056\)。
Code：
#include <bits/stdc++.h>
using namespace std;
#define fi first
#define se second
typedef long long ll;
const int N = 55, mod = 998244353;
int n, k, ans;
int fac[N], inv[N];
map <ll, int> f[N];
int qpow(int x, int y) {
int res = 1;
while (y) {
if (y & 1) res = 1ll * res * x % mod;
x = 1ll * x * x % mod;
y >>= 1;
}
return res;
}
void init(int maxn) {
fac[0] = 1;
for (int i = 1; i <= maxn; ++i) fac[i] = 1ll * fac[i - 1] * i % mod;
inv[maxn] = qpow(fac[maxn], mod - 2);
for (int i = maxn - 1; ~i; --i) inv[i] = 1ll * inv[i + 1] * (i + 1) % mod;
}
int C(int n, int m) {
if (n < 0 || m < 0 || n < m) return 0;
return 1ll * fac[n] * inv[m] % mod * inv[n - m] % mod;
}
void add(int &a, int b) {
a += b;
if (a >= mod) a -= mod;
}
ll lcm(ll x, ll y) { return x / __gcd(x, y) * y; }
int main() {
scanf("%d%d", &n, &k);
init(n);
f[0][1] = 1;
for (int i = 0; i < n; ++i)
for (auto x : f[i])
for (int j = 1; i + j <= n; ++j)
add(f[i + j][lcm(x.fi, j)], 1ll * C(n - i - 1, j - 1) * x.se % mod * fac[j - 1] % mod);
int ans = 0;
for (auto x : f[n]) add(ans, 1ll * x.se * qpow(x.fi, k) % mod);
printf("%d", ans);
return 0;
}