[BZFZ友谊赛]火山喷发

题目

传送门 to CF：这是 $\mathrm{C}\mathrm{F}\_\mathrm{G}\mathrm{Y}\mathrm{M}102798\mathrm{E}$ 一道原题。

题目描述
有 $n$ 个随从，第 $i$ 个的生命值为 $a_i$ 。现在你会连续攻击 $m$ 次，每次等概率地选择一个仍然存活的随从，并对它造成 $1$ 的伤害（使其生命值减小 $1$）。当一个随从的生命值降为 $0$ 时，它立刻死亡（不再存活）。

现在问你，死亡随从数量的期望是多少。对 $998244353$ 取模。

数据范围与提示
$n\le 15,m\le 200,a_i\le 200$，保证 $m\le {\sum }_i{a}_i$ 。

思路壹

我考场上就想的这种做法。但是因为很难打，并且复杂度不太对，就放弃了……

出师不利

经典题目 猎人杀，大家肯定都做过。还是用这种方法：认为死掉的随从仍然可以被攻击，只是它无效而已。

不难发现，这里有两个要素：有效攻击次数；总攻击次数。前者是满足题目条件用的，后者是计算概率用的。

那么我们可以很轻松地写出一个二元生成函数，用 $x$ 代表有效攻击，$y$ 代表总攻击：
$$F_i(x,y)=\sum _{b\ge 0}{x}^{\min (a_i,b)}\cdot \frac{y^b}{b!}\cdot {\left(\frac{1}{n}\right)}^b$$

显然 $y$ 应当是指数型，因为 $y$ 的合并是需要用组合数的。它计算的是排列。

可以稍微化简一下，得到
$$F_i(x,y)=x^{a_i}\left(e^{y/n}-\sum _{j<a_i}\frac{y^j}{n^j\cdot j!}\right)+\sum _{j<a_i}\frac{x^j{y}^j}{n^j\cdot j!}$$

二元生成函数的卷积，大家都会做，就是 $f_{aM+b}$ 为 $x^a{y}^b$ 的系数，然后卷积，其中 $\lfloor \frac{M-1}{2}\rfloor$ 是 $y$ 的次数最大值。考虑到 $x$ 的范围是 $m$ 而 $y$ 的范围是 $\sum a_i\le na$，复杂度是
$$\mathcal{O}(nma\log nma)$$

那么，我们只需要尽力求出
$$P(x,y)=\prod _i{F}_i(x,y)=\sum _{j=0}^n{A}_j(x,y)\cdot e^{yj/n}$$

然后我们求 ${\sum }_{j=0}^{+\infty }[x^m{y}^j]j!\cdot P(x,y)$ 就毫不费力了。毕竟 $e^{yj/n}={\sum }_{t=0}^{+\infty }\frac{(yj{)}^t}{n^t\cdot t!}$ 嘛。枚举 $A_j(x,y)$ 中 $y$ 的次数 $k$，乘 $e^{yj/n}$ 得 ${\sum }_{t=0}^{+\infty }\frac{j^t{y}^{t+k}}{n^t\cdot t!}$，第 $j$ 项乘 $j!$ 求和得到
$$\sum _{t=0}^{+\infty }\frac{j^t\cdot (t+k)!}{n^t\cdot t!}$$

我们有理由相信，这是可以 $\mathcal{O}(\log k)$ 计算的。那么统计答案的总复杂度就是 $\mathcal{O}(n^{2}a\log na)$ 了。

怎么求 $A_j(x,y)$ 呢？考虑直接 $\mathtt{N}\mathtt{T}\mathtt{T}$ 的话，要计算 $n^2$ 组二元卷积，所以复杂度是
$$\mathcal{O}(n^{3}ma\log nma)$$

这是概率，现在要求期望，只需要钦定某个随从必然被干掉。即只取 $b\ge a_i$ 的情况。用算出的 $P(x,y)$ 减去后面的 ${\sum }_{j<a_i}\frac{x^j{y}^j}{n^j\cdot j!}$ 就行了。这个就预处理前缀、后缀，然后三者相乘。

容易发现，复杂度已经受不了了……

工业革命

原因在于 $n,m$ 并不对等，二者数据范围相距悬殊。这么小的 $n$，肯定引导我们往 $2^n\sim 3^n$ 左右去思考。

考虑枚举 $b\ge a_i$ 的集合 $S$，因为此时 ${\sum }_{j<a_i}\frac{x^j{y}^j}{n^j\cdot j!}$ 满足 $x,y$ 次数相同，就不需要使用二元生成函数了。一元生成函数的卷积，复杂度就可以明显优化下来。

