【思路要点】
- 考虑容斥,一个显然的做法是暴力枚举 (k2)\binom{k}{2}(2k) 对点是否相等,将所有相等的点并在一起,可以用简单 O(N∗k)O(N*k)O(N∗k) 的动态规划计算方案数,乘上容斥系数加入答案即可。
- 可以发现,上述做法的动态规划部分的答案只与各个联通块的大小形成的可重集有关,有 O(f(k))O(f(k))O(f(k)) 种可能,其中 f(k)f(k)f(k) 表示 kkk 的整数拆分数。
- 考虑枚举这个可重集,计算可能的连边方案及其容斥系数。
- 假设该可重集为 a1,a2,a3,...,am (a1≤a2≤a3≤...≤am,a1+a2+a3+...+am=k)a_1,a_2,a_3,...,a_m\ (a_1≤a_2≤a_3≤...≤a_m,a_1+a_2+a_3+...+a_m=k)a1,a2,a3,...,am (a1≤a2≤a3≤...≤am,a1+a2+a3+...+am=k) ,可行的点的划分方案应当为 k!∏ci!∏ai!\frac{k!}{\prod c_i!\prod a_i!}∏ci!∏ai!k! ,其中 cjc_jcj 表示等于 jjj 的 aia_iai 个数。
- 并且,我们需要使得各个联通块联通,令一个连边方案的权值为 −1-1−1 的边数次方,记 fif_ifi 表示 iii 个点的完全图,使得图联通的连边方案的权值和, gig_igi 表示所有连边方案的权值和,则有 gi=(1−1)(i2)=[i=1],fi=1,f2=−1,fi=gi−∑j=1i−1(i−1j−1)fj∗gi−j (i≥3)g_i=(1-1)^{\binom{i}{2}}=[i=1],f_i=1,f_2=-1,f_i=g_i-\sum_{j=1}^{i-1}\binom{i-1}{j-1}f_j*g_{i-j}\ (i≥3)gi=(1−1)(2i)=[i=1],fi=1,f2=−1,fi=gi−∑j=1i−1(j−1i−1)fj∗gi−j (i≥3) ,即 fi=(−1)i+1(i−1)!f_i=(-1)^{i+1}(i-1)!fi=(−1)i+1(i−1)! ,那么可以得到可能的连边方案及其容斥系数的和为 k!∏fai∏ci!∏ai!\frac{k!\prod f_{a_i}}{\prod c_i!\prod a_i!}∏ci!∏ai!k!∏fai。
- 最终求得的方案数是无序的,需要将答案除去 k!k!k! 。
- 时间复杂度 O(N∗k∗f(k))O(N*k*f(k))O(N∗k∗f(k)) ,其中 f(40)=37338f(40)=37338f(40)=37338 ,大约可以在 30s30s30s 内得出 N=104,k=40N=10^4,k=40N=104,k=40 的答案。
【代码】
#include<bits/stdc++.h> using namespace std; const int MAXN = 1e4 + 5; const int MAXK = 55; const int P = 998244353; typedef long long ll; typedef long double ld; typedef unsigned long long ull; template <typename T> void chkmax(T &x, T y) {x = max(x, y); } template <typename T> void chkmin(T &x, T y) {x = min(x, y); } template <typename T> void read(T &x) { x = 0; int f = 1; char c = getchar(); for (; !isdigit(c); c = getchar()) if (c == '-') f = -f; for (; isdigit(c); c = getchar()) x = x * 10 + c - '0'; x *= f; } template <typename T> void write(T x) { if (x < 0) x = -x, putchar('-'); if (x > 9) write(x / 10); putchar(x % 10 + '0'); } template <typename T> void writeln(T x) { write(x); puts(""); } int n, k, ans, Max, cnt[MAXN], a[MAXK], dp[MAXK][MAXN]; int coef[MAXK], fac[MAXK], inv[MAXK], binom[MAXK][MAXK]; void update(int &x, int y) { x += y; if (x >= P) x -= P; } int getdp(int depth) { dp[0][0] = 1; for (int i = 1; i <= depth; i++) { memset(dp[i], 0, sizeof(dp[i])); static int tmp[MAXK]; memset(tmp, 0, sizeof(tmp)); for (int j = 0; j <= Max; j++) { update(tmp[j % a[i]], dp[i - 1][j]); dp[i][j] = tmp[j % a[i]]; } } int ans = 1; for (int i = 1; i <= n; i++) ans = 1ll * ans * dp[depth][cnt[i]] % P; return ans; } int getcoef(int depth) { int ans = 1; for (int i = 1; i <= depth; i++) ans = 1ll * ans * inv[a[i]] % P * coef[a[i]] % P; int last = 0; a[depth + 1] = 0; for (int i = 1; i <= depth; i++) if (a[i] != a[i + 1]) { ans = 1ll * ans * inv[i - last] % P; last = i; } return ans; } void work(int lft, int depth, int last) { if (lft == 0) { update(ans, 1ll * getcoef(depth - 1) * getdp(depth - 1) % P); return; } if (last >= 2) work(lft, depth, last - 1); if (lft >= last) { a[depth] = last; work(lft - last, depth + 1, last); } } int power(int x, int y) { if (y == 0) return 1; int tmp = power(x, y / 2); if (y % 2 == 0) return 1ll * tmp * tmp % P; else return 1ll * tmp * tmp % P * x % P; } void calcoef() { int n = 40; for (int i = 0; i <= n; i++) { binom[i][0] = 1; for (int j = 1; j <= i; j++) binom[i][j] = (binom[i - 1][j - 1] + binom[i - 1][j]) % P; } fac[0] = inv[0] = 1; for (int i = 1; i <= n; i++) { fac[i] = 1ll * fac[i - 1] * i % P; inv[i] = power(fac[i], P - 2); } coef[1] = 1, coef[2] = P - 1; for (int i = 3; i <= n; i++) { int tans = 0; update(tans, P - 1ll * binom[i - 1][i - 2] * coef[i - 1] % P); coef[i] = tans; } } void calccnt() { for (int i = 1; i <= n; i++) { int tmp = i; for (int j = 2; j * j <= tmp; j++) while (tmp % j == 0) { tmp /= j; cnt[j]++; } if (tmp != 1) cnt[tmp]++; } for (int i = 1; i <= n; i++) chkmax(Max, cnt[i]); } int main() { freopen("final.in", "r", stdin); freopen("final.out", "w", stdout); read(n), read(k); calccnt(); calcoef(); work(k, 1, k); writeln(ans); return 0; }