UOJ#449. 【集训队作业2018】喂鸽子(期望dp)
题意
有 \(n\) 只鸽子,每只鸽子需要 \(k\) 粒玉米才能喂饱。问每次随意喂给 \(n\) 个鸽子中的一个,期望多久所有鸽子都被喂饱。
对于 \(998244353\) 取模。
数据范围
\(n \le 50, k \le 1000\)
题解
\(\mathcal O(n^2k \log k)\)
题目问的是最晚喂饱的鸽子,我们用 \(\min - \max\) 反演变成对于每个集合问最早被喂饱的鸽子。
不难发现只有集合大小是有用的,我们等价于算:
\[ ans = \sum_{c = 1}^{n} (-1)^{c + 1} {n \choose c} g_c \]
我们只需要算 \(g_c\) 了,即大小为 \(c\) 的集合中最早被喂饱鸽子的期望时间。
我们考虑把期望转成概率,即
\[ g_c = \sum_{s \ge 1} P(x \ge s) \]
我们相当于要算到 \(s - 1\) 时刻,\(c\) 只鸽子都没有被喂饱的概率。
我们为了算这个,辅助设 \(f_{c, s}\) 为 \(c\) 只鸽子,喂了 \(s\) 只玉米还没有被喂饱的概率。
那么就有
\[ 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 \ge 1} {s + t \choose s} p^{i - 1} \]
对于后面的式子,我们不难想到一个经典的生成函数形式,即
\[ (\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}\\ \]
下面我们考虑如何计算 \(f_{c, s}\) ,我们枚举第 \(c\) 个鸽子喂了几颗玉米,那么就有
\[ 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}^{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 \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) << endlusing 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_shadowfreopen ("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 Divinline 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) << endlusing 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_shadowfreopen ("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 Plusinline 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 Divinline 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;}
转载于:https://www.cnblogs.com/zjp-shadow/p/10580673.html
UOJ#449. 【集训队作业2018】喂鸽子(期望dp)相关推荐
- UOJ#449. 【集训队作业2018】喂鸽子
#449. [集训队作业2018]喂鸽子 DP好题 法一:min-max容斥 处理前m个,最快吃饱的鸽子期望的时间 根据期望的定义 考虑每个方案数的概率*期望次数 枚举前m个用了x个,概率都是(1/m ...
- uoj#422. 【集训队作业2018】小Z的礼物
uoj#422. [集训队作业2018]小Z的礼物 题目描述 Solution 所有礼物全部取到的方案数并不好求,因此我们考虑min−maxmin-maxmin−max容斥,转化为第一次取到集合中某一 ...
- 【UOJ#450】【集训队作业2018】复读机(生成函数,单位根反演)
[UOJ#450][集训队作业2018]复读机(生成函数,单位根反演) 题面 UOJ 题解 似乎是\(\mbox{Anson}\)爷的题. \(d=1\)的时候,随便怎么都行,答案就是\(k^n\). ...
- UOJ#418. 【集训队作业2018】三角形
#418. [集训队作业2018]三角形 和三角形没有关系 只要知道儿子放置的顺序,就可以直接模拟了 记录历史最大值 用一个pair(a,b):之后加上a个,期间最大值为增加b个 合并? A1+A2= ...
- UOJ#449. 【集训队作业2018】喂鸽子 min-max容斥,FFT
原文链接www.cnblogs.com/zhouzhendong/p/UOJ449.html 题解 设 f(i) 表示给 i 只鸽子喂食使得至少一只鸽子被喂饱的期望次数,先 min-max容斥 一下. ...
- 【集训队作业2018】喂鸽子
我的计数还是太差了-- 这道题现在知道三种做法. 1. 直接DP 首先显然需要min-max容斥(不知道请百度),不然很难算. 显然对于大小相同的集合答案一样,问题转化为求 \(f_c\) 即 \(c ...
- UOJ#419. 【集训队作业2018】圆形(格林公式)
题面 传送门 题解 首先您得会用格林公式计算圆的面积并 这里只需要动态维护一下圆弧就可以了 时间复杂度\(O(n^2\log n)\) //minamoto #include<bits/stdc ...
- UOJ#450. 【集训队作业2018】复读机 排列组合 生成函数 单位根反演
原文链接https://www.cnblogs.com/zhouzhendong/p/UOJ450.html 题解 首先有一个东西叫做"单位根反演",它在 FFT 的时候用到过: ...
- uoj#420. 【集训队作业2018】矩形(组合数学)
题面 传送门 题解 这辣鸡题目做了咱整整三天--咱果然还是太菜了--好珂怕的推倒啊-- 首先把它变成 \[\left( \sum_{i = 1}^{n} \sum_{j = 1}^{m} F(i, j ...
- uoj#448. 【集训队作业2018】人类的本质(Min_25筛+拉格朗日插值)
题面 传送门 题解 肝了整整一天--膜拜yww和cx巨巨--(虽然它们的题解里我就没看懂几个字) 请备好草稿纸和笔,这种题目就是需要耐心推倒 题目所求是这么一个东西 \[ \begin{aligned ...
最新文章
- .NET笔试题集(一)
- VxWorks中信号量实现任务间通信与同步机制分析
- Git 报错 Push to origin/master was rejected
- 报名|PMCAFF原创专栏作者百人计划
- 新兴机器学习算法:迁移学习
- SystemVerilog声明的位置
- SpringBoot数据访问-------------数据缓存
- IIS安装前已经安装了.NET Framework,安装后如何启用.NETFramework
- SQL Server 用SSMS查看依赖关系有时候不准确,改用代码查
- julia 调用python库_install julia with python
- 关于使用迅雷下载百度云盘文件的教程
- ATmega16 单片机 AVR单片机 自动计价电子秤
- c语言输入一个整数打印出它是奇数还是偶数,1. 编写程序,输入一个整数,打印出它是奇数还是偶数....
- Google File System中文版
- inet_aton和inet_ntoa
- 重力回弹(小球自由落体)
- 多路选择器(multiplexer)简介
- 电源滤波电路的种类、原理及识图技巧
- halcon算子翻译——set_display_font
- FusionCompute8.0.0 实验(3)安装windows虚拟机