@atcoder - AGC038E@ Gachapon

目錄

• @description@
• @solution - 1@
• @accepted code - 1@
• @solution - 2@
• @accepted code - 2@
• @details@

---

@description@

一個隨機數發生器會生成 N 種數。
第 i 種數有參數 Ai，記 SA = ∑Ai，隨機數發生器會有 Ai/SA 的概率生成 i。
每個時刻都會生成一個數，直到某時刻對於所有的 i，第 i 種數被生成的次數 >= Bi 停止。
求停止的期望時刻。

Constraints
1 <= Ai,Bi,N,∑Ai, ∑Bi <= 400。

Input
輸入形式如下：
N
A0 B0
A1 B1
⋮
AN−1 BN−1

Output
輸出期望時刻 mod 998244353。

Sample Input 1
2
1 1
1 1
Sample Output 1
3

@solution - 1@

先對問題進行初步的轉化：
假如在第 t 時刻停止，貢獻為 t。
我們可以將 t 的貢獻看成第 0 時刻沒有停止貢獻 1，第 1 時刻沒有停止貢獻 1，…，第 t-1 時刻沒有停止貢獻 1，總共貢獻 t。
記 p(i) 表示第 i 時刻依然沒有停止的概率，則有停止的期望時刻為 \(\sum_{i=0}p(i)\)。

考慮怎么求 p(i)。
假如給定一個生成的序列 x1, x2, …, xi 表示第 j 時刻生成 xj，則得到該序列的概率應為 A[x1]/SA*A[x2]/SA*…*A[xi]/SA。
一種思路是對於所有長度為 i 的未停止序列計算概率之和，但不好做（但可以做）。
我們逆向思維，用 1 – 已停止序列的概率之和，從而得到 p(i)。

記 d[k] 表示 k 已經出現了 d[k] 次。已停止序列即所有的 1 <= k <= n 都滿足 d[k] >= B[k]。
記 q[i] 表示已停止序列對應的概率，則有：

\[q[i] = \sum_{d[j] \ge B[j]}^{(\sum_{i=1}^{n}d[j])=i}\frac{i!}{\prod_{j=1}^{n}(d[j]!)}*\prod_{j=1}^{n}(\frac{A[j]}{SA})^{d[j]} \]

注意到這個式子里面含有卷積，而且還以階乘為分母。我們考慮將其改寫成指數型生成函數的形式。
記 Q(x) 為 q 的指數型生成函數，則有：

\[Q(x) = \prod_{i=1}^{n}(\sum_{j=B[i]}\frac{(\frac{A[j]}{SA}x)^j}{j!})\\ =\prod_{i=1}^{n}(e^{\frac{A[j]}{SA}x}-\sum_{j=0}^{B[i]-1}\frac{(\frac{A[j]}{SA}x)^j}{j!})\]

設 P(x) 表示 p 的指數型生成函數，得到：

\[P(x) = e^x - Q(x)\\ =e^x - \prod_{i=1}^{n}(e^{\frac{A[j]}{SA}x}-\sum_{j=0}^{B[i]-1}\frac{(\frac{A[j]}{SA}x)^j}{j!})\]

這玩意兒可以暴力 O((∑B)^2*∑A) 展開。
展開后每一項長成 \(c_{i, j}e^{\frac{i}{SA}x}x^{j}\) 的形式，其中 i, j 滿足 0 <= i <= ∑A, 0 <= j <= ∑B。

我們枚舉 i, j，然后算 \(c_{i, j}e^{\frac{i}{SA}x}x^{j}\) 對答案的貢獻。
記 t = i/SA。現暫且不考慮 \(c_{i, j}\)，最后再乘。將 \(e^{tx}x^{j}\) 展開得到：

\[e^{tx}x^j=\sum_{i=0}\frac{t^ix^{i+j}}{i!} \]

因為 \(x^{i+j}\) 對應的系數為 p[i+j]/(i+j)!，所以反過來求 p[i+j] 時需要乘 (i+j)!。因此貢獻為：

\[S(j) = \sum_{i=0}t^i(i+j)^{\underline{j}} \]

