CF24D Broken robot 题解

一、题目：

codeforces原题

洛谷原题

二、思路：

主要用这道题目来讲一下一种新的算法：Tridiagonal matrix algorithm，简称 TDMA，也叫 Thomas algorithm。

首先说一下这道题目如何入手。

若\(m=1\)，设\(a_i\)为第\(i\)行的答案，则显然\(\forall i<n\)，有\(a_i=\dfrac{1}{2}a_i+\dfrac{1}{2}a_{i+1}+1\)，即\(a_i=a_{i+1}+2\)。又有\(a_n=0\)，所以可得\(a_i=2(N-i)\)。

若\(m>1\)，设\(F(i,j)\)为第\(i\)行第\(j\)列的答案。则显然\(\forall j\in[1,m]\)，\(F(n,j)=0\)。

接下来讨论当\(i<n\)的情况。

1. 若\(j=1\)，则有\(F(i,1)=\dfrac{1}{3}F(i,1)+\dfrac{1}{3}F(i,2)+\dfrac{1}{3}F(i+1,1)+1\)，整理后得到\(2F(i,1)-F(i,2)=3+F(i+1,1)\)。
2. 若\(j=m\)，则有\(F(i,m)=\dfrac{1}{3}F(i,m)+\dfrac{1}{3}F(i,m-1)+\dfrac{1}{3}F(i+1,m)+1\)，整理后得到\(-F(i,m-1)+2F(i,m)=3+F(i+1,m)\)。
3. 若\(1<j<m\)，则有\(F(i,j)=\dfrac{1}{4}F(i,j)+\dfrac{1}{4}F(i,j-1)+\dfrac{1}{4}F(i,j+1)+\dfrac{1}{4}F(i+1,j)+1\)，整理后得到\(-F(i,j-1)+3F(i,j)-F(i,j+1)=4+F(i+1,j)\)。

对于第\(i\)行的情况，显然这是一个\(m\)元一次方程组，用高斯消元可以解出来，但是复杂度为\(O(n\times m^3)\)，炸上天去了。

于是需要用到 TDMA。这是一种能处理 Tridiagonal matrix 的问题。那么什么是 Tridiagonal matrix 呢？

我们观察一下我们列出的方程组，把它写成一个矩阵，可以发现这个矩阵的每一行至多只有三个数。类似这样的矩阵，我们可以在\(O(m)\)的时间内解出来方程组的解。

贴一个地址：Tridiagonal matrix algorithm - TDMA (Thomas algorithm) -- CFD-Wiki, the free CFD reference (cfd-online.com)

三、代码：

#include <iostream>
#include <cstdio>
#include <cstring>

using namespace std;

inline int read(void) {
    int x = 0, f = 1; char ch = getchar();
    while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
    while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
    return f * x;
}

const int maxn = 1005;

int n, m, x0, y0;
double F[maxn][maxn], a[maxn][maxn], b[maxn];

inline void init(void) {
    a[1][1] = 2; a[1][2] = -1;
    a[m][m - 1] = -1; a[m][m] = 2;
    for (int i = 2; i <= m - 1; ++ i) a[i][i] = 3, a[i][i - 1] = a[i][i + 1] = -1;
}

inline void work(double x[]) {
    b[1] += 3; b[m] += 3;
    for (int i = 2; i <= m - 1; ++ i) b[i] += 4;
    init();
    for (int i = 1; i < m; ++ i) {
        double k = a[i][i];
        for (int j = i - 1; j <= i + 1; ++ j) a[i][j] /= k;
        b[i] /= k;

        k = a[i + 1][i];
        for (int j = i; j <= i + 2; ++ j) a[i + 1][j] -= k * a[i][j];
        b[i + 1] -= k * b[i];
    }

    x[m] = b[m] / a[m][m];
    for (int i = m - 1; i >= 1; -- i) {
        x[i] = b[i] - a[i][i + 1] * x[i + 1];
    }
}

int main() {
    n = read(); m = read();
    x0 = read(); y0 = read();
    if (m == 1) {
        printf("%.10lf\n", 2.0 * (n - x0));
        return 0;
    }
    for (int j = 1; j <= m; ++ j) F[n][j] = 0.0;
    for (int i = n - 1; i >= 1; -- i) {
        for (int j = 1; j <= m; ++ j) b[j] = F[i + 1][j];
        work(F[i]);
    }
    printf("%.10lf\n", F[x0][y0]);
    return 0;
}