BJOI2018 治疗之雨
作者：互联网
前置知识
- 高斯消元
- 快速幂
Description
共有 \(m + 1\) 个数，第一个初始值为 \(p\)，必须在 \([0, n]\) 的区间修改，剩下 \(m\) 个都是无穷。
每轮操作是：
- 不为最大值的数中等概率随机选择一个把它 \(+1\)
- \(k\) 次找一个不为最小值的数把它 \(-1\)
如果不存在就不操作，现在问期望进行多少次操作，第一个数会变成最小值 \(0\)。
Solution
对于这种数学期望的题，比较套路的就是设状态，根据是否有后效性决定高斯消元/DP，这个状态倒着推最好是当前到终点的期望，否则如果正着推，不太好做。
发现 \(m\) 个无穷的数和目标和答案均无关，可以当做一个出气筒（？。所以对于任意一种状态，我们只需要记录第一个数是什么就可以了。
设置状态
设 \(E_i\) 为第一个数为\(i\) 这个状态走到终点的的期望步数。
已知：\(E_0 = 0\)
目标：\(E_{p}\)
状态转移
一些预处理
考虑在不同的血量情况下，一轮操作可能带来加减多少的概率是一定的，不妨先把这个搞出来。
为了简化公式，可以先预处理一个 \(H_i\) 表示在一轮的 \(k\) 次操作中，受到 \(i\) 的伤害（让这个数 \(-1\)）的概率：
\[H_i = C_k^i \times (\dfrac{1}{m + 1}) ^ i \times (\dfrac{m}{m + 1}) ^ {k - i} \]
形象的理解一下，就是从 \(k\) 轮操作中选出 \(i\) 个要提供伤害的贡献，剩下 \(k - i\) 的滚到负无穷。
这预处理都可以 \(O(N ^ 2)\) 时间解决。
预处理完那个东西，我们可以再搞一个 \(W_{i, j}\) 表示从 \(i\) 这个点走到 \(j\) 的概率：
- \(W_{i, i+1} = \dfrac{1}{m + 1}H_0\)，即加一贡献到第一个数，减贡献全部到无穷
- \(W_{i, j} (n > i \ge j) = \dfrac{m}{m + 1}H_{i - j} + \dfrac{1}{m + 1}H_{i - j + 1}\) 就是枚举 \(+1\) 是否作用在第一个数上，所带来的结果
- \(W_{n, j} (n \ge j) = H_{n - j}\)，无论如何都无法 \(+1\)，所以忽略影响即可。
我们并没有考虑到伤害过量的问题，但实际上这些系数都在 \(E_0 = 0\) 前，所以没有贡献。
真正的转移
考虑一个状态 \(i\) 只有可能转移到 \([1, i + 1]\) 转移过来，枚举即可：
\[E_i = 1 + \sum_{j = 1} ^ {i + 1} E_j \times W_{i, j} \]
这个玩意具有后效性，比如 \(1 \Leftrightarrow 2\) 相互转化，所以只能高斯消元，建立 \(n\) 个数 \(E_{1-n}\)的位置数组成的 \(n\) 元 \(1\) 次方程，解出未知数即可。
高斯消元
常规的高斯消元是 \(O(N^3)\) 的，但是此题的矩阵有特殊的性质，可以做到 \(O(N^2)\)。
可以发现此题的矩阵是一个下三角矩阵，即对于一个 \(i\)，与其有关的只有 \([1, i + 1]\)，剩下系数都是 \(0\)。
对于此类特殊的矩阵，我们有 \(O(N ^ 2)\) 的算法。即枚举到一行 \(i\)，的时候，这行的系数不为 \(0\) 的只有 \(i, i + 1\) 两个位置，消掉下面的行每行都是 \(O(1)\) 的（即只用删 \(3\) 个数）
\[\begin{bmatrix} a_{1, 1} & a_{1, 2} & 0 & 0 & a_{1, 5} \\ a_{2, 1} & a_{2, 2} & a_{2, 3} & 0 & a_{2, 5} \\a_{3, 1} & a_{3, 2} & a_{3, 3} & a_{3, 4} & a_{3, 5} \\a_{4, 1} & a_{4, 2} & a_{4, 3} & a_{4, 4} & a_{4, 5} \end{bmatrix} \]
特判
无解情况：
- \(k = 0\)，永远没法减了
- \(m = 0\) 且 \(n > 1\) 且 \(k = 1\) 每次操作永远保持第一个数不变
时间复杂度
\(O(TN^2)\)
Tips
- 此题貌似有点卡常？LOJ不开 \(O_2\) 跑步过去，但是落谷轻松跑过也许这就是玄学吧
- 注意多测清零！
- 求组合数的时候可以把分母乘积算出来再算逆元，这样可以少一个 \(\log\)
#include <iostream>
#include <cstdio>
#include <cstring>
#define rint register int
using namespace std;
typedef long long LL;
const int N = 1505, P = 1e9 + 7;
int n, p, m, k, a[N][N], H[N];
int inline power(int a, int b) {
int res = 1;
while (b) {
if (b & 1) res = (LL) res * a % P;
a = (LL)a * a % P;
b >>= 1;
}
return res;
}
int inline C(int a, int b) {
int res = 1, s = 1;
for (rint i = a, j = 1; j <= b; j++, i--) {
res = (LL)res * i % P;
s = (LL)s * j % P;
}
res = (LL)res * power(s, P - 2) % P;
return res;
}
void inline build() {
int inv = power(m + 1, P - 2);
for (rint i = 0; i <= min(n + 1, k); i++) {
H[i] = (LL)C(k, i) * power(inv, k) % P * power(m, k - i) % P;
}
for (rint i = 1; i <= n; i++) {
a[i][n + 1] = a[i][i] = P - 1;
for (rint j = 1; j <= min(n, i + 1); j++) {
int w;
if (i + 1 == j) w = (LL)H[0] * inv % P;
else if(i == n) w = H[n - j];
else w = ((LL)m * inv % P * H[i - j] % P + (LL)inv * H[i - j + 1] % P) % P;
(a[i][j] += w) %= P;
}
}
}
void inline gauss() {
for (rint i = 1; i <= n; i++) {
int div = power(a[i][i], P - 2);
a[i][i] = 1, a[i][n + 1] = (LL)a[i][n + 1] * div % P;
if (i < n) a[i][i + 1] = (LL)a[i][i + 1] * div % P;
for (rint j = i + 1; j <= n; j++) {
int v = a[j][i]; a[j][i] = 0;
a[j][i + 1] = (a[j][i + 1] - (LL)v * a[i][i + 1] % P + P) % P;
a[j][n + 1] = (a[j][n + 1] - (LL)v * a[i][n + 1] % P + P) % P;
}
}
for (rint i = n - 1; i; i--) {
a[i][n + 1] = (a[i][n + 1] - (LL)a[i][i + 1] * a[i + 1][n + 1] % P + P) % P;
a[i][i + 1] = 0;
}
}
void inline clear() {
for (rint i = 1; i <= n + 1; i++) H[i] = 0;
for (rint i = 1; i <= n; i++)
for (rint j = 1; j <= n + 1; j++) a[i][j] = 0;
}
int main() {
int T; scanf("%d", &T);
while (T--) {
scanf("%d%d%d%d", &n, &p, &m, &k);
if (k == 0 || (n > 1 && m == 0 && k == 1)) puts("-1");
else {
clear();
build();
gauss();
printf("%d\n", a[p][n + 1]);
}
}
return 0;
}
标签：BJOI2018,治疗,第一个,int,res,之雨,dfrac,预处理,高斯消 来源： https://www.cnblogs.com/dmoransky/p/12545669.html