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 imes (dfrac{1}{m + 1}) ^ i imes (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 imes 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) 个数）

[egin{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
", a[p][n + 1]);
		}
	}
	return 0;
}