UOJ#449. 【集训队作业2018】喂鸽子（期望dp）

题意

有 $n$ 只鸽子，每只鸽子需要 $k$ 粒玉米才能喂饱。问每次随意喂给 $n$ 个鸽子中的一个，期望多久所有鸽子都被喂饱。

对于 $998244353$ 取模。

数据范围

$n \le 50, k \le 1000$

题解

$\mathcal O(n^2k \log k)$

题目问的是最晚喂饱的鸽子，我们用 $\min - \max$ 反演变成对于每个集合问最早被喂饱的鸽子。

不难发现只有集合大小是有用的，我们等价于算：

a n s = \sum_{c = 1}^{n} (- 1)^{c + 1} (\binom{n}{c}) g_{c}

$ans = \sum_{c = 1}^{n} (-1)^{c + 1} {n \choose c} g_c$

我们只需要算 $g_c$ 了，即大小为 $c$ 的集合中最早被喂饱鸽子的期望时间。

我们考虑把期望转成概率，即

g_{c} = \sum_{s \geq 1} P (x \geq s)

$g_c = \sum_{s \ge 1} P(x \ge s)$

我们相当于要算到 $s - 1$ 时刻， $c$ 只鸽子都没有被喂饱的概率。

我们为了算这个，辅助设 $f_{c, s}$ 为 $c$ 只鸽子，喂了 $s$ 只玉米还没有被喂饱的概率。

那么就有

g_{c} = \sum_{i \geq 1} \sum_{s = 1}^{i - 1} (\binom{i - 1}{s}) f_{c, s} (\frac{n - c}{n})^{i - 1}

$g_c = \sum_{i \ge 1} \sum_{s = 1}^{i - 1} {i - 1 \choose s} f_{c, s} (\frac{n - c}{n})^{i - 1}$

对于这种式子，我们通常需要交换和式，为了方便令 $\displaystyle p = \frac{n - c}n$ 那么有

g_{c} = \sum_{s = 1}^{c (k - 1)} f_{c, s} \sum_{t \geq 1} (\binom{s + t}{s}) p^{i - 1}

$g_c = \sum_{s = 1}^{c(k - 1)} f_{c, s} \sum_{t \ge 1} {s + t \choose s} p^{i - 1}$

对于后面的式子，我们不难想到一个经典的生成函数形式，即

(\frac{1}{1 - x})^{a} = \sum_{i \geq 0} (\binom{a + i - 1}{a - 1}) x^{i}

$(\frac 1{1 - x})^a = \sum_{i \ge 0} {a + i - 1 \choose a - 1} x^i$

证明，考虑隔板法或者二项式展开。

那么其实就是

(\frac{1}{1 - p})^{s + 1} = (\frac{n}{c})^{1 + c}

$(\frac{1}{1 - p})^{s + 1} = (\frac{n}{c})^{1 + c}\\$

下面我们考虑如何计算 $f_{c, s}$ ，我们枚举第 $c$ 个鸽子喂了几颗玉米，那么就有

f_{c, s} = \sum_{i = 0}^{min (s, k - 1)} (\binom{s}{i}) f_{c - 1, s - i} \frac{1}{n^{i}}

$f_{c, s} = \sum_{i = 0}^{\min(s, k - 1)} {s \choose i} f_{c - 1, s - i} \frac{1}{n^i}$

直接做是 $\mathcal O(n^2k^2)$ 的，用 $NTT$ 优化就可以做到 $\mathcal O(n^2k \log k)$ 啦。

$\mathcal O(n^2k)$

其实有个更高妙的做法。

称一粒玉米是有效玉米，当且仅当它被投喂给了一只没有饱的鸽子。那么有效玉米序列的长度是固定的 $n k$ 。现在考虑枚举所有的有效玉米序列，计算对答案的贡献。下面记 $r_i$ 表示投喂第 $i$ 粒有效玉米前已经有多少鸽子饱了。

那么对于一个玉米序列的贡献其实就是

(\prod_{i = 1}^{n k} P_{r_{i}}) (\sum_{i = 1}^{n k} E_{r_{i}})

$(\prod_{i = 1}^{nk} P_{r_i}) (\sum_{i = 1}^{nk}E_{r_i})$

其中 $\displaystyle P_x = \frac 1 {n - x}, E_x = \frac n{n - x}$ 前面表示的这个序列的概率（注意每个鸽子是不同的），后一项表示相邻两个有效玉米之间需要投递个数的期望。

