【NOI 专题训练】概率期望
- [CTSC2018] 假面 (dp, 求概率)
- CF865C Gotta Go Fast（dp, 二分答案）
- P3750 [六省联考 2017] 分手是祝愿（简化状态，dp）
- [SCOI2008] 奖励关（dp）
- [HNOI2013] 游走（图上随机游走）
- CF1823F Random Walk（树上随机游走， 递推消元）
- [dmy csp模拟赛] 随机游走（树上随机游走，递推消元， 线段树合并）
- 【PKUWC2018】随机游走（状压dp，树上随机游走， 递推消元）
- CF963E Circles of Waiting（特殊图上随机游走，主元法）
- [JSOI2009] 有趣的游戏（AC自动机上随机游走）
- P4007 小 Y 和恐怖的奴隶主（矩阵快速幂，dp）
- [BJOI2018] 治疗之雨（dp, 递推消元）
概率期望降到了9级，感觉很大可能会考。专题训练一下。
[CTSC2018] 假面 (dp, 求概率)
题意：
有 \(n\) 个敌人，第 \(i\) 个敌人的初始生命值为 \(m_i\)。\(q\) 次操作，操作分两种：
1 x p：选中第 \(x\) 个敌人，有 \(p\) 的概率会使他的生命值减1，如果已经减到0就没有影响。
2 k {a_1,a_2,...,a_k}：选中 \(\{a_i\}\) 这些敌人，如果有 \(c\) 个人生命值不为 \(0\)，那么从这 \(c\) 个人中等概率挑一个"锁定"。你需要回答出这 \(k\) 个人各自被"锁定"的概率。
所有操作结束后你需要回答每个敌人剩余生命值的期望大小。
\(1 \leq n \leq 200, 1 \leq m_i \leq 100, 1 \leq q \leq 2 \times 10^5\)。保证二操作的数量不超过 \(1000\)。
分析：
简单题。
设 \(f_{i, j}\) 表示操作到现在，第 \(i\) 个人剩下 \(j\) 个生命值的概率。那么对于 \(1\) 操作可以 \(O(m_i)\) 的进行转移。转移是很平凡的。
对于操作 \(2\)，发现我们要做的是一个背包：设 \(g_{i}\) 表示有 \(i\) 个人生命值不为 \(0\) 的概率。转移很平凡。怎么对每个人回答呢？？发现我们需要的信息是 去掉这个人后其他人做背包得到的结果。然后发现转移是计数背包，并且每个人加入背包的顺序任意。因此可以退背包。这一部分总复杂度 \(O(n^2)\)。
设二操作数量为 \(C\)，总复杂度 \(O(qm_i + Cn^2)\)。
CODE：
// 复杂度 O(Qm_i + Cn^2)
#include<bits/stdc++.h>
using namespace std;
const int N = 210;
const int M = 105;
typedef long long LL;
const LL mod = 998244353;
inline LL Pow(LL x, LL y) {
LL res = 1, k = x;
while(y) {
if(y & 1) res = res * k % mod;
y >>= 1;
k = k * k % mod;
}
return res;
}
inline LL Inv(int v) {return Pow(v, mod - 2);}
int n, m[N], q, idx[N];
LL f[N][M], g[N], inv[N], inv_f[N]; // f[i][j] 表示 i 单位生命值为 j 的概率 g[i] 表示有 i 个没死的概率
int main() {
scanf("%d", &n);
for(int i = 1; i <= n; i ++ ) scanf("%d", &m[i]);
for(int i = 1; i <= n; i ++ ) f[i][m[i]] = 1, inv[i] = Inv(i);
scanf("%d", &q);
while(q -- ) {
int opt; scanf("%d", &opt);
if(!opt) {
int x, u, v; scanf("%d%d%d", &x, &u, &v);
LL p = 1LL * u * Inv(v) % mod;
for(int i = 0; i <= m[x]; i ++ ) {
if(i != 0) f[x][i] = (f[x][i] * (1 - p) % mod + f[x][i + 1] * p % mod + mod) % mod;
else f[x][i] = (f[x][i] + f[x][i + 1] * p % mod) % mod;
}
inv_f[x] = Inv(f[x][0]);
}
else {
int k; scanf("%d", &k);
memset(g, 0, sizeof g); g[0] = 1;
for(int i = 1; i <= k; i ++ ) {
scanf("%d", &idx[i]);
for(int j = i; j >= 0; j -- )
g[j] = ((j ? g[j - 1] * (1 - f[idx[i]][0]) % mod : 0) + g[j] * f[idx[i]][0] % mod + mod) % mod;
}
for(int i = 1; i <= k; i ++ ) {
// 先退背包 idx[i]
if(f[idx[i]][0] == 0) for(int j = 0; j < k; j ++ ) g[j] = g[j + 1];
else {
for(int j = 0; j <= k; j ++ )
g[j] = ((g[j] - (j ? g[j - 1] * (1 - f[idx[i]][0]) % mod : 0) + mod) % mod) * inv_f[idx[i]] % mod;
}
LL res = 0;
for(int j = 0; j < k; j ++ ) res = (res + g[j] * (1 - f[idx[i]][0]) % mod * inv[j + 1] % mod + mod) % mod;
printf("%lld ", res);
if(f[idx[i]][0] == 0) {
for(int j = k; j >= 1; j -- ) g[j] = g[j - 1];
g[0] = 0;
}
else {
for(int j = k; j >= 0; j -- )
g[j] = ((j ? g[j - 1] * (1 - f[idx[i]][0]) % mod : 0) + g[j] * f[idx[i]][0] % mod + mod) % mod;
}
}
puts("");
}
}
for(int i = 1; i <= n; i ++ ) {
LL res = 0;
for(int j = 0; j <= m[i]; j ++ ) res = (res + f[i][j] * j % mod) % mod;
printf("%lld ", res);
}
return 0;
}
CF865C Gotta Go Fast（dp, 二分答案）
题意：
有 \(n\) 个任务。对于第 \(i\) 个任务，你有 \(p_i\) 的概率花费 \(a_i\) 时间和 \(a_i\) 代价完成它，有 \(1 - p_i\) 的概率花费 \(b_i\) 时间和 \(b_i\) 代价完成它。
你需要依次完成 \(n\) 个任务，并且需要满足总用时不超过 \(m\)。
你可以在完成每个任务后自由选择继续进行还是从头开始，如果选择从头开始，那么过程中的用时清 \(0\)，花费的代价不变。
你想要满足总用时不超过 \(m\) 的前提下，最小化期望花费代价。问最优策略下，期望花费代价是多少。
\(1 \leq n \leq 50, 1 \leq a_i < b_i \leq 100, 80 \leq p_i \leq 99,\sum a_i \leq m \leq \sum b_i\)
分析：
能够比较容易设计一个 \(dp\)：
设 \(dp_{i, j}\) 表示已经完成了前 \(i\) 个任务，当前用时为 \(j\)，往后接着做剩下的任务并满足用时限制的最小期望花费代价是多少。
有转移：
我们发现转移是有环的，并且由于取 \(min\) 操作，我们甚至不会高斯消元。
怎么办呢？
设 \(dp_{0, 0} = c\)，先将它看作已知的常数。那么转移变成了：
我们想让 \(dp_{0, 0}\) 能恰好取到 \(c\)，因此构造 \(f_{i, j} = dp_{i, j} - c\)。那么 \(f\) 的转移为：
那么很容易注意到 随着 \(c\) 的增大， \(f_{i, j}\) 单调不增。并且使 \(f_{i, j}\) 取到 \(0\) 的 \(c\) 肯定有且仅有一个。
因此我们想要找到使 \(f_{0, 0}\) 取到 \(0\) 的 \(c\)，只需要二分答案即可。
复杂度 \(O(nm \log V)\)。
CODE：
// 怎么理解呢？感性的话是考虑 dp_{0, 0} = p1() + (1 - p1)()， 其实变成min(dp(0, 0), p1() + (1 - p1)()) 也是一样的，因为这两个值相同。 然后就是如果将 dp(0, 0) 设出来，那么设的值太小取到右边就偏大，太大取到右边就偏小
// 严谨的证明就是考虑 f_{i, j} = dp_{i, j} - c 的转移， 我们想让 f(0, 0) 恰好等于 0， 并且 f(i, j) 随着 c 的增大单调不增。 那么就可以二分这个值
#include<bits/stdc++.h>
using namespace std;
const int N = 55;
const int M = 10010;
const long double eps = 1e-9;
int n, m, a[N], b[N], l[N], r[N];
long double p[N], f[N][M]; // f[i][j] 表示做完 i 用时 j 的最小期望
inline bool check(long double c) {
for(int i = l[n]; i <= r[n]; i ++ ) f[n][i] = (i <= m ? 0 : c);
for(int i = n - 1; i >= 0; i -- ) {
for(int j = l[i]; j <= r[i]; j ++ ) {
if(i == 0 && j == 0) f[i][j] = p[1] * (f[1][a[1]] + a[1]) + (1 - p[1]) * (f[1][b[1]] + b[1]);
else f[i][j] = min(c, p[i + 1] * (f[i + 1][j + a[i + 1]] + a[i + 1]) + (1 - p[i + 1]) * (f[i + 1][j + b[i + 1]] + b[i + 1]));
}
}
return f[0][0] >= c;
}
int main() {
scanf("%d%d", &n, &m);
for(int i = 1; i <= n; i ++ ) {
scanf("%d%d%Lf", &a[i], &b[i], &p[i]);
l[i] = l[i - 1] + a[i]; r[i] = r[i - 1] + b[i];
p[i] /= 100.0;
}
long double l = 0, r = 4e7, mid, res;
while(r - l > eps) {
mid = (l + r) / 2;
if(check(mid)) res = mid, l = mid;
else r = mid;
}
printf("%.9Lf\n", res);
return 0;
}
总结：概率期望 \(dp\) 的转移很容易出现环。如果遇到转移有环并且没有办法高消的情况，可以考虑能否通过二分答案将一个未知数变成定值，此时转移就无环了。着重思考判定条件，可以通过构造 \(f = dp - mid\) 来加以证明理解。
P3750 [六省联考 2017] 分手是祝愿（简化状态，dp）
题意：
给定一个长度为 \(n\) 的 \(0/1\) 序列 \(\{a_i\}\)。定义操作位置 \(x\) 为：将 \(x\) 所有因数下标的 \(a_i\) 异或 \(1\)。
设 \(f(S)\) 表示将 \(S\) 局面变成全 \(0\) 的最少操作次数。同时给定 \(k\)，你需要求出：
对于初始局面 \(\{a_i\}\)，每次随机一个位置操作，直到当前局面 \(S\) 的 \(f(S) \leq k\)，然后按照最优策略操作 \(f(S)\) 次变为全 \(0\)。求这个过程中操作次数的期望，输出期望值乘以 \(n!\) 的结果。
\(1 \leq n \leq 10^5, 0 \leq k \leq n\)。
分析：
自己想出来了。这题还挺妙的。
刚看题发现一点思路都没有，因为操作的过程很抽象。我们先来思考一个状态 \(S\) 的最少操作次数 \(f(S)\) 是多少，因为这个是解决这道题必须要会的。
发现 \(f(S)\) 的求法实际很简单：从达到小依次扫每一个位置，如果当前位置是 \(1\) 就操作它，否则就跳过。
证明也很容易想：对于最高位的 \(1\) 后面的位置，一旦操作了 \(1\) 次就必须要再操作 \(1\) 次变回来，否则不可能变成全 \(0\)。因此操作它们是没有意义的，一定不会操作。那么最高位的 \(1\) 就必须靠操作最高位消掉，然后最高位变小，变成了一个子问题。
根据这个过程，我们进一步发现我们只关注一个状态 \(S\) 消成全 \(0\) 需要操作哪些位置，设这个位置集合为 \(g(S)\)，那么如果操作了一个不属于 \(g(S)\) 的位置 \(x\)，我们要将此时的 \(S'\) 消成全 \(0\) 就变成了要操作 \(g(S) \cup x\) 里面的所有位置。
那么就可以设计状态了：设 \(dp_{S}\) 表示要恰好只将 \(S\) 中的所有 \(1\) 都操作一次的期望次数，那么转移就是分别考虑 \(S\) 中的一个位置被操作和 \(S\) 外的一个位置被操作进行转移。
也由此我们发现我们只关心 \(|S|\)，因此将状态变成 \(f_i\) 表示要将恰好 \(i\) 个位置都操作一次的期望。转移有： \(f_{i} = \frac{i}{n}f_{i - 1} + \frac{n - i}{n}f_{i + 1}\)。
\(k\) 有什么用？相当于初始 \(0 \leq i \leq k, f_{i} = i\)。
不用高斯消元，移项化简递推即可。 最后需要求出 \(\{a_i\}\) 需要操作的位置大小，直接暴力枚举因数即可。这一部分不是难点。
复杂度 \(O(n \sqrt{n})\)。
#include<bits/stdc++.h>
using namespace std;
const int N = 1e5 + 10;
typedef long long LL;
const LL mod = 1e5 + 3;
int n, k, a[N];
LL f[N], v[N], inv[N], fac[N];
inline LL Pow(LL x, LL y) {
LL res = 1, k = x;
while(y) {
if(y & 1) res = res * k % mod;
y >>= 1;
k = k * k % mod;
}
return res;
}
int main() {
scanf("%d%d", &n, &k);
for(int i = 1; i <= n; i ++ ) scanf("%d", &a[i]);
for(int i = 0; i <= k; i ++ ) f[i] = i;
fac[0] = 1; v[n] = 1;
for(int i = 1; i <= n; i ++ ) inv[i] = Pow(i, mod - 2), fac[i] = fac[i - 1] * i % mod;
for(int i = n - 1; i > k; i -- ) v[i] = ((n - i) * inv[i] % mod * v[i + 1] % mod + n * inv[i] % mod) % mod;
for(int i = k + 1; i <= n; i ++ ) f[i] = (f[i - 1] + v[i]) % mod;
int cnt = 0;
for(int i = n; i >= 1; i -- ) {
if(a[i]) {
cnt ++;
for(int j = 1; j * j <= i; j ++ )
if(i % j == 0) {a[j] ^= 1; if(j * j != i) a[i / j] ^= 1;}
}
}
printf("%lld\n", f[cnt] * fac[n] % mod);
return 0;
}
[SCOI2008] 奖励关（dp）
题意：
比较简单就不写了。
分析：
最优策略下的期望值，实际上就是要最大化期望得分。
很自然想到设 \(dp_{i, S}\) 表示还剩下 \(i\) 次，当前有的物品状态为 \(S\) 还能获得的分数的最大期望值。
然后转移物品枚举。复杂度 \(O(nK2^n)\)。代码不放了。
[HNOI2013] 游走（图上随机游走）
题意：
给定一张 \(n\) 个点， \(m\) 条边的无向联通图。你需要给每条边赋一个边权，需要保证所有边权构成 \(1 \sim m\) 的排列。然后从 \(1\) 号点开始随机游走，到 \(n\) 号点结束，每次等概率从当前点的所有连边中选一条经过。
你需要最小化游走过程中经过边权和的期望，问最优赋值策略下的期望值是多少。
\(1 \leq n \leq 500, 1 \leq m \leq \frac{n(n - 1)}{2}\)，保证是简单图。
分析：
首先如果边权固定，那么可以设 \(f_{i}\) 表示从 \(i\) 号点开始游走，到 \(n\) 号点的边权和的期望。然后可以列方程 高斯消元。
边权没有确定，如果我们把边权也设成未知数高斯消元，那么复杂度为 \(O(n^6)\)，无法接受。
注意到上述做法最后得到每个边权前面的系数就是 每条边被经过的期望次数。因此我们只要求出 \(E_i\) 表示第 \(i\) 条边被经过的期望次数即可。
\(n\) 的范围比较小，尝试转化成求点的信息。
设 \(f_x\) 表示 \(x\) 号点期望被经过多少次。 观察可知 \(E_{i} = f_{u_i} \times \frac{1}{deg_{u_i}} + f_{v_i} \times \frac{1}{deg_{v_i}}\)。这个可以感性理解。
问题变成了怎么求 \(f_{x}\)，类别 点 \(\to\) 边 的思路，我们也可以列出 点 \(\to\) 点 的方程。
\(f_{x} = \sum\limits_{u \in E(x)} f_{u} \times \frac{1}{deg_{u}}\)
其中 \(E(x)\) 表示 \(x\) 的邻域。
需要特别注意起点和终点：由于最开始就在 \(1\) 号点，因此他的次数要加 \(1\)。对于终点，由于到它就会立刻停止，因此 \(f_{n} = 0\)。
\(f_{1} = 1 + \sum\limits_{u \in E(1)} f_{u} \times \frac{1}{deg_{u}}\)
\(f_{n} = 0\)
高斯消元即可，复杂度 \(O(n^3)\)。
CODE：
#include<bits/stdc++.h>
#define pb emplace_back
using namespace std;
const int N = 505;
const long double eps = 1e-9;
int n, m, u[N * N], v[N * N], deg[N], odr[N * N];
vector< int > E[N];
long double f[N][N], g[N * N]; // 计算 E[i] 表示从 1 随即游走到 n， i 号点期望被经过几次
inline bool cmp(int x, int y) {return g[x] > g[y];}
int main() {
scanf("%d%d", &n, &m);
for(int i = 1; i <= m; i ++ ) {
scanf("%d%d", &u[i], &v[i]);
E[u[i]].pb(v[i]); E[v[i]].pb(u[i]);
deg[u[i]] ++; deg[v[i]] ++;
}
for(int i = 1; i <= n; i ++ ) {
if(i == n) f[i][n] = 1, f[i][0] = 0;
else {
f[i][i] = 1;
for(auto u : E[i]) f[i][u] = -(1.0 / (long double)(deg[u]));
if(i == 1) f[i][0] += 1.0;
}
}
for(int i = 1; i <= n; i ++ ) {
int idx = i;
for(int j = i; j <= n; j ++ ) {
if(abs(f[j][i]) > eps) idx = j;
}
for(int j = 0; j <= n; j ++ ) swap(f[i][j], f[idx][j]);
long double tmp = f[i][i];
for(int j = 0; j <= n; j ++ ) f[i][j] /= tmp;
for(int j = 1; j <= n; j ++ ) {
if(j == i) continue;
long double tmp = f[j][i];
for(int k = 0; k <= n; k ++ ) f[j][k] -= f[i][k] * tmp;
}
}
for(int i = 1; i <= m; i ++ ) {
g[i] = (1.0 / (long double)(deg[u[i]])) * f[u[i]][0] + (1.0 / (long double)(deg[v[i]])) * f[v[i]][0];
odr[i] = i;
}
sort(odr + 1, odr + m + 1, cmp);
long double ret = 0;
for(int i = 1; i <= m; i ++ ) ret += (long double)(1.0 * i) * g[odr[i]];
printf("%.3Lf\n", ret);
return 0;
}
CF1823F Random Walk（树上随机游走， 递推消元）
题意：
给定一棵 \(n\) 个点的树，从 \(s\) 开始随机游走，每次等概率到达当前点的一个邻域，走到 \(t\) 停止。求所有点被经过的期望次数。
\(2 \leq n \leq 2 \times 10^5\)
分析：
树上随机游走的套路。
首先容易想到以 \(s\) 为根，那么 \(t\) 子树内除 \(t\) 外的所有点的答案都是 \(0\)。
设 \(f_{i}\) 表示随机游走过程中 \(i\) 号点被经过的期望次数。与上一题类似，我们可以列出 \(f_i\) 之间的方程：
但是暴力高斯消元复杂度 \(O(n^3)\)，显然 GG。
但是我们注意到这是一棵树，而不是一张图。因此可以优化。
具体的，我们有结论：
- 自下而上化简 \(f_i\) 的表达式，除了根节点外所有结点 \(x\)，\(f_x\) 都可以化简成 \(k_x \times f_{fa_x} + b_x\) 的形式。其中 \(k_x, b_x\) 是可以在化简过程中求出的常数。
证明：归纳即可。
那么我们就可以自下向上维护 \(k_x, b_x\)，然后到根节点 \(root\) 就可以解出 \(f_{root}\)，最后再自上到下还原每个 \(f_i\) 即可。
时间复杂度线性，但是需要快速幂求逆元。
CODE：
#include<bits/stdc++.h>
#define pb emplace_back
using namespace std;
typedef long long LL;
const int N = 2e5 + 10;
const LL mod = 998244353;
int n, s, t, deg[N];
LL f[N], g[N], inv[N];
vector< int > E[N];
inline LL Pow(LL x, LL y) {
LL res = 1, k = x % mod;
while(y) {
if(y & 1) res = res * k % mod;
y >>= 1;
k = k * k % mod;
}
return res;
}
void dfs(int x, int fa) {
if(x == t || (deg[x] == 1 && x != s)) {
f[x] = inv[deg[fa]]; g[x] = 0;
return ;
}
for(auto v : E[x]) {
if(v == fa) continue;
dfs(v, x);
if(v != t) {
f[x] = (f[x] + f[v] * inv[deg[v]] % mod) % mod;
g[x] = (g[x] + g[v] * inv[deg[v]] % mod) % mod;
}
}
if(x != s) {
g[x] = g[x] * Pow(1 - f[x] + mod, mod - 2) % mod;
f[x] = inv[deg[fa]] * Pow(1 - f[x] + mod, mod - 2) % mod;
}
else f[x] = (g[x] + 1) * Pow(1 - f[x] + mod, mod - 2) % mod;
}
void Dfs(int x, int fa) {
if(x == t) return ;
for(auto v : E[x]) {
if(v == fa) continue;
f[v] = (f[v] * f[x] % mod + g[v]) % mod;
Dfs(v, x);
}
}
int main() {
scanf("%d%d%d", &n, &s, &t);
for(int i = 1; i <= n; i ++ ) inv[i] = Pow(i, mod - 2);
for(int i = 1; i < n; i ++ ) {
int u, v; scanf("%d%d", &u, &v);
E[u].pb(v); E[v].pb(u);
deg[u] ++; deg[v] ++;
}
if(s == t) {
for(int i = 1; i <= n; i ++ )
printf("%d ", (i == s));
}
else {
dfs(s, 0);
Dfs(s, 0);
for(int i = 1; i <= n; i ++ )
printf("%lld ", f[i]);
}
return 0;
}
[dmy csp模拟赛] 随机游走（树上随机游走，递推消元， 线段树合并）
题意：
给定一棵 \(n\) 个点的有根有向树，\(1\) 号点为根， 所有树边方向都是父亲指向儿子。
同时添加了 \(m\) 条非树边，对于每一条非树边 \(u \to v\)，满足：
- \(u\) 不是叶子结点
- \(v\) 是 \(u\) 的祖先
从 \(1\) 号点开始随机游走，每次等概率从当前点的所有出边选择一条经过，走到任意一个叶子停止。问期望步数。
\(2 \leq n \leq 2.5 \times 10^5, 0 \leq m \leq 2.5 \times 10^5\)
分析：
是上一道题的加强。
首先要明白，虽然添加了一些非树边，但是得到的图依然是有很多性质的，不能够按照一般图随机游走处理。
设 \(f_i\) 表示从 \(i\) 出发，游走到一个叶子的期望步数，那么有转移：
\(L\) 表示叶子集合， \(E(i)\) 表示 \(i\) 的所有出边到达的点。
我们注意到添加的非树边都是 返祖边，因此可以将一般的树上随机游走结论进行拓展：
- 自下向上化简 \(f_i\)，所有结点 \(x\) 的 \(f_x\) 都可以化简成 \(f_x = b_x + \sum\limits_{u}( k_{x, u} \times f_{u})\) 的形式。其中 \(u\) 是 \(x\) 的祖先，并且存在一条返祖边 \(y \to u\)， 满足 \(y\) 在 \(x\) 的子树中。
证明：归纳即可。
那么只需要能不断向上维护出 \(b_x, k_{x, u}\)，到根结点就可以解出 \(f_1\)。
发现我们要支持的是 将所有系数乘上某个数，合并同一个未知数的系数。相当于 全局乘，合并两个数组。线段树合并就可以维护这个过程，并且只会在一条非树边的深度较大端点出插入一个系数，空间复杂度就是 \(O(m \log n)\)。
时间复杂度 \(O((n + m) \log n)\)。
并且只要我们将线段树合并的过程倒过来，应该也是可以自上到下确定每个 \(f_i\)。因此实际上我们也能对每个点求出以它为起点期望走多少步。 但是思考了一下好像因为是有向边所以不太能计算每个点期望被经过多少次。
CODE：
#include<bits/stdc++.h>
#define pb push_back
using namespace std;
typedef long long LL;
const LL mod = 998244353;
const int N = 2.5e5 + 10;
int n, m, fa[N], u, v, dep[N];
int rt[N], tot;
LL inv[N * 2];
vector< int > E[N];
inline LL Pow(LL x, LL y) {
LL res = 1LL, k = x;
while(y) {
if(y & 1) res = (res * k) % mod;
y >>= 1;
k = k * k % mod;
}
return res;
}
struct SegmentTree {
int ls, rs; LL mul, num;
#define ls(x) t[x].ls
#define rs(x) t[x].rs
#define mul(x) t[x].mul
#define num(x) t[x].num
}t[N * 2 * 25];
int build() {
tot ++; mul(tot) = 1;
return tot;
}
void spread(int p) {
if(p && mul(p) != 1) {
mul(ls(p)) = mul(ls(p)) * mul(p) % mod;
mul(rs(p)) = mul(rs(p)) * mul(p) % mod;
num(ls(p)) = num(ls(p)) * mul(p) % mod;
num(rs(p)) = num(rs(p)) * mul(p) % mod;
mul(p) = 1;
}
}
void ins(int p, int lp, int rp, int pos, LL c, int T) {
if(lp == rp) {
if(!T) num(p) = c;
else num(p) = (num(p) + c) % mod;
return ;
}
spread(p);
int mid = (lp + rp >> 1);
if(pos <= mid) {
if(!ls(p)) ls(p) = build();
ins(ls(p), lp, mid, pos, c, T);
}
else {
if(!rs(p)) rs(p) = build();
ins(rs(p), mid + 1, rp, pos, c, T);
}
}
LL ask(int p, int lp, int rp, int pos) {
if(lp == rp) return num(p);
spread(p);
int mid = (lp + rp >> 1);
if(pos <= mid) return ask(ls(p), lp, mid, pos);
else return ask(rs(p), mid + 1, rp, pos);
}
int Merge(int p, int q, int lp, int rp) {
if(!p || !q) return (p ^ q);
if(lp == rp) {
num(p) = (num(p) + num(q)) % mod;
return p;
}
spread(p); spread(q);
int mid = (lp + rp >> 1);
ls(p) = Merge(ls(p), ls(q), lp, mid);
rs(p) = Merge(rs(p), rs(q), mid + 1, rp);
return p;
}
void dfs(int x, int fa) {
rt[x] = build();
dep[x] = dep[fa] + 1;
int sz = E[x].size();
if(!sz) {
ins(rt[x], 0, n, 0, 0, 0);
return ;
}
for(auto v : E[x]) {
if(dep[v] == 0) continue;
if(dep[v] < dep[x]) ins(rt[x], 0, n, v, inv[sz], 1);
}
for(auto v : E[x]) {
if(dep[v] == 0) {
dfs(v, x);
mul(rt[v]) = mul(rt[v]) * inv[sz] % mod; // 区间乘
rt[x] = Merge(rt[x], rt[v], 0, n);
}
}
LL p = ask(rt[x], 0, n, x); // 单点查
p = ((1LL - p) % mod + mod) % mod;
ins(rt[x], 0, n, x, 0, 0);
ins(rt[x], 0, n, 0, 1, 1); // 加入1
mul(rt[x]) = mul(rt[x]) * Pow(p, mod - 2LL) % mod;
}
int main() {
scanf("%d%d", &n, &m);
for(int i = 1; i <= n + m; i ++ ) inv[i] = Pow(1LL * i, mod - 2LL);
for(int i = 2; i <= n; i ++ ) {
scanf("%d", &fa[i]);
E[fa[i]].pb(i);
}
for(int i = 1; i <= m; i ++ ) {
scanf("%d%d", &u, &v); E[u].pb(v);
}
dfs(1, 0);
printf("%lld\n", ask(rt[1], 0, n, 0));
return 0;
}
【PKUWC2018】随机游走（状压dp，树上随机游走， 递推消元）
不太懂为啥要 \(Min-Max\) 容斥， 不太懂为啥是黑。。。
题意：
给定一棵 \(n\) 个结点的树，从 \(x\) 出发，每次等概率从当前点选择一条邻边经过。
\(Q\) 次询问，每次询问给定一个点集 \(S\)，你需要求出从 \(x\) 出发随机游走，直到 \(S\) 中的点都至少经过一次的期望步数。
\(1 \leq n \leq 18, 1 \leq Q \leq 5000\)
分析：
设 \(dp_{p, S}\) 表示从 \(p\) 出发开始随机游走，直到 \(S\) 中的点都至少经过一次的期望步数，同时需要满足 \(p \notin S\)。
那么就可以列出一些方程，但是高斯消元复杂度肯定 gg。
发现如果经过 \(S\) 中的一个点， 那么 \(S\) 就会变小，考虑按照 \(S\) 从小到大的顺序确定 \(dp_{p, S}\)。
假设当前要求 \(dp_{p, S}\)，不妨认为对于 \(S\) 的所有子集 \(S'\)，所有 \(dp_{p, S'}\) 的值已知。那么我们发现整个过程就是一个普通的树上游走！
因此还是可以自下到上将每个点的期望用父亲的期望表示，然后解出根的期望以后自上到下求出每个点的期望。
复杂度 \(O(n \times 2^n)\)。
CODE：
#include<bits/stdc++.h>
#define pb emplace_back
using namespace std;
const int N = 19;
typedef long long LL;
const LL mod = 998244353;
int n, q, root, deg[N];
vector< int > E[N];
LL f[N][1 << 18], k[N], b[N]; // f[i][S] 表示当前在 i，还要经过 S 中的点的期望步数(i \notin S)
inline LL Pow(LL x, LL y) {
LL res = 1, k = x;
while(y) {
if(y & 1) res = res * k % mod;
y >>= 1;
k = k * k % mod;
}
return res;
}
inline LL Inv(LL x) {return Pow(x, mod - 2);}
void dfs(int x, int fa, int S) {
for(auto v : E[x]) {
if(v == fa) continue;
dfs(v, x, S);
}
if((S >> x - 1) & 1) k[x] = 0, b[x] = f[x][S ^ (1 << x - 1)];
else {
k[x] = 0; b[x] = 0; LL p = Inv(deg[x]);
for(auto v : E[x]) {
if(v != fa) k[x] = (k[x] + k[v] * p) % mod, b[x] = (b[x] + b[v] * p) % mod;
}
b[x] = (b[x] + 1) % mod;
if(x != root) b[x] = b[x] * Inv(mod + 1 - k[x]) % mod, k[x] = p * Inv(mod + 1 - k[x]) % mod;
else f[x][S] = b[x] * Inv(mod + 1 - k[x]) % mod;
}
}
void Dfs(int x, int fa, int S) {
for(auto v : E[x]) {
if(v == fa) continue;
f[v][S] = (k[v] * f[x][S] % mod + b[v]) % mod;
Dfs(v, x, S);
}
}
inline void solve(int S) {
if(!S) for(int i = 1; i <= n; i ++ ) f[i][S] = 0;
else {
dfs(root, 0, S);
Dfs(root, 0, S);
}
}
int main() {
scanf("%d%d%d", &n, &q, &root);
for(int i = 1; i < n; i ++ ) {
int u, v; scanf("%d%d", &u, &v);
E[u].pb(v); E[v].pb(u);
deg[u] ++; deg[v] ++;
}
for(int S = 0; S < (1 << n); S ++ ) solve(S);
for(int i = 1; i <= q; i ++ ) {
int k, s = 0; scanf("%d", &k);
for(int j = 1; j <= k; j ++ ) {
int p; scanf("%d", &p);
if(p != root) s |= (1 << p - 1);
}
printf("%lld\n", f[root][s]);
}
return 0;
}
CF963E Circles of Waiting（特殊图上随机游走，主元法）
学到了 主元法。
题意：
直角坐标系上有一个点， 初始在 \((0, 0)\)。这个点会在平面上随机游走，设当前在 \((x, y)\)，那么下一步有四种移动方式：
- 以 \(p_1\) 的概率移动到 \((x - 1, y)\)。
- 以 \(p_2\) 的概率移动到 \((x, y - 1)\)。
- 以 \(p_3\) 的概率移动到 \((x + 1, y)\)。
- 以 \(p_4\) 的概率移动到 \((x, y + 1)\)。
求这个点移动到某个距离原点大于 \(R\) 的点的期望步数。(这里的距离为欧几里德距离，即满足 \(\sqrt{x^2 + y^2} > R\))
\(0 \leq R \le 50\)， \(R\) 为整数。
分析：
首先容易想到的是将园内所有整点的期望设出来， \(f_{x, y}\) 表示从 \((x, y)\) 开始游走走到圆外一点的期望步数。有用的整点数量为 \(O(R^2)\)， 暴力列方程高斯消元复杂度 \(O(R^6)\)。
注意到图是网格图比较特殊，但是尝试后发现我们根本无法仿照 树上随机游走 的思路：每次用某个下一步转移的量来表示当前量，最终合并成一个方程解出一个量。
这时我们引入 主元法：
对于一些特殊图上随机游走问题，或是发现高斯消元列出的方程比较特殊的情况：思路是找到一些 核心变量 设出来，需要满足这些变量能线性表示出所有其余变量，设这些 核心变量 的数量为 \(c\)，那么再找到 \(c\) 个不等价方程进行高斯消元，就可以在 \(O(c^3)\) 内解出所有核心变量。最后对所有其余变量求出表达式的值即可。
可以发现 核心变量 的数量越少，这个算法的复杂度就越低。
对于这道题：
首先将圆心移动到 \((R, R)\) 的位置。然后对于 \(y \in [0, 2R]\) 的每个 \(y\)，找到园内最小的 \(x\)，将所有这些 \((x, y)\) 的期望设出来(一共有 \(2R + 1\) 个)，那么就可以表示园内所有整点的期望了。
为什么呢？
归纳证明：
按照 \(x\) 从小到大依次确定所有点的表达式，假设当前 \(x' \in [1, x]\) 对应的所有 \((x', y)\) 都已经被表示出来了。
那么有：
化简得：
右边都是已知的，所以左边可以被表示出来。那么除了 \(x + 1\) 出现的在 \(x\) 中没有的 \(y\)，其余的 \((x + 1, y)\) 都能被表示出来。这些新出现的 \(y\) 一定在 \(x + 1\) 取到最小横坐标因此已经被设出来了所以不用管。
\(2R + 1\) 个方程在哪呢？
一直往右递推，当转移到圆外的一点时，我们知道它的期望等于 \(0\)，因此每一个 \(y\) 上都能得到一个方程，恰好时 \(2R + 1\) 个方程。
最后总复杂度就是 \(O(R^3)\)。 这个过程可以用 优先队列 维护，每次先转移 \(x\) 小的点。
CODE：
// 主元法
#include<bits/stdc++.h>
using namespace std;
const int N = 110;
typedef long long LL;
const LL mod = 1e9 + 7;
int R, lim;
inline LL Pow(LL x, LL y) {
LL res = 1, k = x;
while(y) {
if(y & 1) res = res * k % mod;
y >>= 1;
k = k * k % mod;
}
return res;
}
struct equation { // 方程
LL f[N];
friend equation operator + (equation a, equation b) {
equation c;
for(int i = 0; i <= lim; i ++ ) c.f[i] = (a.f[i] + b.f[i]) % mod;
return c;
}
friend equation operator * (LL a, equation b) {
equation c;
for(int i = 0; i <= lim; i ++ ) c.f[i] = (a * b.f[i]) % mod;
return c;
}
};
equation A[N][N], B[N]; // A[x][y] 表示 (x, y) 的表示， B[i] 表示第 i 个方程
LL a[4], p[4], inv_p[N];
struct node {
int x, y;
friend bool operator < (node a, node b) {return a.x > b.x;}
};
priority_queue< node > q;
inline void guass() {
for(int i = 1; i <= lim; i ++ ) {
int id = i;
for(int j = i; j <= lim; j ++ ) {
if(B[j].f[i]) {id = j; break;}
}
for(int j = 0; j <= lim; j ++ ) swap(B[id].f[j], B[i].f[j]);
LL inv = Pow(B[i].f[i], mod - 2);
for(int j = 0; j <= lim; j ++ ) B[i].f[j] = (B[i].f[j] * inv) % mod;
for(int j = 1; j <= lim; j ++ ) {
if(j == i) continue;
LL v = B[j].f[i];
for(int k = 0; k <= lim; k ++ ) B[j].f[k] = (B[j].f[k] - B[i].f[k] * v % mod + mod) % mod;
}
}
}
inline int Find(int y) {
int res = R;
for(int i = R; i >= 0; i -- ) {
if((i - R) * (i - R) + (y - R) * (y - R) <= R * R) res = i;
else break;
}
return res;
}
int main() {
cin >> R; lim = 2 * R + 1; LL sum = 0;
for(int i = 1; i <= 4; i ++ ) cin >> a[i], sum += a[i];
for(int i = 1; i <= 4; i ++ ) p[i] = a[i] * Pow(sum, mod - 2) % mod, inv_p[i] = Pow(p[i], mod - 2);
for(int i = 0; i <= 2 * R; i ++ ) { // 对于每个 y 将最靠左的 (x, y) 设出来
int x, y;
A[x = Find(i)][y = i].f[i + 1] = 1;
q.push((node) {x, y});
}
while(!q.empty()) {
node u = q.top(); q.pop();
int x = u.x, y = u.y;
if((x - R) * (x - R) + (y - R) * (y - R) > R * R) {B[y + 1] = A[x][y]; B[y + 1].f[0] = (mod - B[y + 1].f[0]) % mod; continue;}
equation tmp = (A[x][y] + ((mod - p[4]) * A[x][y + 1]));
if(x > 0) tmp = tmp + ((mod - p[1]) * A[x - 1][y]);
if(y > 0) tmp = tmp + ((mod - p[2]) * A[x][y - 1]);
tmp.f[0] = (tmp.f[0] + mod - 1) % mod;
A[x + 1][y] = inv_p[3] * tmp;
q.push((node) {x + 1, y});
}
guass(); // (2R + 1)^3 高斯消元
LL res = 0;
for(int i = 1; i <= lim; i ++ ) res = (res + B[i].f[0] * A[R][R].f[i] % mod) % mod;
res = (res + A[R][R].f[0]) % mod;
cout << res << endl;
return 0;
}
这里还有一道 主元法 的题目：传送门
[JSOI2009] 有趣的游戏（AC自动机上随机游走）
题意：
给定 \(n\) 个仅由前 \(m\) 个大写字母构成的长度为 \(l\) 的字符串 \(s_1,\dots,s_n\)，初始有一个空字符串 \(T\)，每次以 \(p_i\) 的概率在 \(T\) 的末尾加入第 \(i\) 个大写字母，保证 \(\sum\limits_{i = 1}^{m}p_i = 1\)。
你需要对每个 \(s_i\) 求出：
- \(s_i\) 先于其它所有字符串在 \(T\) 里出现的概率。
\(1 \leq n,m,l \leq 10\)。
分析：
比较一眼的题。
多模式串匹配肯定要考虑 \(AC\) 自动机。对 \(s_i\) 建出 \(ACAM\)，每个 \(s_i\) 都对应了 \(trie\) 树上的一个叶子。那么问题就变成了从 \(0\) 号点开始随机游走，对于每个叶子求出它比其它所有叶子都先经过的概率。
枚举叶子 \(p\) 计算答案，设 \(f_i\) 表示从 \(i\) 号点开始随机游走，第一次经过的叶子是枚举的叶子的概率，那么有：
然后暴力高斯消元即可。复杂度 \(O(n \times (nl)^3)\)。
CODE：
#include<bits/stdc++.h>
using namespace std;
const int N = 11;
const double eps = 1e-9;
double p[N], q[N], f[N * N][N * N];
int n, m, l, leaf[N];
char s[N];
struct ACAM {
int tr[N * N][N], fail[N * N], tot;
bool et[N * N];
inline void ins(char *s, int idx) {
int len = strlen(s + 1);
int p = 0;
for(int i = 1; i <= len; i ++ ) {
if(!tr[p][s[i] - 'A']) tr[p][s[i] - 'A'] = ++ tot;
p = tr[p][s[i] - 'A'];
}
et[p] = 1; leaf[idx] = p;
}
inline void build() {
queue< int > q;
for(int i = 0; i < m; i ++ )
if(tr[0][i]) q.push(tr[0][i]);
while(!q.empty()) {
int u = q.front(); q.pop();
for(int i = 0; i < m; i ++ ) {
if(tr[u][i]) {
fail[tr[u][i]] = tr[fail[u]][i];
q.push(tr[u][i]);
}
else tr[u][i] = tr[fail[u]][i];
}
}
}
} AC;
inline double calc(int x) {
int o = AC.tot + 1;
for(int i = 0; i < o; i ++ )
for(int j = 0; j <= o; j ++ ) f[i][j] = 0;
for(int i = 1; i <= n; i ++ ) {
if(i == x) f[leaf[i]][leaf[i]] = 1, f[leaf[i]][o] = 1;
else f[leaf[i]][leaf[i]] = 1, f[leaf[i]][o] = 0;
}
for(int i = 0; i < o; i ++ ) {
if(!AC.et[i]) {
f[i][i] = 1;
for(int j = 0; j < m; j ++ ) {
f[i][AC.tr[i][j]] -= p[j];
}
f[i][o] = 0;
}
}
for(int i = 0; i < o; i ++ ) {
int idx = i;
for(int j = i; j < o; j ++ )
if(fabs(f[j][i]) > eps) {idx = j; break;}
for(int j = 0; j <= o; j ++ ) swap(f[i][j], f[idx][j]);
double tmp = f[i][i];
for(int k = 0; k <= o; k ++ ) f[i][k] /= tmp;
for(int j = 0; j < o; j ++ ) {
if(j == i) continue;
double tmp = f[j][i];
for(int k = 0; k <= o; k ++ ) f[j][k] -= f[i][k] * tmp;
}
}
return f[0][o];
}
int main() {
scanf("%d%d%d", &n, &l, &m);
for(int i = 0; i < m; i ++ ) {
scanf("%lf%lf", &p[i], &q[i]);
p[i] /= q[i];
}
for(int i = 1; i <= n; i ++ ) {
scanf("%s", s + 1);
AC.ins(s, i);
}
AC.build();
for(int i = 1; i <= n; i ++ ) printf("%.2lf\n", calc(i));
return 0;
}
P4007 小 Y 和恐怖的奴隶主（矩阵快速幂，dp）
题意：
有一个血量无限的 \(boss\)，初始时只有一个小兵。
给定小兵的数量上限 \(k\) 和满血生命值 \(m\)，小兵有特殊技能：如果你攻击到一个小兵并且它被攻击后没有 死亡，那么如果当前存活的小兵数量小于 \(k\)，它会召唤一个满血的小兵。(死亡指被攻击后血量 \(\leq 0\))
现在你要攻击 \(n\) 次，每次从存活的所有小兵和 \(boss\) 中随机选择一个攻击，被攻击的对象生命值减少 \(1\)。你需要求出 \(n\) 次攻击后 \(boss\) 血量减少值的期望。
\(T\) 次询问，每次给定 \(n\)，你需要每次求出答案。
\(1 \leq T \leq 10^3, 1 \leq n \leq 10^{18}, 1 \leq m \leq 3, 1 \leq k \leq 8\)。
分析：
注意到 \(m, k\) 很小，对于每个局面，我们只关心每种生命值下还有多少个小兵，因此设状态 \(f_{i, a, b, c}\) 表示当前有 \(a\) 个小兵生命值为 \(1\)， \(b\) 个为 \(2\)， \(c\) 个为 \(3\)，你还能攻击 \(i\) 次的期望减少血量。转移很显然是线性变换，因此可以矩阵加速。
当 \(m = 3, k = 8\) 时状态数最多，打表发现是 \(165\) 个，但是由于还要加上常数 \(1\) 那一项，因此有 \(166\) 个状态。那么我们就得到了一个 \(O(T \times \log_2n \times (165)^3)\)，这显然太大了。
考虑 向量乘矩阵 的复杂度是 \(O(size^2)\)，因此我们可以 \(O(\log_2n \times (165)^3)\) 预处理倍增矩阵，然后每次询问复杂度做到 \(O(T \times \log_2 n \times (165)^2)\)。
还是有点大，需要卡常，考虑将 \(2\) 进制矩阵快速幂改成 \(K\) 进制矩阵快速幂。那么我们需要预处理 \(i \in [0, \log_K n]，T^{1 \times K^i},\dots T^{(K - 1) \times K^i}\)，这个复杂度是 \(O(\log_K n \times (K - 1) \times (165)^3)\)，查询复杂度就变成了 \(O(T \times \log_K n \times (165)^2)\)，测试发现 \(K\) 取 \(4\) 比较优。
然后还有一个卡常是把 \(n\) 离线下来升序求解，每次继承前面的 \(dp\) 数组。然后就过了。
CODE：
#include<bits/stdc++.h>
#define pb emplace_back
using namespace std;
const int N = 166; // 多一个是常数
const int M = 1010;
typedef long long LL;
const int mod = 998244353;
inline int add(int x, int y) {return x + y >= mod ? x + y - mod : x + y;}
inline int mul(int x, int y) {return 1LL * x * y % mod;}
inline int Pow(int x, int y) {
int res = 1, k = x;
while(y) {
if(y & 1) res = mul(res, k);
y >>= 1;
k = mul(k, k);
}
return res;
}
int T, m, K, odr[M], S, lim;
int ans[M], inv[N];
struct matrix {
int mt[N][N];
friend matrix operator * (matrix a, matrix b) {
matrix c; memset(c.mt, 0, sizeof c.mt);
for(int i = 0; i <= lim; i ++ )
for(int j = 0; j <= lim; j ++ )
for(int k = 0; k <= lim; k ++ )
c.mt[i][j] = add(c.mt[i][j], mul(a.mt[i][k], b.mt[k][j]));
return c;
}
friend matrix operator / (matrix a, matrix b) {
matrix c; memset(c.mt, 0, sizeof c.mt);
for(int i = 0; i <= lim; i ++ )
for(int k = 0; k <= lim; k ++ )
c.mt[0][i] = add(c.mt[0][i], mul(a.mt[0][k], b.mt[k][i]));
return c;
}
} Mt[32][4];
LL n[M];
map< vector< int >, int > idx;
map< int, vector< int > > msk;
vector< int > g;
void dfs(int x, int r) {
if(x == m + 1) {
idx[g] = lim; msk[lim] = g; lim ++;
return ;
}
for(int i = 0; i <= r; i ++ ) g.pb(i), dfs(x + 1, r - i), g.pop_back();
}
inline void pre() {
dfs(1, K);
for(int i = 0; i < lim; i ++ ) {
g = msk[i]; int cnt = 0;
for(int j = 0; j < m; j ++ ) cnt += g[j];
cnt ++;
// 攻击 bos
Mt[0][1].mt[i][i] = add(Mt[0][1].mt[i][i], inv[cnt]);
Mt[0][1].mt[lim][i] = add(Mt[0][1].mt[lim][i], inv[cnt]);
for(int j = 0; j < m; j ++ ) { // 枚举攻击什么生命的人
vector< int > tmp = g;
if(tmp[j] > 0) {
if(j == 0) tmp[j] --;
else tmp[j] --, tmp[j - 1] ++, tmp[m - 1] += (cnt - 1 < K);
Mt[0][1].mt[idx[tmp]][i] = add(Mt[0][1].mt[idx[tmp]][i], mul(g[j], inv[cnt]));
}
}
}
for(int i = 0; i < lim; i ++ ) Mt[0][1].mt[i][lim] = 0;
Mt[0][1].mt[lim][lim] = 1;
for(int i = 0; i <= 29; i ++ ) {
for(int j = 2; j <= 3; j ++ ) Mt[i][j] = Mt[i][j - 1] * Mt[i][1];
if(i != 29) Mt[i + 1][1] = Mt[i][3] * Mt[i][1];
}
g.clear(); for(int i = 0; i < m; i ++ ) g.pb(i == m - 1 ? 1 : 0);
S = idx[g];
}
inline bool cmp(int x, int y) {return n[x] < n[y];}
int main() {
ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
cin >> T >> m >> K;
for(int i = 1; i <= K + 1; i ++ ) inv[i] = Pow(i, mod - 2);
pre();
for(int i = 1; i <= T; i ++ ) cin >> n[i], odr[i] = i;
sort(odr + 1, odr + T + 1, cmp);
matrix f; memset(f.mt, 0, sizeof f.mt); f.mt[0][lim] = 1;
LL lst = 0;
for(int i = 1; i <= T; i ++ ) {
int o = odr[i];
LL cnt = n[o] - lst;
int p = 0;
while(cnt) {
if(cnt % 4 != 0) f = f / Mt[p][cnt % 4];
p ++; cnt /= 4;
}
ans[o] = f.mt[0][S], lst = n[o];
}
for(int i = 1; i <= T; i ++ ) cout << ans[i] << endl;
return 0;
}
[BJOI2018] 治疗之雨（dp, 递推消元）
题意：
给定 \(n, m, p, k\)，表示现在有 \(m+ 1\) 个变量，第一个变量的初始值为 \(p\)，上界为 \(n\)，下界为 \(0\)，剩余 \(m\) 个变量上界为 \(\infty\)，下界为 \(-\infty\)。你需要进行若干轮以下的操作，直到第一个变量的值变为 \(0\)。一轮操作定义为：
- 从所有 未达到上界 的变量里 等概率 选择一个，并把它 \(+1\)。如果没有跳过则这一步。
- 先后执行 \(k\) 次：从 未达到下界 的变量里 等概率 选择一个，并把它 \(-1\)。如果没有则跳过这一步。
\(T\) 组询问，每次给出 \(n,m,p, k\)。你需要求出不断进行上述过程直到第一个变量的值变为 \(0\) 的期望操作轮数。
\(1 \leq T \leq 100, 1 \leq p \leq n \leq 1500, 0 \leq m, k \leq 10^9\)。
分析：
设 \(f_i\) 表示第一个变量的值为 \(i\)，变成 \(0\) 的期望操作轮数。注意到进行一轮操作后第一个变量最多增大 \(1\)。因此可以列出方程：
其中 \(p(i, j)\) 表示第一个变量经过一轮操作从 \(i\) 变成 \(j\) 的概率。
发现这好像可以类似树上随机游走的方程进行 递推消元：\(i\) 从大到小递推，每个 \(f_i\) 都可以表示成 \(j < i\) 的状态 \(f_j\) 乘系数加常数的形式。那么递推到 \(f_1\) 就可直接求值，然后再从小到大递推回来即可。
模拟这个过程对状态的系数进行修改更新的时间复杂度为 \(O(n^2)\)，可以接受。
那么现在的问题就在于如何求 \(p(i, j)\)。
将 \(p(i, j)\) 分为三类： \(j = 0\)， \(i = n 且 j \ne 0\)，\(i \ne n 且 j \ne 0\)。
对于第二类，发现一轮操作的第一步加 \(1\) 是没有影响的，因此有式子：
\(p(n , i) = \binom{k}{n -i}(\frac{1}{m + 1})^{n - i}(\frac{m}{m + 1})^{k - (n - i)}\)。
对于第三类，发现我们只关心 \(i - j\) 的值，设 \(i - j = t\)，那么有：
\(p(i, j) = \frac{1}{m + 1} \binom{k}{t + 1} (\frac{1}{m + 1})^{t + 1}(\frac{m}{m + 1})^{k - t - 1} + \frac{m}{m + 1}\binom{k}{t} (\frac{1}{m + 1})^{t}(\frac{m}{m + 1})^{k - t}\)。
只需要预处理 \(\binom{k}{1 \sim n}\) 和 \((\frac{1}{m + 1})^{1 \sim n}\)，\((\frac{m}{m + 1})^{k - 1 \sim n}\) 就可以 \(O(n)\) 求出这两类概率。
第一类比较难，因为好像一旦把第一个变量减成 \(0\)，那么接下来的概率就不能用 \(\frac{1}{m + 1}\) 或 \(\frac{m}{m + 1}\) 表示，而只能用 \(1\) 来表示了。
难道我们要枚举在哪一次把第一个变量减成 \(0\) 吗？这个枚举量是 \(O(k)\) 并且尝试后发现我不会化简。
冷静思考你会发现，这等价于 第一个变量减到 \(0\) 后还可以再减，操作完后它的值小于等于 \(1\) 的概率。
那么就很简单了，设 \(r_i\) 表示进行 \(k\) 次减 \(1\) 操作后，第一个变量减少了 \(i\) 的概率，那么有：
\(r_i = \binom{k}{i}(\frac{1}{m + 1})^{i}(\frac{m}{m + 1})^{k - i}\)。
然后求第一个变量最后小于等于 \(0\)的概率 等于 \(1\) 减去最后大于 \(0\) 的概率，这样枚举量就是 \(O(n)\) 的。
总复杂度 \(O(Tn^2)\)。
CODE：
// 好像还是递推消元
#include<bits/stdc++.h>
using namespace std;
const int N = 1510;
typedef long long LL;
const LL mod = 1e9 + 7;
int T, n, m, p, k;
inline LL Pow(LL x, LL y) {
LL res = 1, k = x;
while(y) {
if(y & 1) res = res * k % mod;
y >>= 1;
k = k * k % mod;
}
return res;
}
inline LL Inv(LL x) {return Pow(x, mod - 2);}
LL fac[N], inv[N], C[N], t[N], q1[N], q2[N]; // t[i] 表示从 i 开始， 每次以 1 / (m + 1) 的概率 -1， 以 m / (m + 1) 的概率 不减， 最后减到 <= 0 的概率
LL g1[N], g2[N], g3[N], h[N][N]; // h[i][j] 表示初始为 i 减到 0 的期望轮数, 对于每个变量 j 的期望。 其中 n + 1 的值是常数
inline void solve() {
cin >> n >> p >> m >> k;
if(k == 0 || (m == 0 && k == 1)) {puts("-1"); return ;} // 无解
else {
memset(C, 0, sizeof C); memset(t, 0, sizeof t); memset(q1, 0, sizeof q1); memset(q2, 0, sizeof q2);
memset(g1, 0, sizeof g1); memset(g2, 0, sizeof g2); memset(g3, 0, sizeof g3);
q1[0] = 1; LL v = Inv(m + 1); for(int i = 1; i <= n; i ++ ) q1[i] = q1[i - 1] * v % mod;
v = v * m % mod; q2[min(n, k)] = Pow(v, k - min(n, k)); for(int i = min(n, k) - 1; i >= 0; i -- ) q2[i] = q2[i + 1] * v % mod;
for(int i = 0; i <= min(n, k); i ++ ) { // C(k, i)
LL tp = 1;
for(int j = k; j > k - i; j -- ) tp = tp * j % mod;
C[i] = tp * inv[i] % mod;
}
for(int i = 0; i <= min(n, k); i ++ ) {
t[i] = 1; LL tp = 0;
for(int j = 0; j < i; j ++ ) { // 枚举减了几次(不合法)
tp = (tp + C[j] * q1[j] % mod * q2[j] % mod) % mod;
}
t[i] = (t[i] - tp + mod) % mod;
}
v = Inv(m + 1);
for(int i = 0; i <= n; i ++ ) { // g1[i] 从 i 开始随机一个过程降到 0 的概率
if(i == 0) g1[i] = 1;
else {
if(i == n) g1[i] = t[i];
else g1[i] = (v * t[i + 1] % mod + (1 - v + mod) * t[i] % mod) % mod;
}
}
for(int i = n; i > 0; i -- ) { // g2[i] 表示从 n 开始随机一个过程变成 i 的概率
g2[i] = C[n - i] * q1[n - i] % mod * q2[n - i] % mod;
}
for(int i = -1; i <= n; i ++ ) { // 降低 i 个， 第一步可以增大
if(i == -1) g3[i + 1] = v * q2[0] % mod; // 增大
else g3[i + 1] = (v * C[i + 1] % mod * q1[i + 1] % mod * q2[i + 1] % mod + (1 - v + mod) * C[i] % mod * q1[i] % mod * q2[i] % mod) % mod;
}
for(int i = n; i >= 1; i -- ) {
for(int j = 0; j <= n + 1; j ++ ) h[i][j] = 0;
h[i][n + 1] = 1;
if(i != n) { // 先加上 i + 1 的贡献
for(int j = 0; j <= i; j ++ ) h[i][j] = (h[i][j] + h[i + 1][j] * g3[0] % mod) % mod;
h[i][n + 1] = (h[i][n + 1] + h[i + 1][n + 1] * g3[0] % mod) % mod;
}
// 然后加上自己的系数
for(int j = 0; j <= i; j ++ ) h[i][j] = (h[i][j] + (j == 0 ? g1[i] : (i == n ? g2[j] : g3[i - j + 1]))) % mod;
LL v = Inv(1 - h[i][i] + mod);
for(int j = 0; j < i; j ++ ) h[i][j] = h[i][j] * v % mod;
h[i][n + 1] = h[i][n + 1] * v % mod;
}
for(int i = 2; i <= n; i ++ ) {
for(int j = 1; j < i; j ++ ) h[i][n + 1] = (h[i][n + 1] + h[j][n + 1] * h[i][j] % mod) % mod;
}
printf("%lld\n", h[p][n + 1]);
}
}
int main() {
ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
cin >> T;
fac[0] = 1; for(int i = 1; i < N; i ++ ) fac[i] = fac[i - 1] * i % mod;
inv[N - 1] = Inv(fac[N - 1]); for(int i = N - 2; i >= 0; i -- ) inv[i] = inv[i + 1] * (i + 1) % mod;
while(T -- ) solve();
return 0;
}
浙公网安备 33010602011771号