具体而言，我们只需要求
$$[y^{m-{\sum }_{j\in S}{a}_j}]\prod _{j\notin S}\left(\sum _{b<a_j}\frac{y^b}{n^b\cdot b!}\right)$$

这是被 $x$ 的次数限制住了的。然后再求
$$\prod _{j\in S}\left(e^{y/n}-\sum _{b<a_j}\frac{y^b}{n^b\cdot b!}\right)$$

注意这里 $y$ 的次数，有一个很强劲的限制：它是钦定了被干掉的随从，所以 $\sum a\le m$，所以 $y$ 的最高次数是 $m$ 而不是 $na$ 。

但是我们有一个很高妙的做法：要知道，系数为 $1$ 的指数型生成函数，求导 是不会大变样的。我们试试，对 $y$ 求导：
$${\left(e^{y/n}-\sum _{b<a_j}\frac{y^b}{n^b\cdot b!}\right)}^{\prime }=\frac{e^{y/n}}{n}-\sum _{b=1}^{a_j-1}\frac{y^{b-1}}{n^b\cdot (b-1)!}$$

记左式（原式）为 $F_j(y)$，右式则为 $F_j^{\prime }(y)$，不难发现
$$\begin{gathered} F_j(y)=n\cdot F_j^{\prime }(y)-\frac{y^{a_j-1}}{n^{a_j-1}\cdot (a_j-1)!} \\ \Rightarrow F_j^{\prime }(y)=\frac{F_j(y)}{n}+\frac{y^{a_j-1}}{n^{a_j}\cdot (a_j-1)!} \end{gathered}$$

那么，我们将目标式子 ${\prod }_{j\in S}{F}_j(y)$ 求导试一试：
$$\begin{array}{rl}& {\left[\prod _{j\in S}{F}_j(y)\right]}^{\prime }=\sum _{x\in S}{F}_x^{\prime }(y)\prod _{j\in S}^{j\ne x}{F}_j(y)\\ & =\sum _{x\in S}\left[\frac{F_x(y)}{n}+\frac{y^{a_x-1}}{n^{a_x}\cdot (a_x-1)!}\right]\prod _{j\in S}^{j\ne x}{F}_j(y)\\ & =\prod _{j\in S}{F}_j(y)+\sum _{x\in S}\frac{y^{a_x-1}}{n^{a_x}\cdot (a_x-1)!}\prod _{j\in S}^{j\ne x}{F}_j(y)\end{array}$$

怎么从导数反推回原式呢？考虑到
$${\left[A_j(y)\cdot e^{yj/n}\right]}^{\prime }=A_j(y)\cdot \frac{j\cdot e^{yj/n}}{n}+A_j^{\prime }(y)\cdot e^{yj/n}$$

那么，当我们知道左式，又知道 $A_j(y)$ 的 $k$ 次项时，我们就可以推出 $A_j^{\prime }(y)$ 的 $k$ 次项，进而知道 $A_j(y)$ 的 $k+1$ 次项。

对应到这个式子上，我们对于每个 $e^{yj/n}$ 单独求系数。此时，我们考虑去计算右式的 $k$ 次项，这当然是 $\mathcal{O}(|S|)$ 的，然后就求出了 $k+1$ 次项。所以总复杂度
$$\mathcal{O}(n^{2}m\cdot 2^n)$$

新能源

这是另一个求解 $\prod (e^{y/n}-\sum \frac{y^b}{n^b\cdot b!})$ 的方法，由 $\mathsf{T}\mathsf{i}\mathsf{w}\text{-}\mathsf{A}\mathsf{i}\mathsf{r}\text{-}\mathsf{O}\mathsf{A}\mathsf{O}$ 提供。
它同时也解决了第一个式子 $\prod \sum \frac{y^b}{n^b\cdot b!}$，可能更高妙。

考虑再枚举一次 $e^{y/n}-{\sum }_{b<a_j}\frac{y^b}{n^b\cdot b!}$ 中，究竟是二者的哪一个，就可以直接得知 $e$ 的指数了。那么我们得到的就是
$$\sum _{T⫅S}(-1{)}^{|T|}{e}^{y(|S|-|T|)/n}\prod _{j\in T}\left(\sum _{b<a_j}\frac{y^b}{n^b\cdot b!}\right)$$

你以为它退化成了 $3^n$ 。注意到 $|T|$ 相等时，$e$ 的次数是固定的，考虑合并它们。定义
$$Z_{S,k}(y)=\sum _{T⫅S}^{|T|=k}\prod _{j\in T}\left(\sum _{b<a_j}\frac{y^b}{n^b\cdot b!}\right)$$