怎么求 S(j) 呢？利用下降冪的性質 \((i+1)^{\underline{j}} - i^{\underline{j}} = j*i^{\underline{j-1}}\)，以及推導等比數列求和的過程，得到：

\[S(j) = \sum_{i=0}t^i(i+j)^{\underline{j}}\\ t*S(j) = \sum_{i=1}t^i(i+j-1)^{\underline{j}}\\ (1-t)*S(j) = j! + j*\sum_{i=1}t^i(i+j-1)^{\underline{j-1}} = j*S(j-1) \]

因此得到 \(S(j) = \frac{j}{1-t}S(j-1)\)。
又因 j = 0 時就是等比數列求和，即 \(S(0) = \frac{1}{1-t}\)，所以 \(S(j) = (\frac{1}{1-t})^{j+1}*j!\)。

總時間復雜度 O((∑B)^2*∑A)。

@accepted code - 1@

#include<cstdio>
const int MAXN = 400;
const int MOD = 998244353;
int inv[MAXN + 5];
struct mint{
	int x;
	mint(int _x=0) : x(_x) {}
	friend mint operator + (const mint &a, const mint &b) {return a.x + b.x >= MOD ? a.x + b.x - MOD : a.x + b.x;}
	friend mint operator - (const mint &a, const mint &b) {return a.x - b.x < 0 ? a.x - b.x + MOD : a.x - b.x;}
	friend mint operator * (const mint &a, const mint &b) {return 1LL * a.x * b.x % MOD;}
	friend mint operator / (const mint &a, const int &b) {return a * inv[b];}
};
mint pow_mod(mint b, int p) {
	mint ret = 1;
	while( p ) {
		if( p & 1 ) ret = ret*b;
		b = b*b;
		p >>= 1;
	}
	return ret;
}
void init() {
	inv[1] = 1;
	for(int i=2;i<=MAXN;i++)
		inv[i] = MOD - 1LL*inv[MOD%i]*(MOD/i)%MOD;
}
mint coef[MAXN + 5][MAXN + 5];
int A[MAXN + 5], B[MAXN + 5], N, SA, SB;
void read() {
	scanf("%d", &N);
	for(int i=0;i<N;i++)
		scanf("%d%d", &A[i], &B[i]), SA += A[i], SB += B[i];
}
mint tmp[MAXN + 5][MAXN + 5];
void get() {
	coef[0][0] = 1;
	for(int i=0;i<N;i++) {
		for(int j=0;j<=SA;j++)
			for(int k=0;k<=SB;k++)
				tmp[j][k] = coef[j][k], coef[j][k] = 0;
		mint p = mint(A[i]) / SA, f = MOD - 1;
		for(int l=0;l<B[i];l++,f=f/l*p)
			for(int j=0;j<=SA;j++)
				for(int k=l;k<=SB;k++)
					coef[j][k] = coef[j][k] + f*tmp[j][k-l];
		for(int j=A[i];j<=SA;j++)
			for(int k=0;k<=SB;k++)
				coef[j][k] = coef[j][k] + tmp[j-A[i]][k];
	}
	for(int i=0;i<=SA;i++)
		for(int j=0;j<=SB;j++)
			coef[i][j] = 0 - coef[i][j];
	coef[SA][0] = coef[SA][0] + 1;
}
/*
let t = e^(x/S)
t^S - (t^Ai - (x*Ai/S)^j/j!)
*/
int solve() {
	mint ret = 0, f = 1;
	for(int j=0;j<=SB;j++,f=f*j)
		for(int i=0;i<=SA;i++) {
			mint p = mint(i) / SA;
			p = pow_mod(1 - p, MOD - 2 - j);
			ret = ret + coef[i][j]*p*f;
		}
	return ret.x;
}
int main() {
	init(), read(), get();
	printf("%d\n", solve());
}

@solution - 2@

這道題會很容易令人想起 uoj - 449 那道題。
我們嘗試使用 min-max 容斥能否解決。

對於某一集合 S = {x1, x2, ..., xm}，我們想要求它最早達到 Bi 的限制的期望時刻。
還是轉化一下問題，求第 t 時刻沒有任何元素達到 Bi 的概率 p(t)，將 p(t) 求和即是期望。
因此，這個集合的貢獻如下所示：

