标签:put i++ name build read 消元 代码 -- pre
高斯消元新玩法
一眼期望\(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;
}
标签:put i++ name build read 消元 代码 -- pre
原文地址:https://www.cnblogs.com/Hs-black/p/12266961.html