还用上面的法子，求导，不难发现
F T ′ ( y ) = n F T ( y ) − ∑ x ∈ T y a x − 1 n a x − 1 ⋅ ( a x − 1 ) ! F T − { x } F_T'(y)=nF_T(y)-\sum_{x\in T}\frac{y^{a_x-1}}{n^{a_x-1}\cdot(a_x-1)!}F_{T-\{x\}} FT′​(y)=nFT​(y)−x∈T∑​nax​−1⋅(ax​−1)!yax​−1​FT−{x}​

这里 $F_T$ 就是你想的那个意思。进而
$$Z_{S,k}^{\prime }(y)=\sum _{T⫅S}^{|T|=k}\left[n{F}_T(y)-\sum _{x\in T}\cdots \right]$$

（由于后面那一坨太长，我就不写出来啦。直接用上面的式子替换得到的嘛。）
= n Z S , k ( y ) − ∑ x ∈ S y a x − 1 n a x − 1 ⋅ ( a x − 1 ) ! Z S − { x } , k − 1 ( y ) =nZ_{S,k}(y)-\sum_{x\in S}\frac{y^{a_x-1}}{n^{a_x-1}\cdot(a_x-1)!}Z_{S-\{x\},k-1}(y) =nZS,k​(y)−x∈S∑​nax​−1⋅(ax​−1)!yax​−1​ZS−{x},k−1​(y)

同理，这个是 $\mathcal{O}(n)$ 转移，一共 $\mathcal{O}(nm{2}^n)$ 个状态，复杂度仍然是
$$\mathcal{O}(n^{2}m\cdot 2^n)$$

应该也可以过吧？但不是最优解。

思路贰

考场上被我直接否定的做法：“概率的分母在变，还是要不得。”——《论 $\mathsf{O}\mathsf{n}\mathsf{e}\mathsf{I}\mathsf{n}\mathsf{D}\mathsf{a}\mathsf{r}\mathsf{k}$ 为什么弱》

留坑待填

考虑 $k$ 个死掉的随从，按照时间顺序是第 $x_1,x_2,x_3,\dots ,x_k$ 次攻击时被打死的。

那么这个方案的概率一定是
$${\left(\frac{1}{n-k}\right)}^{m-x_k}\prod _{i=0}^{k-1}{\left(\frac{1}{n-i}\right)}^{x_{i+1}}$$

与每次具体是哪个随从受伤无关。有了这一点，我们考虑先勾勒出这个轮廓，然后再确定究竟是哪些随从。

用 $f(i,S)$ 表示，考虑了前 $i$ 次攻击，集合 $S$ 中的随从死掉了，所有这种方案的概率和。如果下一步没有某个随从的死亡，直接乘 $\frac{1}{n-|S|}$ 就行；如果有随从死亡，考虑在前面分配给这个随从 $a_j-1$ 次攻击（因为这一次攻击是肯定打到它的）。所以
f ( i + 1 , S ) = f ( i , S ) n − ∣ S ∣ + ∑ j ∈ S ( i − ∑ x ∈ S x ≠ j a x a j − 1 ) f ( i , S − { j } ) n − ∣ S ∣ + 1 f(i+1,S)=\frac{f(i,S)}{n-|S|}+\sum_{j\in S}{i-\sum_{x\in S}^{x\ne j}a_x\choose a_j-1}\frac{f(i,S-\{j\})}{n-|S|+1} f(i+1,S)=n−∣S∣f(i,S)​+j∈S∑​(aj​−1i−∑x∈Sx​=j​ax​​)n−∣S∣+1f(i,S−{j})​

显然这是 $\mathcal{O}(nm\cdot 2^n)$ 的。这就比 思路壹 少了一个 $n$ 。

折半填坑

但是 $f$ 数组是答案吗？当然不是。因为它还留着很多洞。

再定义 $g(i,S)$ 表示，把 $i$ 次攻击分配给 $S$ 集合中的随从，使得每个人都没有死亡，方案数的和。这个其实就是一个背包。真正的答案是
$$\sum _{S⫅\mathbb{U}}|S|\cdot f(m,S)\cdot g\left(m-\sum _{x\in S}{a}_x,\mathbb{U}-S\right)$$

直接求这个背包，复杂度是 $\mathcal{O}(m^2{2}^n)$，上天了！国集大佬 $\mathsf{C}\mathsf{h}\mathsf{a}\mathsf{r}\mathsf{l}\mathsf{o}\mathsf{t}\mathsf{t}{\mathsf{e}}^{\prime }\mathsf{s}\mathsf{S}\mathsf{p}\mathsf{a}\mathsf{c}\mathsf{e}$ 写了 $\mathtt{N}\mathtt{T}\mathtt{T}$ 优化，还是过不去。

