UOJ390UNR #3百鸽笼容斥，EGF

Posted 2021-06-25 AThousandMoons

根据老套路，把操作改为可以鞭尸，\\(\\forall k\\)，求最后把 \\(a_k\\) 变为非正数的概率。

对操作序列算 EGF：

\\[F(x)=\\frac{x^{a_k-1}}{(a_k-1)!}\\prod_{i\\ne k}(e^x-\\sum_{j=0}^{a_k-1}\\frac{x^i}{i!}) \\\\ Ans=\\sum_{i\\ge 0}\\frac{i![x^i]F(x)}{n^{i+1}} \\]

爆展开的每一项形如 \\(x^de^{jx}\\)，考虑计算其贡献。

\\[\\begin{aligned}[] [x^de^{jx}]Ans&=n^{-d-1}\\sum_{i\\ge 0}\\frac{j^i(d+i)!}{i!n^i} \\\\ &=\\frac{d!}{n^{d+1}}\\sum_{i\\ge 0}(\\frac jn)^i\\binom{d+i}i \\\\ &=\\frac{d!}{(n-j)^{d+1}} \\end{aligned} \\]

直接暴力算是 6 次方的，不能过，但根据老套路，\\(\\prod_{i\\ne k}\\) 可以"先乘后除"，就是 \\(5\\) 次方的了。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 31, M = N*N, mod = 998244353;
template<typename T>
void rd(T &x){
	int ch = getchar(); x = 0;
	for(;ch < \'0\' || ch > \'9\';ch = getchar());
	for(;ch >= \'0\' && ch <= \'9\';ch = getchar()) x = x * 10 + ch - \'0\';
}
int ksm(int a, int b){
	int r = 1;
	for(;b;b >>= 1, a = (LL)a * a % mod)
		if(b & 1) r = (LL)r * a % mod;
	return r;
}
void qmo(int &x){x += x >> 31 & mod;}
int n, m, a[N], res[N], fac[M], inv[M], f[N][M], g[N][M], pw[N][M];
int main(){
	memset(res, -1, sizeof res); rd(n);
	for(int i = 0;i < n;++ i) rd(a[i]);
	fac[0] = 1;
	for(int i = 1;i < M;++ i) fac[i] = (LL)fac[i-1] * i % mod;
	inv[M-1] = ksm(fac[M-1], mod-2);
	for(int i = M-1;i;-- i) inv[i-1] = (LL)inv[i] * i % mod;
	f[0][0] = 1;
	for(int i = 0;i < n;++ i){
		for(int j = i;~j;-- j)
			for(int k = m;~k;-- k) if(f[j][k])
				for(int l = 0;l < a[i];++ l)
					f[j+1][k+l] = (f[j+1][k+l] + (LL)f[j][k] * inv[l]) % mod;
		m += a[i];
	}
	for(int i = 1, x;i <= n;++ i){
		x = (LL)inv[i] * fac[i-1] % mod; pw[i][0] = 1;
		for(int j = 1;j < M;++ j) pw[i][j] = (LL)pw[i][j-1] * x % mod;
	}
	for(int i = 0;i < n;++ i){
		if(~res[a[i]]){printf("%d ", res[a[i]]); continue;}
		int ans = 0; memcpy(g, f, sizeof g);
		for(int j = 0;j < n;++ j)
			for(int k = 0;k < m;++ k) if(g[j][k]){
				for(int l = 0;l < a[i];++ l)
					qmo(g[j+1][k+l] -= (LL)g[j][k] * inv[l] % mod);
				int tmp = (LL)g[j][k] * fac[k+a[i]-1] % mod * inv[a[i]-1] % mod * pw[j+1][k+a[i]] % mod;
				if(j & 1) qmo(ans -= tmp); else qmo(ans += tmp-mod);
			}
		printf("%d ", res[a[i]] = ans);
	}
}