[CF1477F] Nezzar and Chocolate Bars
Problem
给定 \(n\) 个长度为 \(a_i\) 的巧克力，每次以正比于 \(a_i\) 的概率取得一个巧克力，然后在 \((0,a_i)\) 中随机选择一个实数 \(r\) 并将其分成 \(r,a_i-r\) 两个部分放回。
计算使得所有巧克力的长度均小于 \(K\) 的期望操作次数。
\(n\le 50,k,\sum a_i\le 2000\)。
Sol
这里将每块巧克力看成线段，长度为 \(L_i\)，线段长度总和为 \(L\)。
n=1
首先，我们来考虑 \(n=1\) 的情形。
问题变成：给你一条长度为 \(L\)，每次等概率在线段上随机一个点切断，求期望次数，使得第一次所有线段长度 \(\le K\)。
现在，假设切 \(n\) 次得到 \(n\) 个点 \(X_{(1)},X_{(2)},\cdots,X_{(n)}\)，并且按照从小到大排序得到 \(X_1\le X_2\le\cdots X_n\)。令 \(X_0=0\)，按题意，必须满足
令 \(w=\frac KL\)，设 \(z_1=X_1-X_0\)，\(z_2=X_2-X_1\)，\(\cdots\)，\(z_n=X_n-X_1\)，我们要求概率
考虑求后面的积分，发现其可以差分，记
则之前的积分可以写成
请注意：这个差分当 \(\frac 1w=1\) 时需要修正，直接导致后面代码里的特判。
现在我们来求 \(F(x)\)。其实 \(\sum_{i=1}^nz_i\) 遵循 Irwin-Hall 分布，其 累积分布函数 满足
我们从 OI 的角度试图不严谨地证明它。首先假设 \(0\le x\le1\)，当 \(n=1\) 时，\(F(x)=x\)；当 \(n=2\) 时，\(F(x)\) 是二维平面坐标系中的直角三角形面积，为 \(\frac{x^2}2\)；当 \(n=3\) 时，是直角四面体体积，为 \(\frac{x^3}6\cdots\) 不难发现，当 \(x\le 1\) 时，\(F(x)=\frac{x^n}{n!}\)。当 \(x\ge 1\) 时，我们可以通过容斥，每次枚举哪些 \(z_i\ge 1\)，不断扣掉，直接推出来上式成立。
记 \(p_n\) 表示带回原式
这里 \(x=\lfloor\frac1w\rfloor\)。
最终答案根据期望定义，有
利用公式
即可求解。
n 任意
这里定义第 \(i\) 条线段，切 \(m\) 次得到结果非法的概率为 \(q_{i,m}\)，由上面得到
设 \(p_m\) 表示切 \(m\) 刀 不合法 的概率，则有
记 \(Q_i(x)\) 表示后面部分的 EGF，则
记 \(P(x)\) 表示 \(p_n\) 的 EGF，则
现在详细推导 \(Q_i(x)\)
由于
所以有
求出来的 \(P(x)\) 是 EGF，而设 \(P^*(x)\) 是 \(p_m\) 的 OGF。通过转化，我们用 \(\sum_{t,r}a_{t,r}x^{t-r}\exp\left(\left(1-\dfrac KLt\right)x\right)\) 的形式表示 \(P(x)\)，其中 \(0\le t\le \frac LK, 0\le r\le\min(n,t)\)。现在，我们想要 \(P(x)\rightarrow P^*(x)\)，就要着手考虑 \(x^k\exp(Cx)\) 的形式，不难发现
转变成 OGF：
通过这样我们就能得到 \(P^*(x)\)。最后答案为 \(P^*(1)\)，直接代值即可。
复杂度 \(\mathcal O(nL^2)\)，使用 NTT 可使复杂度降成 \(\mathcal O(nL\log n\log nL)\)。
#include <bits/stdc++.h>
#define pb push_back
using std::vector;
typedef vector<int> poly;
typedef vector<poly> poly2;
const int N = 55, M = 2005, P = 998244353;
int n, K, L, l[N], len;
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[M * 2], inv[M], fac[M], ifac[M];
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 < i; j++)
W[i + j] = 1LL * W[i + j - 1] * Wn % P;
inv[1] = fac[0] = ifac[0] = 1;
for (int i = 2; i <= n; i++)
inv[i] = 1LL * (P - P / i) * inv[P % i] % P;
for (int i = 1; i <= n; i++)
fac[i] = 1LL * fac[i - 1] * i % P, ifac[i] = 1LL * ifac[i - 1] * inv[i] % P;
}
void fft(poly &a, int n, int op) {
a.resize(n);
static int rev[M * 4] = {0};
for (int i = 1; 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.begin() + 1, a.end());
}
int getsz(int n) { int x = 1; while (x < n) x <<= 1; return x; }
struct Poly {
poly2 a; int n, m;
poly &operator [] (int i) { return a[i]; }
Poly(poly2 &po, int nn, int mm) : a(po), n(nn), m(mm) {}
};
vector<Poly> F;
Poly operator * (Poly a, Poly b) {
poly2 res;
int n = a.n + b.n - 1, m = a.m + b.m - 1, ln = getsz(n), lm = getsz(m);
for (int i = 0; i < a.n; i++) fft(a[i], lm, 1);
for (int i = 0; i < b.n; i++) fft(b[i], lm, 1);
res.resize(n);
for (int i = 0; i < lm; i++) {
poly l, r;
for (int j = 0; j < a.n; j++) l.pb(a[j][i]);
for (int j = 0; j < b.n; j++) r.pb(b[j][i]);
fft(l, ln, 1), fft(r, ln, 1);
for (int j = 0; j < ln; j++) l[j] = 1LL * l[j] * r[j] % P;
fft(l, ln, 0);
for (int j = 0; j < n; j++) res[j].pb(l[j]);
}
for (int i = 0; i < n; i++) fft(res[i], lm, 0), res[i].resize(m);
return Poly(res, n, m);
}
int main() {
scanf("%d%d", &n, &K);
for (int i = 1; i <= n; i++)
scanf("%d", &l[i]), L += l[i];
len = L + 1; prework(len);
for (int i = 1; i <= n; i++) {
poly tmp0, tmp1; poly2 tmp;
for (int j = 0; j <= l[i] / K; j++) {
int t0, t1;
t0 = 1LL * (j & 1 ? P - 1 : 1) * ifac[j] % P * qpow(1LL * (l[i] - j * K) * inv[L] % P, j) % P;
if (j) t1 = 1LL * (j & 1 ? P - 1 : 1) * ifac[j - 1] % P * qpow(1LL * (l[i] - j * K) * inv[L] % P, j - 1) % P; else t1 = 0;
tmp0.pb(t0), tmp1.pb(t1);
}
if (l[i] == K) tmp1[1] = inc(tmp1[1], 1); // 由于 F(1/w)-F(0) 出错导致的修正
F.pb(Poly(tmp = {tmp0, tmp1}, 2, l[i] / K + 1));
}
for (int i = 1; i < n; i <<= 1)
for (int j = 0; j < n - i; j += i * 2)
F[j] = F[j] * F[j + i];
int ans = 0;
for (int i = 0; i < F[0].n; i++)
for (int j = 0; j < F[0].m; j++) {
if (j < i) continue;
ans = (ans + 1LL * fac[j - i] * qpow(L, j - i + 1) % P * qpow(qpow(1LL * j * K % P, P - 2), j - i + 1) % P * F[0][i][j]) % P;
}
printf("%d\n", ans ? P - ans : 0);
return 0;
}
