题解 UOJ #99. 【集训队互测2015】普罗达科特
想了一个上午，看题解发现方向错了/kk
整个问题有两个子问题，我们分三部分叙述做法，分别为将两个子问题递归到一个「基本问题」的方法以及「基本问题」的做法。
其中「基本问题」为给出一个 \(k\) 的整数拆分 \(\lambda \vdash k\)，并计算将 \(n\) 个带标号元素划分到带标号集合中，恰有 \(\lambda_i\) 个大小为 \(i\) 的集合，要求集合内部元素相同，集合之间元素无限制的方案数。
第一部分 集合个数
先假装任意元素带标号，最后只需要除以 \(k!\)。考虑容斥。条件变成所有元素都互不相同，朴素容斥即枚举 \(2^{\binom k2}\) 种可能的关系分别计算。
考虑优化。发现两种连通块大小集合相同的关系方案数相同，于是剩余的问题无非是类似「基本问题」的形式，但是我们需要对每个整数拆分算一下方案数。
考虑一个整数拆分 \(\lambda\vdash k\)。首先是将 \(k\) 个元素分到所有集合中的方案数，即 \(\frac{k!}{\prod_i (i!)^{\lambda_i}}\)。而「基本问题」的集合是有标号的，实际的集合要求是无标号的，所以还需要除掉一个 \(\prod_i \lambda_i!\)。剩下的是考虑这个整数拆分在朴素容斥中的系数和。
发现集合之间独立，只需要分别计算然后乘起来即可。于是相当于计算所有大小为 \(n\) 的带标号连通图 \((V,E)\) 的 \((-1)^{|E|}\) 的和。容易联想到 \(\mathcal O(n^4)\) 的带边数的连通图计数的做法。考虑设 \(f_{i,j}\) 表示 \(i\) 个连通块，\(j\) 条边的方案数，\(g_{i,j}\) 为钦定出 \(i\) 个内部无要求，之间不相连的集合，\(j\) 条边的方案数，有等式：
反演即有
要求的即
考虑到后面那个和式是每个无要求的集合的和式的乘积，则如果有大小 \(>1\) 的连通块则贡献一定是 \(0\)，也就是说有贡献的只有 \(\lang1,1,\cdots,1\rang\) 这一种拆分，也就是说上面那个柿子其实就是 \((n-1)!(-1)^{n-1}\)。
于是我们的答案为
其中 \(f(\lambda)\) 为 \(\lambda\) 对应的「基本问题」的答案。
第二部分 可重集个数
容斥原理在可重集上的推广不是非常方便。我们考虑另一种做法。发现我们要求的即所有「在进行一次 \(S_k\) 中的置换后等价的方案」的个数。这启发我们考虑 Burnside 引理。
考虑枚举轨道的划分，于是一个大小为 \(n\) 的轨道对应恰 \((n-1)!\) 种排列，而把元素划分到轨道中的方案是和上一个部分类似的。于是不难得到：
这个柿子和上面的柿子只相差一个 \(-1\) 的系数，这当中可能有更加本质的关系。但仅就这道题而言，这个式子已经足够了。
第三部分 「基本问题」的答案
首先考虑容斥，假设我们钦定了整数拆分中的大小为 \(i\) 的集合有 \(c_i\) 个超过了限制，则对应的容斥系数为
然后质因子独立，分别计算然后乘起来即可。这部分是
只需要对每一对 \((a,b)\) 计算上述柿子然后乘起来即可。
朴素做法是暴力线性递推，但是复杂度爆炸了。我们来考虑先把上面的部分乘开变成 \(F(x^{b_i+1})\)，然后后面的形式的部分是形如
的形式。根据 dls 给的结论，记 \(L=\mathrm{lcm}(t_1,t_2,\cdots, t_m)\)，这个柿子的 \(kL+r(0\leq r<L)\) 处的系数是关于 \(k\) 的 \(m\) 次多项式。于是只需要暴力求出前 \(L(k+1)\) 项的系数即可单次 \(\mathcal O(k)\) 查询这个问题。
以下是对上述问题的复杂度分析。
首先是枚举 \(\lambda\) 和 \(c\) 的复杂度，这部分一共有 \(\sum_{\lambda\vdash k} \prod_i(\lambda_i+1)\) 种状态。写个暴搜跑一下发现只有 \(Q_{25}=129512\) 种状态。
然后是计算那个查询的复杂度，对于 \(Q_k\) 种状态需要把分子乘开，这部分复杂度是 \(\mathcal O(Q_kk^2)\)。而后面的查询注意到只和 \(\lambda\) 有关系，于是可以对每个 \(\lambda\) 暴力计算，然后对于 \(Q_k\) 种状态 \(\mathcal O(nk)\) 查询。而对于 \(\lambda\) 的那部分的复杂度是 \(k^2\sum_{\lambda\vdash k}\mathrm{lcm}(\lambda)=k^2R_k\)，写个暴搜发现这部分只有 \(R_{25}=123292\) 次计算。
当然还有插值的复杂度，这部分是 \(nk^2P_k\)，其中 \(P_k\) 为划分数。
于是总复杂度为 \(\mathcal O((k^2+nk)Q_k+k^2R_k+nk^2P_k)\)，其中 \(Q_{25}=129512,R_{25}=123292,P_{25}=1958\)，算一下发现能过。
#include <bits/stdc++.h>
using ll = long long;
using poly = std::vector<ll>;
constexpr int maxn = 55;
constexpr int maxk = 30;
constexpr ll mod = 1E+9 + 21;
inline ll fsp(ll a, ll b, ll res = 1) {
for(a %= mod; b; a = a * a % mod, b >>= 1)
b & 1 && (res = res * a % mod); return res;
}
inline ll sgn(int x) { return x & 1 ? -1 : 1; }
inline poly mul(poly F, poly G) {
poly res; res.resize(F.size() + G.size() - 1);
for(int i = 0; i < F.size(); ++i)
for(int j = 0; j < G.size(); ++j)
(res[i + j] += F[i] * G[j]) %= mod;
return res;
}
inline poly inv(poly F, int deg) {
poly res;
res.push_back(fsp(F[0], mod - 2));
for(int i = 1; i <= deg; ++i) {
ll tmp = 0;
for(int j = 1; j < F.size() && j <= i; ++j)
(tmp += F[j] * res[i - j]) %= mod;
tmp = -res[0] * tmp % mod, res.push_back(tmp);
}
return res;
}
int n, k;
ll a[maxn], b[maxn], ans1, ans2;
ll Fac[maxk], Inv[maxk], val[maxn][maxk];
inline ll binom(int n, int k) {
if(k < 0 || n < k) return 0;
return Fac[n] * Inv[k] % mod * Inv[n - k] % mod;
}
int siz[maxk], cnt[maxk];
inline void getcoef() {
poly P{ 1 };
for(int i = 1; i <= k; ++i) {
poly tmp; tmp.resize(i + 1);
tmp[0] = 1, tmp[i] = -1;
for(int j = 1; j <= cnt[i]; ++j) P = mul(P, tmp);
}
ll res = 1;
for(int i = 1; i <= n; ++i) {
ll tmp = 0;
for(int j = 0; j < P.size(); ++j)
(tmp += P[j] * val[i][j]) %= mod;
res = res * tmp % mod;
}
ll coef1 = 1, coef2 = 1;
for(int i = 1; i <= k; ++i) {
for(int j = 1; j <= siz[i]; ++j) {
res = res * Inv[i] % mod;
coef1 = coef1 * Fac[i - 1] % mod * sgn(i - 1);
coef2 = coef2 * Fac[i - 1] % mod;
}
res = res * Inv[cnt[i]] % mod * Inv[siz[i] - cnt[i]] % mod * sgn(cnt[i]);
}
(ans1 += res * coef1) %= mod, (ans2 += res * coef2) %= mod;
}
inline void DFS2(int p) {
if(p == k + 1) return getcoef();
for(int i = 0; i <= siz[p]; ++i) cnt[p] = i, DFS2(p + 1);
}
ll pre[maxk], suf[maxk];
inline ll getval(poly val, ll x) {
x %= mod;
if(x < val.size()) return val[x];
int n = val.size() - 1;
pre[0] = 1, suf[n] = 1;
for(int i = 1; i <= n; ++i) {
pre[i] = pre[i - 1] * (x - i + 1) % mod;
suf[n - i] = suf[n - i + 1] * (x - n + i - 1) % mod;
}
ll res = 0;
for(int i = 0; i <= n; ++i)
(res += Inv[i] * Inv[n - i] % mod * pre[i] % mod * suf[i] % mod * sgn(n - i) * val[i]) %= mod;
return res;
}
inline void calc() {
poly G{ 1 }; int lcm = 1;
for(int i = 1; i <= k; ++i) {
poly tmp; tmp.resize(i + 1);
tmp[0] = 1, tmp[i] = -1;
for(int j = 1; j <= siz[i]; ++j) G = mul(G, tmp);
if(siz[i]) lcm = lcm * i / std::__gcd(lcm, i);
}
G = inv(G, lcm * (k + 1) - 1);
for(int i = 1; i <= n; ++i)
for(int j = 0; j <= k; ++j) {
if(a[i] / (b[i] + 1) < j) continue;
ll x = a[i] - j * (b[i] + 1), r = x % lcm;
poly vec;
for(int h = 0; h <= k; ++h)
vec.push_back(G[h * lcm + r]);
val[i][j] = getval(vec, (x - r) / lcm);
}
DFS2(1);
}
inline void DFS(int p, int rem) {
if(p == k + 1) {
if(!rem) calc();
return;
}
for(int i = 0; i * p <= rem; ++i)
siz[p] = i, DFS(p + 1, rem - i * p);
}
int main() {
scanf("%d%d", &n, &k);
Fac[0] = 1;
for(int i = 1; i <= k; ++i) Fac[i] = Fac[i - 1] * i % mod;
Inv[k] = fsp(Fac[k], mod - 2);
for(int i = k; i >= 1; --i) Inv[i - 1] = Inv[i] * i % mod;
for(int i = 1; i <= n; ++i) scanf("%lld", &a[i]);
for(int i = 1; i <= n; ++i) scanf("%lld", &b[i]);
DFS(1, k);
printf("%lld %lld\n", (ans1 + mod) % mod, (ans2 + mod) % mod);
}