直接 $dp$ 似乎不太方便。因为无法确定下一粒玉米投喂后是否会是一只鸽子吃饱。注意到贡献只和 $r_i$ 有关，而一只鸽子吃饱前是不会对 $r_i$ 产生影响的。所以可以认为一只鸽子吃饱前其有效玉米都是 **“白玉米” **。我们只在一只鸽子吃饱的时侯把白玉米染色。

这样就可以 $dp$ 了,先强制鸽子吃饱的顺序是 $1$ 到 $n$ ，最后乘 $n!$ 。设 $f_{m, c}$ 表示投喂了 $m$ 粒有效玉米,前 $c$ 只鸽子已经饱了的贡献之和。 $g_{m, c}$ 表示概率之和。

推一下式子，那么有

\sum \prod_{i \leq m} P_{r_{i}} P_{r_{m + 1}} (\sum_{i \leq m} E_{r_{i}} + E_{r_{m} + 1}) = P_{r_{m + 1}} (\sum \prod_{i \leq m} P_{r_{i}} \sum_{i \leq m} E_{r_{i}}) + P_{r_{m + 1}} E_{r_{m + 1}} (\sum \prod_{i \leq m} P_{r_{i}}) = P_{r_{m + 1}} f_{m, c} + P_{r_{m + 1}} E_{r_{m + 1}} g_{m, c}

$\sum \prod_{i \le m} P_{r_i} P_{r_{m + 1}} (\sum_{i \le m} E_{r_i} + E_{r_m + 1})\\ = P_{r_{m + 1}}(\sum \prod_{i \le m} P_{r_i} \sum_{i \le m} E_{r_i}) + P_{r_{m + 1}} E_{r_{m + 1}} (\sum \prod_{i \le m} P_{r_i})\\ = P_{r_{m + 1}}f_{m, c} + P_{r_{m + 1}}E_{r_{m + 1}} g_{m, c}$

显然 $r_{m+1} = c$ 。而新的概率之和只要简单地乘个 $P_{r_{m+1}}$ 就行了。

接下来有两种转移，第一种是加入一粒白玉米，这种直接做。另一种是在 $m + 1$ 处有一只鸽子吃饱了，这种转移需要乘上 $\displaystyle {m−ck \choose k - 1}$ 表示给白玉米染色的方案数。最后有一只鸽子吃饱了 $f_{nk, n} · n!$ 就是答案。

代码

$\mathcal O(n^2k \log k)$

#include <bits/stdc++.h>

#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl

using namespace std;

template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }

inline int read() {
	int x(0), sgn(1); char ch(getchar());
	for (; !isdigit(ch); ch = getchar()) if (ch == '-') sgn = -1;
	for (; isdigit(ch); ch = getchar()) x = (x * 10) + (ch ^ 48);
	return x * sgn;
}

void File() {
#ifdef zjp_shadow
	freopen ("449.in", "r", stdin);
	freopen ("449.out", "w", stdout);
#endif
}

const int N = 55, K = 1010;

namespace Computation {
	const int Mod = 998244353, g = 3;

	inline int fpm(int x, int power) {
		int res = 1;
		for (; power; power >>= 1, x = 1ll * x * x % Mod)
			if (power & 1) res = 1ll * res * x % Mod;
		return res;
	}

	inline void add(int &x, int y) { if ((x += y) >= Mod) x -= Mod; }
	inline void sub(int &x, int y) { if ((x -= y) < 0) x += Mod; }
	inline int mul(int x, int y) { return 1ll * x * y % Mod; }
#define div Div
	inline int div(int x, int y) { return 1ll * x * fpm(y, Mod - 2) % Mod; }

	int fac[N * K], ifac[N * K];
	void Fac_Init(int maxn) {
		fac[0] = ifac[0] = 1;
		For (i, 1, maxn) fac[i] = mul(fac[i - 1], i);
		ifac[maxn] = fpm(fac[maxn], Mod - 2);
		Fordown (i, maxn - 1, 1) ifac[i] = mul(ifac[i + 1], i + 1);
	}
	inline int comb(int n, int m) {
		if (n < 0 || m < 0 || n < m) return 0;
		return mul(mul(fac[n], ifac[n - m]), ifac[m]);
	}
}

namespace Poly {

	using namespace Computation;

	const int Maxn = 1 << 20;

	int powg[Maxn], invpowg[Maxn];

