[洛谷]P5401 [CTS2019] 珍珠 题解
题意概述
有 \(D\) 种珍珠，每种有无限颗，现在等概率的从 \(D\) 种珍珠中抽 \(n\) 次珍珠，每次抽 \(1\) 个珍珠，记第 \(i\) 种珍珠最后一共抽到 \(cnt_i\) 个，求有要多少抽法使得 \(\sum_{i=1}^{n}\lfloor \frac{cnt_i}{2}\rfloor \ge m\) 的方案数。两种方案不同当且仅当存在一轮抽取，一种方案抽到第 \(i\) 种珍珠，另一种方案抽到第 \(j\) 种珍珠， \(i\not=j\) 。
\(n\le 1e9,m\le 1e9,D\le1e5\)
题目分析
容易将题意转化，转化为求 \(res=\sum_{i=1}^{n}[cnt_i\equiv 1\mod 2] \le n-2\times m\) ，记 \(lim=\min(D,n-2m)\) ，那么容易构造 \(EGF\) ：
\[\begin{aligned} ans&=\sum_{k=0}^{mx}\binom{D}{k}\times n!\times [x^n](\sum_{i\ge0}x^{2i+1})^k\times(\sum_{i\ge 0}x^{2i})^{D-k}\\ &=n!\times [x^n]\sum_{k=0}^{mx}\binom{D}{k}\times (\frac{e^x-e^{-x}}{2})^k\times(\frac{e^x+e^{-x}}{2})^{D-k}\\ &=\frac{n!}{2^D}\times[x^n]e^{-Dx}\times\sum_{k=0}^{mx}\binom{D}{k}\times (e^{2x}-1)^k\times (e^{2x}+1)^{D-k}\\ &=\frac{n!}{2^D}\times[x^n]e^{-Dx}\times\sum_{k=0}^{mx}\binom{D}{k}\times \sum_{i=0}^{k}(-1)^{k-i}\binom{k}{i}e^{2ix}\times \sum_{j=0}^{D-k}\binom{D-k}{j}e^{2jx}\\ \end{aligned} \]这样做的话也是可以往下做的，但是因为我数学太弱了，并没有做出来（但下文会继续写到如何按照这种方法做）。
可以发现这样做不好做的困难主要在于对于 \(F(x)\) 每一项的系数出现了两个不固定的 \(EGF\) 相乘，所以说导致无法出现一个单一的和式（即将 \((e^{2x}-1)^k(e^{2x}+1)^{D-k}\) 二项式展开之后会出现 \(\sum_{i=0}^{k}\cdots \times \sum_{j=0}^{D-k}\cdots\)），无法直接取出 \(EGF\) 的 \(x^n\) 项系数，所以其实可以考虑改变一下求解的方式。
我们可以尝试构造一个多项式 \(F(x)\) ，使得 \(F(x)\) 的 \(x^i\) 的系数为 \(res= i\) 的方案数。
然后设 \(f(x)\) 的 \(x^i\) 项系数表示 \(res\) 至少为 \(i\) 的方案数，那么这样就有：
\[\begin{aligned} &F_i(x)=\sum_{j=i}^{D}(-1)^{j-i}\binom{j}{i}f_j(x)\\ &f_i(x)=\binom{D}{i}\times n!\times [x^n] (\frac{e^x-e^{-x}}{2})^i\times (e^x)^{D-i} \end{aligned} \]这样求解 \(F(x)\) 只需要求出 \(f(x)\) 然后卷积即可，而求 \(f(x)\) 的难度也较低：
\[\begin{aligned} f_i(x)&=\frac{\binom{D}{i}\times n!}{2^i}[x^n]e^{-ix}(e^{2x}-1)^i\times e^{(D-i)x}\\ &=\frac{\binom{D}{i}\times n!}{2^i}[x^n]e^{(D-2i)x}\times (e^{2x}-1)^i\\ &=\frac{\binom{D}{i}\times n!}{2^i}[x^n]\sum_{j=0}^{i}(-1)^{i-j}\binom{i}{j}e^{(D-2i)x}\times e^{2jx}\\ &=\frac{D!\times n!}{2^i\times i!\times (D-i)!}[x^n]\sum_{j=0}^{i}(-1)^{i-j}\frac{i!}{j!(i-j)!}e^{(D-2(i-j))x}\\ &=\frac{D!}{2^i\times (D-i)!}\sum_{j=0}^{i}(-1)^{i-j}\times\frac{1}{j!(i-j)!}\times(D-2(i-j))^n\\ &=\frac{D!}{2^i\times (D-i)!}\sum_{j=0}^{i}\frac{1}{j!}\times \frac{(-1)^{i-j}\times (D-2(i-j))^n}{(i-j)!} \end{aligned} \]对于求 \(F(x)\) 的卷积就比较套路，令 \(G_i(x)=F_{D-i}(x),\ g_j(x)=g_{D-j}(x)\) ，则：
\[\begin{aligned} &G_{D-i}(x)&=&\sum_{j=i}^{D}(-1)^{j-i}\frac{j!}{i!(j-i)!}g_{D-j}(x)\\ &G_{i}(x)&=&\sum_{j=0}^{i}(-1)^{i-j}\frac{(D-j)!}{(D-i)!(i-j)!}g_{j}(x)\\ & &=&\frac{1}{(D-i)!}\sum_{j=0}^{i}[(D-j)!\times g_j(x)]\times \frac{(-1)^{i-j}}{(i-j)!} \end{aligned} \]综上，两次卷积就可以求出答案，\(ans=\sum_{i=0}^{\min(D,n-2m)} F_i(x)\) ，时间复杂度 \(O(D\log_2n)\) ，下面附上代码：
#include <bits/stdc++.h>
using namespace std;
const int N = 3e5 + 30;
const int M = 3e5;
const int mod = 998244353;
const int pr = 3;
const int ig = 332748118;
inline int add(int x, int y) {
x += y;
return x >= mod ? x - mod : x;
}
inline int del(int x, int y) {
x -= y;
return x < 0 ? x + mod : x;
}
inline int read() {
char ch = getchar(); int x = 0;
while(!isdigit(ch)) ch = getchar();
while(isdigit(ch)) {x = (x << 3) + (x << 1) + ch - 48; ch = getchar();}
return x;
}
int qpow(int a, int b) {
int res = 1;
while(b) {
if(b & 1) res = 1ll * res * a % mod;
a = 1ll * a * a % mod;
b >>= 1;
}
return res;
}
struct node {
int n, mx = 1, f[N], g[N];
void NTT(int *a, int len) {
int x, y, g, pw;
for(int j = len >> 1; j >= 1; j >>= 1) {
g = qpow(pr, (mod - 1) / (j << 1));
for(int i = 0; i < len; i += (j << 1)) {
pw = 1;
for(int k = 0; k < j; ++k, pw = 1ll * pw * g % mod) {
x = a[i + k]; y = a[i + j + k];
a[i + k] = add(x, y);
a[i + j + k] = 1ll * pw * del(x, y) % mod;
}
}
}
}
void INTT(int *a, int len) {
int x, y, g, pw;
for(int j = 1; j < len; j <<= 1) {
g = qpow(ig, (mod - 1) / (j << 1));
for(int i = 0; i < len; i += (j << 1)) {
pw = 1;
for(int k = 0; k < j; ++k, pw = 1ll * pw * g % mod) {
x = a[i + k]; y = 1ll * pw * a[i + j + k] % mod;
a[i + k] = add(x, y);
a[i + j + k] = del(x, y);
}
}
}
int Inv = qpow(len, mod - 2);
for(int i = 0; i < len; ++i) a[i] = 1ll * a[i] * Inv % mod;
}
void mul(int *a, int *b, int len) {
NTT(a, len); NTT(b, len);
for(int i = 0; i < len; ++i) a[i] = 1ll * a[i] * b[i] % mod;
INTT(a, len);
}
void init() {while(mx <= n + n) mx <<= 1;}
void clear() {for(int i = n + 1; i < mx; ++i) f[i] = g[i] = 0;}
}F;
int D, n, m, k, fac[N], inv[N], ip[N], ans;
int main() {
D = read(); n = read(); m = read(); k = min(n - 2 * m, D);
if(k < 0) return puts("0"), 0;
fac[0] = ip[0] = 1;
for(int i = 1; i <= D; ++i) ip[i] = 499122177ll * ip[i - 1] % mod;
for(int i = 1; i <= D; ++i) fac[i] = 1ll * fac[i - 1] * i % mod;
inv[D] = qpow(fac[D], mod - 2);
for(int i = D; i >= 1; --i) inv[i - 1] = 1ll * inv[i] * i % mod;
F.n = D; F.init();
for(int i = 0; i <= F.n; ++i) F.g[i] = inv[i];
for(int i = 0; i <= F.n; ++i) F.f[i] = 1ll * ((i & 1) ? mod - 1 : 1) * qpow(del(D, 2 * i), n) % mod * inv[i] % mod;
F.mul(F.f, F.g, F.mx); F.clear();
for(int i = 0; i <= F.n; ++i) F.f[i] = 1ll * fac[D] * ip[i] % mod * F.f[i] % mod * fac[i] % mod * inv[D - i] % mod;
reverse(F.f, F.f + F.n + 1);
for(int i = 0; i <= F.n; ++i) F.g[i] = 1ll * ((i & 1) ? mod - 1 : 1) * inv[i] % mod;
F.mul(F.f, F.g, F.mx); F.clear(); reverse(F.f, F.f + F.n + 1);
for(int i = 0; i <= k; ++i) F.f[i] = 1ll * F.f[i] * inv[i] % mod, ans = add(ans, F.f[i]);
printf("%d\n", ans);
return 0;
}
接下来再说另一种做法，也就是本文最开始就提到的做法：
\[ans=\frac{n!}{2^D}\times[x^n]e^{-Dx}\times\sum_{k=0}^{mx}\binom{D}{k}\times (e^{2x}-1)^k\times (e^{2x}+1)^{D-k}\\ \]其实主要就是对于 \(\sum_{k=0}^{mx}\binom{D}{k}\times (e^{2x}-1)^k\times (e^{2x}+1)^{D-k}\) 无法化简，展开。现在我们设 \(G(x)=\sum_{k=0}^{mx}\binom{D}{k}(x-1)^k\times (x+1)^{D-k}\) ，那么只要我们知道了 \(G(x)\) 的各项系数，这个式子就是可做的， \(ans\) 可化为：
\[\begin{aligned} ans&=\frac{n!}{2^D}\times[x^n]e^{-Dx}\times G(e^{2x})\\ &=\frac{n!}{2^D}\times[x^n]\sum_{i=0}^{D}g_i\times (e^{2x})^i\times e^{-Dx}\\ &=\frac{n!}{2^D}\times[x^n]\sum_{i=0}^{D}g_i\times e^{(2i-D)x}\\ &=\frac{n!}{2^D}\times\sum_{i=0}^{D}g_i\times \frac{(2i-D)^n}{n!}\\ &=\frac{1}{2^D}\times\sum_{i=0}^{D}g_i\times (2i-D)^n \end{aligned} \]这一步可以 \(O(D+\frac{D}{\log D}\times \log n)=O(D\log_D n)\) 通过欧拉筛筛 \(n\) 次幂求解。
所以现在的问题就是如何求 \(G(x)\) 的系数。这一步可以考虑求导。
\[\begin{aligned} G'(x)&=\sum_{k=0}^{mx}\frac{D}{k!(D-k)!}\times (k\times (x-1)^{k-1}\times(x+1)^{D-k}+(D-k)\times (x-1)^k\times (x+1)^{D-k-1})\\ &=\sum_{k=0}^{mx}\frac{D\times (D-1)!}{k!(D-k)!}\times k(x-1)^{k-1}(x+1)^{D-k}+\frac{D\times (D-1)!}{k!(D-k)!}(D-k)(x-1)^k(x+1)^{D-k-1}\\ &=\sum_{k=0}^{mx}D\times (\frac{(D-1)!}{(k-1)!(D-k)!}\times (x-1)^{k-1}(x+1)^{D-k}+\frac{(D-1)!}{k!(D-k-1)!}(x-1)^k(x+1)^{D-k-1})\\ &=\sum_{k=0}^{mx}D\times (\binom{D-1}{k-1}\times (x-1)^{k-1}(x+1)^{D-k}+\binom{D-1}{D-k-1}(x-1)^k(x+1)^{D-k-1})\\ \end{aligned} \]恰好 \(\binom{D-1}{k-1}+\binom{D-1}{D-k-1}=\binom{D}{k}\) ，所以尝试将 \(G(x)\) 变一下形：
\[G(x)=\sum_{k=0}^{mx}(\binom{D-1}{k-1}(x-1)^k\times (x+1)^{D-k}+\binom{D-1}{D-k-1}(x-1)^k\times (x+1)^{D-k}) \]这样 \(G(x)\) 和 \(G'(x)\) 形式就相近了，然后可以作差之类的找到 \(G(x)\) 与 \(G'(x)\) 之间的系数关系，也即找到了 \(G(x)\) 的 \(x^i\) 项和 \(x^{i-1}\) 项之间的关系，所以可以凑一下式子， \(G(x)\) 差一个常数 \(D\) ， \(G'(x)\) 次数低了 \(1\) ，于是考虑用 \(D\times G(x)-x\times G'(x)\) ：
\[\begin{aligned} DG(x)-xG'(x)&=\sum_{k=0}^{mx}D\times (-\binom{D-1}{k-1}(x-1)^{k-1}(x+1)^{D-k}+\binom{D-1}{D-k-1}(x-1)^{k}(x+1)^{D-k-1})\\ &=D\times (-\sum_{k=0}^{mx}\binom{D-1}{k-1}(x-1)^{k-1}(x+1)^{D-k}+\sum_{k=0}^{mx}\binom{D-1}{k}(x-1)^{k}(x+1)^{D-k-1})\\ &=D\times (\sum_{k=1}^{mx+1}\binom{D-1}{k-1}(x-1)^{k-1}(x+1)^{D-k}-\sum_{k=0}^{mx}\binom{D-1}{k-1}(x-1)^{k-1}(x+1)^{D-k})\\ &=D\times \binom{D-1}{mx}(x-1)^{mx}(x+1)^{D-mx-1} \end{aligned} \]这样就会变得简洁很多，并且由于 \([x^i](DG(x)-xG'(x))=(D-n)[x^i]G(x)\) ，所以可以通过求 \(D\times \binom{D-1}{mx}(x-1)^{mx}(x+1)^{D-mx-1}\) 的系数来得出 \(G(x)\) 的系数，这个可以直接卷积，也可以再次求导降低复杂度 ，下面是求系数方法：
令 \(H(x)=(x-1)^a(x+1)^b\) ，则 ：
\[\begin{aligned} &H'(x)=a(x-1)^{a-1}(x+1)^{b}+b(x-1)^{a}(x+1)^{b-1}\\ &\frac{H'(x)}{H(x)}=\frac{a}{x-1}+\frac{b}{x+1}\\ &x^2H'(x)-H'(x)=(a+b)xH(x)+(a-b)H(x)\\ &(i-1)h_{i-1}-(i+1)h_{i+1}=(a+b)h_{i-1}+(a-b)h_i\\ &h_{i+1}=-\frac{1}{i+1}((a+b-i+1)h_{i-1}+(a-b)h_i) \end{aligned} \]这样，令 \(a=mx,b=D-mx-1\) 然后就可以递推出 \(D\times \binom{D-1}{mx}(x-1)^{mx}(x+1)^{D-mx-1}\) 的系数，但有几个要注意的地方：
1.根据 \(H(x)=(x-1)^a(x+1)^b\) 这个式子可知 \(h_0=(-1)^a,h_1=(-1)^a\times b+(-1)^{a-1}\times a\) ；
2.观察递推式可以发现 \(h_D\) 通过递推得到为 \(0\) ，但显然 \(h_D\) 不为 \(0\) ，所以要单独特殊计算，即 \(h_D=\sum_{i=0}^{a}\binom{D}{i}\) 。
递推是 \(O(D)\) 的，所以总的复杂度为 \(O(D\log_D n)\) 。代码如下：
#include <bits/stdc++.h>
using namespace std;
const int N = 3e5 + 30;
const int M = 3e5;
const int mod = 998244353;
const int pr = 3;
const int ig = 332748118;
inline int add(int x, int y) {
x += y;
return x >= mod ? x - mod : x;
}
inline int del(int x, int y) {
x -= y;
return x < 0 ? x + mod : x;
}
inline int read() {
char ch = getchar(); int x = 0;
while(!isdigit(ch)) ch = getchar();
while(isdigit(ch)) {x = (x << 3) + (x << 1) + ch - 48; ch = getchar();}
return x;
}
int qpow(int a, int b) {
int res = 1;
while(b) {
if(b & 1) res = 1ll * res * a % mod;
a = 1ll * a * a % mod;
b >>= 1;
}
return res;
}
int D, n, m, k, fac[N], inv[N], ans, f[N], g[N], a, b, mem, Inv[N], prime[N], pw[N], opt = 1;
bool vis[N];
inline int C(int x, int y) {return 1ll * fac[x] * inv[y] % mod * inv[x - y] % mod;}
void init() {
pw[1] = 1; pw[0] = 0;
for(int i = 2; i <= D; ++i) {
if(!vis[i]) {
prime[++prime[0]] = i;
pw[i] = qpow(i, n);
}
for(int j = 1; j <= prime[0] && i * prime[j] <= D; ++j) {
vis[i * prime[j]] = true;
pw[i * prime[j]] = 1ll * pw[i] * pw[prime[j]] % mod;
if(i % prime[j] == 0) break;
}
}
}
int main() {
D = read(); n = read(); m = read(); k = min(n - 2 * m, D);
if(n & 1) opt = mod - 1;
if(k == D) return printf("%d\n", qpow(D, n)), 0;
if(k < 0) return puts("0"), 0;
fac[0] = 1; Inv[1] = 1; init();
for(int i = 1; i <= D; ++i) fac[i] = 1ll * fac[i - 1] * i % mod;
inv[D] = qpow(fac[D], mod - 2);
for(int i = D; i >= 1; --i) inv[i - 1] = 1ll * inv[i] * i % mod;
for(int i = 2; i <= D; ++i) Inv[i] = 1ll * (mod - mod / i) * Inv[mod % i] % mod;
a = D - k - 1; b = k; mem = C(D - 1, k);
f[0] = (b & 1) ? (mod - 1) : 1;
f[1] = add(1ll * ((b & 1) ? mod - 1 : 1) * a % mod, 1ll * (((b - 1) & 1) ? mod - 1 : 1) * b % mod);
for(int i = 1; i < D; ++i)
f[i + 1] = 1ll * (mod - 1) * del(add(1ll * a * del(f[i - 1], f[i]) % mod, 1ll * b * add(f[i - 1], f[i]) % mod), 1ll * (i - 1) * f[i - 1] % mod) % mod * Inv[i + 1] % mod;
for(int i = 0; i < D; ++i)
g[i] = 1ll * f[i] * D % mod * mem % mod * Inv[D - i] % mod;
for(int i = 0; i <= k; ++i)
g[D] = add(g[D], C(D, i));
for(int i = 0; i <= D; ++i)
if(2 * i >= D) ans = add(ans, 1ll * g[i] * pw[2 * i - D] % mod) % mod;
else ans = add(ans, 1ll * g[i] * pw[D - 2 * i] % mod * opt % mod);
ans = 1ll * ans * qpow(2, mod - 1 - D) % mod;
printf("%d\n", ans);
return 0;
}