CF1477F

CF1477F [* easy]

给定 \(n\) 个长度为 \(a_i\) 的巧克力，每次以正比于 \(a_i\) 的概率取得一个巧克力，然后在 \((0,a_i)\) 中随机选择一个实数 \(r\) 并将其分成 \(r,a_i-r\) 两个部分放回。

计算使得所有巧克力的长度均小于 \(k\) 的期望操作次数。

\(n\le 50,\sum a_i,k\le 2000\)

---

对于 \(n=1\) 的情况简单推导可以给出答案的算式（令 \(m=a_1\)）：

\[\sum_{j\ge 1}^{m/k} (\frac{m-jk}{m})^{j-1}(-1)^{j-1} (\frac{m}{jk})^{j+1} \]

具体的：

仅考虑 \(n=1\)，可以发现题目描述的概率可以简单的描述为在长度为 \(m\) 的数轴上切割若干刀，使得相邻两个间隙的距离最大值小于等于 \(k\) 的期望刀数。

一般的，我们转换为期望并计数，答案可以写为 \(\sum P(X\ge i)\)，此条件下有这个问题与 \(i+1\) 个随机均匀实数变量 \(\{x_j\}\)，满足 \(\sum x_j=m\)，其中 \(\max x_i\ge k\) 的概率是一致，转换为 \(1\) 减去 \(x_i\in [0,k]\) 的概率即可。

这个部分只需要考虑 \(n\) 个随机变量满足 \(\sum x_i=m\) 的 “超体积” （不知道准不准确，我大概是这样理解的/yun）即可，是 \(\frac{m^{n-1}}{(n-1)!}\)，对于 \([0,k]\) 的限制容斥掉，不难给出答案的算式：

\[\begin{aligned} &\sum_{i}\bigg(1-\sum_{jk\le m} \binom{i}{j}(-1)^{j}\frac{(m-jk)^{i-1}/(i-1)!}{m^{i-1}/(i-1)!}\bigg) \\&=\sum_{j\ge 1}^{jk\le m} \sum_{i\ge j} \binom{i}{j}(-1)^{j-1}(\frac{m-jk}{m})^{i-1} \end{aligned}\]

考虑到 \(\frac{1}{(1-x)^{j+1}}\) 的展开形式为 \(\sum \binom{i+j}{j} x^i\)，我们可以给出答案为：

\[\sum_{j\ge 1}^{jk\le m} (-1)^{j-1}(\frac{m-jk}{m})^{j-1} \frac{1}{(1-\frac{m-jk}{m})^{j+1}} \]

---

考虑 \(n\ne 1\) 的情况，一般化的，我们仍然考虑枚举 \(i\) 并计算 \(P(X\ge i)\)

答案就是：

\[\sum (1-P(X\ge i)) \]

我们假定每个段分别被操作了 \(i_1,i_2...i_m\) 次，最终被操作了 \(i=i_1+...i_m\) 次，此时对应的概率为 \(\frac{i!}{\prod i_j! (\frac{a_j}{S})^{i_j-1}}\)

对于下半部分，我们考虑对于每个 \(i\) 建立一个关于 \(a_i\) 的指数型生成函数，我们先假定我们将无穷项都存储了起来，那么接下来只需要将这些多项式乘起来即可得到 \(P(X\ge i)\)，于是方便起见设 \(m=a_i\)，此处与 \(n=1\) 的时候相似：

\[\sum f_i \frac{(\frac{m}{S}x)^i}{i!} \]

其中 \(f_i\) 表示操作了 \(i\) 次（对应着 \(i+1\) 个随机变量的情况）但最大值小于等于 \(k\) 的概率。