	void NTT_Init() {
		for (int i = 2; i < Maxn; i <<= 1)
			invpowg[i] = fpm(powg[i] = fpm(g, (Mod - 1) / i), Mod - 2);
	}

	int len, rev[Maxn];

	void NTT(int *P, int opt) {
		Rep (i, len) if (i < rev[i]) swap(P[i], P[rev[i]]);
		for (int i = 2, p = 1; i <= len; p = i, i <<= 1) {
			int Wi = opt == 1 ? powg[i] : invpowg[i];
			for (int j = 0; j < len; j += i)
				for (int k = 0, x = 1; k < p; ++ k) {
					int u = P[j + k], v = mul(x, P[j + k + p]);
					P[j + k] = (u + v) % Mod; 
					P[j + k + p] = (u - v + Mod) % Mod; 
					x = mul(x, Wi);
				}
		}
		if (!~opt) {
			int inv = fpm(len, Mod - 2);
			Rep (i, len) P[i] = mul(P[i], inv);
		}
	}

	int A[Maxn], B[Maxn], C[Maxn];
	void Mult(int *a, int *b, int *c, int na, int nb) {
		int nc = na + nb, bit = 0;
		for (len = 1; len <= nc; len <<= 1) ++ bit;
		Rep (i, len) A[i] = B[i] = 0;

		For (i, 0, na) A[i] = a[i];
		For (i, 0, nb) B[i] = b[i];

		Rep (i, len) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
		NTT(A, 1); NTT(B, 1);
		Rep (i, len) C[i] = mul(A[i], B[i]);
		NTT(C, -1);
		For (i, 0, nc) c[i] = C[i];
	}

}

using namespace Computation;

int n, k, f[N][N * K * 2], invn[N * K];

int a[N * K * 2], b[N * K * 2];

int main () {

	File();

	n = read(); k = read();

	Fac_Init(n * k);

	invn[0] = 1; invn[1] = fpm(n, Mod - 2);
	For (i, 2, k) invn[i] = mul(invn[i - 1], invn[1]);

	int ans = 0;

	Poly :: NTT_Init();

	f[0][0] = 1;
	For (c, 1, n) {
		For (s, 0, (c - 1) * (k - 1)) 
			a[s] = mul(f[c - 1][s], ifac[s]);
		For (i, 0, k - 1)
			b[i] = mul(invn[i], ifac[i]);

		Poly :: Mult(a, b, f[c], (c - 1) * (k - 1), k - 1);
		For (s, 0, c * (k - 1))
			f[c][s] = mul(f[c][s], fac[s]);
	}

	For (c, 1, n) {
		int res = 0, base = div(n, c), coef = base;
		For (s, 0, c * (k - 1)) 
			add(res, mul(f[c][s], coef)), coef = mul(coef, base);
		add(ans, mul(comb(n, c), mul(c & 1 ? 1 : Mod - 1, res)));
	}
	printf ("%d\n", ans);

	return 0;

}

$\mathcal O(n^2k)$

#include <bits/stdc++.h>

#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl

using namespace std;

template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }

inline int read() {
	int x(0), sgn(1); char ch(getchar());
	for (; !isdigit(ch); ch = getchar()) if (ch == '-') sgn = -1;
	for (; isdigit(ch); ch = getchar()) x = (x * 10) + (ch ^ 48);
	return x * sgn;
}

void File() {
#ifdef zjp_shadow
	freopen ("449.in", "r", stdin);
	freopen ("449.out", "w", stdout);
#endif
}

const int N = 55, K = 1010;

namespace Computation {
	const int Mod = 998244353;

	inline int fpm(int x, int power) {
		int res = 1;
		for (; power; power >>= 1, x = 1ll * x * x % Mod)
			if (power & 1) res = 1ll * res * x % Mod;
		return res;
	}

	inline void add(int &x, int y) { if ((x += y) >= Mod) x -= Mod; }
	inline void sub(int &x, int y) { if ((x -= y) < 0) x += Mod; }
#define plus Plus
	inline int plus(int x, int y) { return (x += y) >= Mod ? x - Mod : x; }
	inline int mul(int x, int y) { return 1ll * x * y % Mod; }
#define div Div
	inline int div(int x, int y) { return 1ll * x * fpm(y, Mod - 2) % Mod; }