\[(-1)^{m+1}*\frac{SA}{\sum A[xi]}*\sum_{0\le di < B[xi]}\frac{(\sum di)!}{\prod(di!)}*\prod(\frac{A[xi]}{\sum A[xi]})^{di} \]

前面的 \((-1)^{m+1}\) 為容斥系數，\(\frac{SA}{\sum A[xi]}\) 為在集合內抽中一個數的期望。

顯然我們不能枚舉每個集合然后算。考慮往集合內加入一個元素 y，貢獻會怎么變：

\[(-1)^{m+2}*\frac{SA}{A[y] + \sum A[xi]}*\sum_{0\le di < B[xi]}\sum_{0\le dy < B[y]}\frac{(dy + \sum di)!}{dy!*\prod(di!)}*\prod(\frac{A[xi]}{A[y] + \sum A[xi]})^{di}*(\frac{A[y]}{A[y] + \sum A[xi]})^{dy} \]

通過觀察式子，假如我們限制以前的集合中 \(\sum A[xi] = p\)，\(\sum di = q\)；並在加入 y 的時候枚舉 dy，則貢獻是不受其他條件影響的。
記 f[i][p][q] 表示已經考慮了 n 個數中的前 i 個數，當前集合 \(\sum A[xi] = p, \sum di = q\)，所有與 p, q 無關的值的乘積之和。
與 p, q 無關的值實際上就是 \((-1)*\frac{A[xi]^{di}}{di!}\)，轉移時乘一下這個即可。

最后再把 f[n][p][q] 拿出來算貢獻即可。
總共轉移不超過 O(∑Bi) 次，時間復雜度為 O((∑Bi)^2*Ai)。

@accepted code - 2@

#include<cstdio>
const int MOD = 998244353;
const int MAXN = 400;
int inv[MAXN + 5];
void init() {
	inv[1] = 1;
	for(int i=2;i<=MAXN;i++)
		inv[i] = MOD - 1LL*inv[MOD%i]*(MOD/i)%MOD;
}
struct modint{
	int x;
	modint(int _x=0):x(_x) {}
	friend modint operator + (modint x, modint y) {return (x.x + y.x) >= MOD ? x.x + y.x - MOD : x.x + y.x;}
	friend modint operator - (modint x, modint y) {return (x.x - y.x) < 0 ? x.x - y.x + MOD : x.x - y.x;}
	friend modint operator * (modint x, modint y) {return 1LL*x.x*y.x%MOD;}
	friend modint operator / (modint x, int k) {return x*inv[k];}
};
modint pow_mod(modint b, int p) {
	modint ret = 1;
	while( p ) {
		if( p & 1 ) ret = ret*b;
		b = b*b;
		p >>= 1;
	}
	return ret;
}
modint f[2][MAXN + 5][MAXN + 5];
int A[MAXN + 5], B[MAXN + 5], SA, SB;
int main() {
	init();
	int n; scanf("%d", &n);
	for(int i=0;i<n;i++)
		scanf("%d%d", &A[i], &B[i]), SA += A[i], SB += B[i];
	f[0][0][0] = MOD - 1;
	for(int i=0;i<n;i++) {
		for(int j=0;j<=SA;j++)
			for(int k=0;k<=SB;k++)
				f[1][j][k] = f[0][j][k];
		modint p = 1;
		for(int l=0;l<B[i];l++,p=p*A[i]/l)
			for(int j=A[i];j<=SA;j++)
				for(int k=l;k<=SB;k++)
					f[0][j][k] = f[0][j][k] - f[1][j-A[i]][k-l]*p;
	}
	modint ans = 0;
	for(int j=1;j<=SA;j++) {
		modint p = 1, q = modint(SA)/j;
		for(int k=0;k<=SB;k++,p=p*k/j)
			ans = ans + f[0][j][k]*p*q;
	}
	printf("%d\n", ans.x);
}

@details@

主要是式子的推導，其實代碼實現本身沒有太大難度。
這個題涉及了好多組合數學以及期望概率的套路。