CF24D Broken robot(高斯消元)

CF24D Broken robot(高斯消元)

高斯消元新玩法

一眼期望\(dp\), 考虑逆推因为第\(n\)层的期望是确定的(都是\(0\)), \(F[x][y]\)表示从第\(x\)行第\(y\)列开始到第\(n\)层的期望步数

转移方程:

\[F[x][y] = (F[x][y] + F[x][y+1] + F[x][y-1] + F[x+1][y]) / 4 + 1 ~~(y ~!= 1 \&y~!=m)\\F[x][y] = (F[x][y] + F[x][y+1] + F[x+1][y])~ /~ 3 + 1 ~~(y = 1)\\F[x][y] = (F[x][y] + F[x][y-1] + F[x+1][y]) ~/~ 3 + 1 ~~(y = m)\\ \]

应该很好理解⑧, 从左边/右边/下边转移

因为这个\(dp\)有后效性, 考虑用高斯消元解决, 方程如下:

\[3F[x][i] - F[x][i+1] - F[x][i-1] = 4 + F[x+1][i] ~~(y ~!= 1 \&y~!=m)\\2F[x][i] - F[x][i+1] = F[x+1][i] + 3 ~~(y = 1)\\2F[x][i] - F[x][i-1] = F[x+1][i] + 3 ~~(y = m)\\ \]

因为上一行对下一行没有影响, 所以从下至上每行进行高斯消元, 时间复杂度\(\Theta(n^4)\)

仔细观察系数矩阵, 发现很稀疏, 应该有更快的消元

\[\left[ \begin{matrix} 2 & -1 & 0 & 0 & 0\\ -1 & 3 & -1 & 0 & 0\\ 0 & -1 & 3 & -1 & 0\\ 0 & 0 & -1 & 3 & -1\\0 & 0 & 0 & -1 & 2\end{matrix} \right] \]

没错, 每次消元时只是上一行把下一行消掉即可, 只需\(n\)次就可以笑成稀疏的上三角矩阵, 消元时不必整行都消一遍, 因为大部分都是零, 只枚举有数的即可, 最后消成上三角的时候每行只有两个系数, 回代一下就可以了

代码:

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define ll long long
using namespace std;

template <typename T>
void read(T &x) {
    x = 0; bool f = 0;
    char c = getchar();
    for (;!isdigit(c);c=getchar()) if (c=='-') f=1;
    for (;isdigit(c);c=getchar()) x=x*10+(c^48);
    if (f) x=-x;
}

template <typename T>
void write(T x) {
    if (x < 0) putchar('-'), x = -x;
    if (x >= 10) write(x / 10);
    putchar('0' + x % 10);
}

const int N = 1005;
double F[N][N];
int n, m, x, y;

/*

4f[x][i] = (f[x][i] + f[x][i+1] + f[x][i-1] + f[x+1][i]) + 4;
3f[x][i] - f[x][i+1] - f[x][i-1] = 4 + f[x+1][i]

3f[x][i] = (f[x][i] + f[x][i+1] + f[x+1][i]) + 3
2f[x][i] - f[x][i+1] = f[x+1][i] + 3


*/

double M[N][N];
void build(int lin) {
	M[1][1] = M[m][m] = 2, M[1][2] = M[m][m-1] = -1;
	M[1][m+1] = F[lin+1][1] + 3, M[m][m+1] = F[lin+1][m] + 3;
	for (int i = 2;i < m; i++) 
		M[i][i-1] = M[i][i+1] = -1, 
		M[i][i] = 3, M[i][m+1] = F[lin+1][i] + 4;
}

void work(int lin) {
	for (int i = 1;i < m; i++) {
		double k = M[i+1][i] / M[i][i];
		M[i+1][i] = 0, M[i+1][i+1] -= M[i][i+1] * k;
		M[i+1][m+1] -= M[i][m+1] * k;
	}
	
	F[lin][m] = M[m][m+1] / M[m][m];
	for (int i = m - 1;i >= 1; i--) 
		F[lin][i] = (M[i][m+1] - F[lin][i+1] * M[i][i+1]) / M[i][i];
}

int main() {
	read(n), read(m), read(x), read(y);
	
	if (m == 1) return cout << (n - x) * 2 << endl, 0;
	
	for (int i = n - 1;i >= x; i--) 
		build(i), work(i);
	
	printf ("%.10lf", F[x][y]);
	return 0;
}