	int fac[N * K], ifac[N * K];
	void Fac_Init(int maxn) {
		fac[0] = ifac[0] = 1;
		For (i, 1, maxn) fac[i] = mul(fac[i - 1], i);
		ifac[maxn] = fpm(fac[maxn], Mod - 2);
		Fordown (i, maxn - 1, 1) ifac[i] = mul(ifac[i + 1], i + 1);
	}
	inline int comb(int n, int m) {
		if (n < 0 || m < 0 || n < m) return 0;
		return mul(mul(fac[n], ifac[n - m]), ifac[m]);
	}
}

using namespace Computation;

int n, k, f[N * K][N], g[N * K][N];

int P[N], E[N], inv[N];

int main () {

	File();

	n = read(); k = read();

	Fac_Init(n * k);

	inv[1] = 1;
	For (i, 2, n)
		inv[i] = mul(inv[Mod % i], (Mod - Mod / i));
	For (i, 0, n)
		P[i] = inv[n - i], E[i] = mul(n, inv[n - i]);

	f[0][0] = 0; g[0][0] = 1;
	Rep (i, n * k) For (j, 0, i / k) if (g[i][j]) {
		int coefg = mul(g[i][j], P[j]),
			coeff = plus(mul(P[j], f[i][j]), mul(mul(P[j], E[j]), g[i][j])),
			coef = comb(i - j * k, k - 1);
		add(f[i + 1][j], coeff);
		add(g[i + 1][j], coefg);
		add(f[i + 1][j + 1], mul(coef, coeff));
		add(g[i + 1][j + 1], mul(coef, coefg));
	}
	printf ("%d\n", mul(f[n * k][n], fac[n]));

	return 0;

}

__EOF__

本文作者：zjp_shadow
本文链接：https://www.cnblogs.com/zjp-shadow/p/10580673.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！

posted @ 2019-03-22 20:03 zjp_shadow 阅读(888) 评论(1) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

阅读排行：
· 地球OL攻略 —— 某应届生求职总结
· 周边上新：园子的第一款马克杯温暖上架
· Open-Sora 2.0 重磅开源！
· 提示词工程——AI应用必不可少的技术
· .NET周刊【3月第1期 2025-03-02】

公告

UOJ#449. 【集训队作业2018】喂鸽子（期望dp）

发表于 2019-03-22 20:03阅读：888评论：1推荐：5

期望概率 dp NTT

关注

跳至底部

昵称： zjp_shadow
园龄： 7年8个月
粉丝： 274
关注： 21

+加关注

随笔分类 (11)

随笔档案 (198)

文章档案 (1)

2018年3月(1)

zjp_shadow

世界是个回音谷，念念不忘，必有回响，你大声喊唱，山谷雷鸣，音传千里，一叠一叠，一浪一浪，彼岸世界都收到了。凡事念念不忘，必有回响。因为它在传递你心间的声音，绵绵不绝，遂相印于心。

UOJ#449. 【集训队作业2018】喂鸽子（期望dp）

题意

数据范围

题解

$\mathcal O(n^2k \log k)$

$\mathcal O(n^2k)$

代码

$\mathcal O(n^2k \log k)$

$\mathcal O(n^2k)$

公告

zjp_shadow

UOJ#449. 【集训队作业2018】喂鸽子（期望dp）

搜索

常用链接

最新随笔

我的标签

积分与排名

随笔分类 (11)

随笔档案 (198)

文章档案 (1)

相册 (41)

友情链接

阅读排行榜

推荐排行榜

最新评论

zjp_shadow

世界是个回音谷，念念不忘，必有回响，你大声喊唱，山谷雷鸣，音传千里，一叠一叠，一浪一浪，彼岸世界都收到了。凡事念念不忘，必有回响。因为它在传递你心间的声音，绵绵不绝，遂相印于心。

UOJ#449. 【集训队作业2018】喂鸽子（期望dp）

题意

数据范围

题解

O(n2klogk)O(n2klog⁡k)\mathcal O(n^2k \log k)

O(n2k)O(n2k)\mathcal O(n^2k)

代码

O(n2klogk)O(n2klog⁡k)\mathcal O(n^2k \log k)

O(n2k)O(n2k)\mathcal O(n^2k)

公告

zjp_shadow

UOJ#449. 【集训队作业2018】喂鸽子（期望dp）

搜索

常用链接

最新随笔

我的标签

积分与排名

随笔分类 (11)

随笔档案 (198)

文章档案 (1)

相册 (41)

友情链接

阅读排行榜

推荐排行榜

最新评论

$\mathcal O(n^2k \log k)$

$\mathcal O(n^2k)$

$\mathcal O(n^2k \log k)$

$\mathcal O(n^2k)$