BJOI2018 治疗之雨

落谷、Loj。

前置知识

1. 高斯消元
2. 快速幂

Description

共有 \(m + 1\) 个数，第一个初始值为 \(p\)，必须在 \([0, n]\) 的区间修改，剩下 \(m\) 个都是无穷。

每轮操作是：

1. 不为最大值的数中等概率随机选择一个把它 \(+1\)
2. \(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

1. 此题貌似有点卡常？LOJ不开 \(O_2\) 跑步过去，但是落谷轻松跑过也许这就是玄学吧
2. 注意多测清零！
3. 求组合数的时候可以把分母乘积算出来再算逆元，这样可以少一个 \(\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;
}