Problem
有 (n) 个范围在 ([1,D]) 内的整数均匀随机变量，求至少能选出 (m) 个瓶子，使得存在一种方案，选择一些变量，并把选出来的每一个变量放到一个瓶子中，满足每个瓶子都恰好装两个值相同的变量的概率。答案乘上 (D^n) 对 (998244353) 取模。
(Dle 10^5)，(n,mle 10^9)。
Sol
首先根据奇偶性，写出来关于每种颜色的多项式（(x) 表示选出来的数目，(y) 表示能组合的 pair ( imes 2)）
[F(x,y)=frac{e^{xy}+e^{-xy}}2+frac{e^{xy}-e^{-xy}}{2y}
]
答案为
[Ans=n![x^ny^{ge2m}]F^D(x,y)
]
[Rightarrow n![x^ny^{ge2m-n}]left(frac{e^x+e^{-x}}2+frac{e^x-e^{-x}}{2y}
ight)^D
]
[Rightarrow n![x^ny^{ge 2m-n+D}]left(frac{e^x+e^{-x}}2y+frac{e^x-e^{-x}}2
ight)^D
]
无必要的 (y) 消掉了。
展开右式
[sum_{i=0}^D{Dchoose i}y^ileft(frac{e^x+e^{-x}}2
ight)^ileft(frac{e^x-e^{-x}}2
ight)^{D-i}
]
答案变成
[Ans=n![x^n]sum_{i=a}^D{Dchoose i}left(frac{e^x+e^{-x}}2
ight)^ileft(frac{e^x-e^{-x}}2
ight)^{D-i}
]
这里 (a=2m-n+D)。
继续考察右式
[egin{aligned}
&sum_{i=a}^D{Dchoose i}left(frac{e^x+e^{-x}}2
ight)^ileft(frac{e^x-e^{-x}}2
ight)^{D-i}\
=&2^{-D}e^{-Dx}sum_{i=a}^D{Dchoose i}(e^{2x}+1)^i(e^{2x}-1)^{D-i}
end{aligned}
]
这里重点来看后半部分。暂时用 (x) 替换 (e^{2x})。
[egin{aligned}
sum_{i=a}^D{Dchoose i}(x+1)^i(x-1)^{D-i}
end{aligned}
]
它一定可以变成
[sum_{i=0}^Da_ix^i
]
这样子就可以解决问题了，现在我们来尝试化简。这里令 (z=x+1)，则式子变成
[egin{aligned}
&sum_{i=a}^D{Dchoose i}z^i(z-2)^{D-i}\
=&sum_{i=a}^D{Dchoose i}sum_{j=0}^{D-i}{D-ichoose j}(-2)^{D-i-j}z^{i+j}\
=&sum_{i=a}^D(-2)^{D-i}z^isum_{j=0}^{i-a}{Dchoose i-j}{D-(i-j)choose j}\
=&sum_{i=a}^D(-2)^{D-i}z^i{Dchoose i}sum_{j=0}^{i-a}{ichoose j}
end{aligned}
]
发现 (sum_{j=0}^{i-a}{ichoose j}) 是可以递推的！于是我们可以在 (mathcal O(D)) 的时间求出上面的多项式，对 (z) 做平移得到原来的多项式，问题得解。最后得到
[Ans=n![x^n]sum_{i=0}^Df_ie^{(2i-D)x}
]
直接解系数即可。复杂度 (mathcal O(Dlog n))。
附平移多项式推导
[egin{aligned}
f(x+c)&=sum_{i=0}^na_i(x+c)^i\
&=sum_{i=0}^na_isum_{j=0}^i{ichoose j}x^jc^{i-j}\
&=sum_{j=0}^nfrac{x^j}{j!}sum_{i=j}^ni!a_ifrac{c^{i-j}}{(i-j)!}
end{aligned}
]
NTT 即可。
Code
#include <bits/stdc++.h>
using std::max;
const int N = 100005, P = 998244353;
int D, n, m, a[N * 4], b[N * 4], s;
int inc(int a, int b) { return (a += b) >= P ? a - P : a; }
int qpow(int a, int b) {
int t = 1;
for (; b; b >>= 1, a = 1LL * a * a % P)
if (b & 1) t = 1LL * t * a % P;
return t;
}
int W[N * 4], inv[N], fac[N], ifac[N];
int C(int n, int m) {
if (n < 0 || m > n) return 0;
return 1LL * fac[n] * ifac[m] % P * ifac[n - m] % P;
}
void prework(int n) {
for (int i = 1; i < n; i <<= 1)
for (int j = W[i] = 1, Wn = qpow(3, (P - 1) / i >> 1); j < n; j++)
W[i + j] = 1LL * W[i + j - 1] * Wn % P;
}
void fft(int *a, int n, int op) {
static int rev[N * 4] = {0};
for (int i = 0; i < n; i++)
if ((rev[i] = rev[i >> 1] >> 1 | (i & 1 ? n >> 1 : 0)) > i) std::swap(a[i], a[rev[i]]);
for (int q = 1; q < n; q <<= 1)
for (int p = 0; p < n; p += q << 1)
for (int i = 0, t; i < q; i++)
t = 1LL * W[q+i] * a[p+q+i] % P, a[p+q+i] = inc(a[p+i], P-t), a[p+i] = inc(a[p+i], t);
if (op) return;
for (int i = 0, inv = qpow(n, P - 2); i < n; i++)
a[i] = 1LL * a[i] * inv % P;
std::reverse(a + 1, a + n);
}
int getsz(int n) { int x = 1; while (x < n) x <<= 1; return x; }
int main() {
scanf("%d%d%d", &D, &n, &m);
prework((D + 1) * 2);
inv[1] = fac[0] = ifac[0] = 1;
for (int i = 2; i <= D; i++)
inv[i] = 1LL * (P - P / i) * inv[P % i] % P;
for (int i = 1; i <= D; i++)
fac[i] = 1LL * fac[i - 1] * i % P, ifac[i] = 1LL * ifac[i - 1] * inv[i] % P;
s = 2 * m - n + D;
int tmp = 1;
for (int i = max(0, s); i <= D; i++) {
a[i] = 1LL * qpow(P - 2, D - i) * C(D, i) % P * tmp % P * fac[i] % P;
tmp = (2LL * tmp + C(i, i - s + 1)) % P;
}
for (int i = 0; i <= D; i++)
b[i] = ifac[i];
std::reverse(b, b + D + 1);
int l = getsz(2 * (D + 1));
fft(a, l, 1), fft(b, l, 1);
for (int i = 0; i < l; i++) a[i] = 1LL * a[i] * b[i] % P;
fft(a, l, 0);
int ans = 0;
for (int i = 0; i <= D; i++) {
a[i + D] = 1LL * a[i + D] * ifac[i] % P;
ans = (ans + 1LL * a[i + D] * qpow((2 * i - D + P) % P, n)) % P;
}
ans = 1LL * ans * qpow((P + 1) / 2, D) % P;
printf("%d
", ans);
return 0;
}