\[\begin{aligned} &\sum_{i} \frac{(\frac{m}{S}x)^i}{i!}\sum_{jk\le m} \binom{i+1}{j}(-1)^j (\frac{m-jk}{m})^{i} \\&=\sum_{jk\le m} (-1)^j\sum_{i} \frac{(\frac{m-jk}{S}x)^i}{i!}\binom{i+1}{j} \\&=\sum_{jk\le m} (-1)^j\frac{1}{j!}\sum_{i} \frac{(i+1)(\frac{m-jk}{S}x)^i}{(i+1-j)!} \\&=\sum_{jk\le m} (-1)^j\frac{1}{j!}\frac{S}{m-jk}((\frac{m-jk}{S}x)^{j}e^{(\frac{m-jk}{S}x)})' \\&=\sum_{jk\le m} (-1)^j\frac{1}{j!}\bigg(j((\frac{m-jk}{S})x)^{j-1}e^{(\frac{m-jk}{S}x)}+(\frac{m-jk}{S})^{j}x^{j}e^{(\frac{m-jk}{S}x)}\bigg) \end{aligned}\]

对于每个 \(i\)，记上述多项式为 \(G_i(x)\)，令

\[F(x)=\prod G_i(x) \]

答案为：

\[\sum_{i} (1-F(x)[x^i]i!) \]

我们可以观察到答案中必然存在一个 \(e^x\) 的项数，其展开形式恰好抵消了 \(+1\) 的影响。

我们考虑 \(f_i x^ie^{zx}\) 在求和式中的贡献可以写为 \(f_i\sum \frac{(i+j)!}{j!}z^j\)，等价于 \(f_i\times i!\sum \binom{i+j}{j}z^j\) 这与 \(n=1\) 的形式是一致的，即 \(f_i\times i!\times \frac{1}{(1-z)^{i+1}}\)

其次，可以观察到 \(j\) 与 \(\frac{m-jk}{S}\) 中的 \(j\) 差值至多为 \(1\)（来源于 \(x^{j-1}e^{\frac{m-jk}{S}}\) 这一项），总差值不超过 \(n\) 对此记录并 dp 可以做到 \(\mathcal O(\frac{nL^2}{k})\)

可以使用 NTT 优化至 \(\mathcal O(nL\log L)\) 。

tips：可以注意到 \(e\) 的指标一定是 \(e^{1-\frac{jk}{S}}\)，提取系数时的答案刚好是 \(\frac{S^{i+1}}{(jk)^{i+1}}\)

\(Code:\)

#include<bits/stdc++.h>
using namespace std ;
#define Next( i, x ) for( register int i = head[x]; i; i = e[i].next )
#define rep( i, s, t ) for( register int i = (s); i <= (t); ++ i )
#define drep( i, s, t ) for( register int i = (t); i >= (s); -- i )
#define re register
#define mp make_pair
#define pi pair<int, int>
#define pb push_back
#define vi vector<int>
#define int long long
int gi() {
	char cc = getchar() ; int cn = 0, flus = 1 ;
	while( cc < '0' || cc > '9' ) {  if( cc == '-' ) flus = - flus ; cc = getchar() ; }
	while( cc >= '0' && cc <= '9' )  cn = cn * 10 + cc - '0', cc = getchar() ;
	return cn * flus ;
}
const int N = 52 ; 
const int M = 4000 + 5 ; 
const int P = 998244353 ;
int n, m, L, S, a[N], f[N][M], g[N][M] ; 
int inv[M], fac[M], ifac[M] ; 
int fpow(int x, int k) {
	int ans = 1, base = x ;
	while(k) {
		if(k & 1) ans = 1ll * ans * base % P ;
		base = 1ll * base * base % P, k >>= 1 ;
	} return ans ;
}
void inc(int &x, int y) { ((x += y) >= P) && (x -= P) ; }
void dec(int &x, int y) { ((x -= y) < 0) && (x += P) ; }
signed main()
{
	n = gi(), m = gi(), fac[0] = 1 ; 
	rep( i, 1, 2000 ) inv[i] = fpow(i, P - 2), fac[i] = fac[i - 1] * i % P ; 
	rep( i, 0, 2000 ) ifac[i] = fpow(fac[i], P - 2) ; 
	rep( i, 1, n ) a[i] = gi(), S += a[i] ; 
	f[0][0] = 1 ; 
	rep( i, 1, n ) {
		int lim = (a[i] / m) ; L += lim ; 
		rep( j, 0, i ) rep( k, 0, L ) g[j][k] = f[j][k], f[j][k] = 0 ; 
		rep( j, 0, lim ) {
			int z = a[i] - j * m, o = z * inv[S] % P, v = ifac[j] * fpow(P - 1, j) % P ; 
			int zz = fpow(o, j + P - 2) * j % P * v % P ;
			rep( k, 0, i ) rep( l, 0, L ) 
			inc(f[k + 1][l + j], g[k][l] * zz % P) ; 
			int zt = fpow(o, j) * v % P ; 
			rep( k, 0, i ) rep( l, 0, L ) 
			inc( f[k][l + j], g[k][l] * zt % P) ; 
		}
	} int ans = 0 ; 
	rep( o, 0, n ) rep( j, 0, L ) {
		if((!o) && (!j)) continue ; 
		if(!f[o][j]) continue ; 
		int i = j - o ; 
		int d = f[o][j] * fac[i] % P * fpow(S * fpow(j * m, P - 2) % P, i + 1) % P ; 
		dec(ans, d) ; 
	}
	cout << ans << endl ; 
	return 0 ;
}