注意到这样一个事实：$g(i,S)$ 其实是 $S$ 集合内的元素对应的指数型生成函数的乘积的 $i$ 次项，即 $g(i,S)=i!\cdot [y^i]{\prod }_{j\in S}({\sum }_{b<a_j}\frac{y^b}{b!})$，然后我们求的是 $G(\mathbb{U}-S)$ 的某一项。

两个多项式相乘，只取一项时，复杂度是 $\mathcal{O}(m)$ 的。于是乎，我们 $\mathtt{m}\mathtt{e}\mathtt{e}\mathtt{t}\text{-}\mathtt{i}\mathtt{n}\text{-}\mathtt{m}\mathtt{i}\mathtt{d}\mathtt{d}\mathtt{l}\mathtt{e}$，就把复杂度优化到了
$$\mathcal{O}(m^2\sqrt{2^n}+m{2}^n)$$

另：此处也可以使用 $\mathsf{T}\mathsf{i}\mathsf{w}\text{-}\mathsf{A}\mathsf{i}\mathsf{r}\text{-}\mathsf{O}\mathsf{A}\mathsf{O}$ 的求导大法。它更通用。

代码

这里就给出 思路贰 的代码了，毕竟好写。

#include <cstdio>
#include <iostream>
#include <cstring>
using namespace std;
inline int readint(){
	int a = 0; char c = getchar(), f = 1;
	for(; c<'0'||c>'9'; c=getchar())
		if(c == '-') f = -f;
	for(; '0'<=c&&c<='9'; c=getchar())
		a = (a<<3)+(a<<1)+(c^48);
	return a*f;
}

const int MaxM = 200;
const int MaxN = 16;
const int Mod = 998244353;
int dp[MaxM+1][1<<MaxN];
int inv[MaxN+1], c[MaxM+1][MaxM+1];
int a[1<<MaxN], cnt[1<<MaxN];

const int HalfN = 8;
int g[2][1<<HalfN][MaxM+1];
int getG(int S,int m){
	int L = S&((1<<HalfN)-1);
	int R = S>>HalfN, ans = 0;
	for(int i=0; i<=m; ++i)
		ans = (ans+1ll*g[0][L][m-i]
			*g[1][R][i]%Mod*c[m][i])%Mod;
	return ans;
}

int main(){
	int n = readint(), m = readint();
	for(int i=(inv[1]=1)+1; i<=n; ++i)
		inv[i] = (0ll+Mod-Mod/i)*inv[Mod%i]%Mod;
	for(int i=0; i<=m; ++i)
	for(int j=c[i][0]=1; j<=i; ++j)
		c[i][j] = (c[i-1][j-1]+c[i-1][j])%Mod;
	for(int i=0; i<n; ++i)
		a[1<<i] = readint();
	for(int S=1; S!=(1<<n); ++S){
		a[S] = a[S^(S&-S)]+a[S&-S];
		cnt[S] = cnt[S^(S&-S)]+1;
	}
	dp[0][0] = 1; // nothing ever happens
	for(int i=1; i<=m; ++i)
	for(int S=0; S<(1<<n); ++S){
		if(a[S] > i) continue; // impossible
		dp[i][S] = 1ll*dp[i-1][S]*inv[n-cnt[S]]%Mod;
		for(int j=0; j<n; ++j) if(S>>j&1)
			dp[i][S] = (dp[i][S]+1ll*dp[i-1][S^(1<<j)]
				*c[i-1-a[S^(1<<j)]][a[1<<j]-1]%Mod
				*inv[n-cnt[S]+1])%Mod;
	}
	g[0][0][0] = g[1][0][0] = 1;
	for(int d=0; d<2; ++d)
	for(int S=1,id; S<(1<<HalfN); ++S){
		for(id=0; !(S>>id&1); ++id);
		g[d][S][0] = 1; // for sure
		for(int i=1; i<=m; ++i)
		for(int j=0; j<a[1<<(id+d*HalfN)]&&j<=i; ++j)
			g[d][S][i] = (g[d][S][i]+1ll*c[i][j]*
				g[d][S^(1<<id)][i-j])%Mod;
	}
	int ans = 0;
	for(int S=1; S<(1<<n); ++S){
		if(a[S] > m) continue; // impossible
		ans = (ans+1ll*cnt[S]*dp[m][S]%Mod
			*getG((1<<n)-1-S,m-a[S]))%Mod;
	}
	printf("%d\n",ans);
	return 0;
}