hdu 3359 Kind of a Blur (高斯消元 浮点型)

题目链接

题意：

H * W (W,H <= 10) 的矩阵A的某个元素A[i][j]，从它出发到其他点的曼哈顿距离小于等于D的所有值的和S[i][j]除上可达点的数目，构成了矩阵B。给定矩阵B，求矩阵A。

分析：

将所有矩阵A的元素看成自变量，一共有H*W个变量，每个矩阵B的元素是由这些变量组合而成的，对于固定的B[i][j]，曼哈顿距离在D以内的A[x][y]的系数为1，其它为0，这样就变成了求H*W个变量和H*W个方程的线性方程组，高斯消元求解。这题数据量比较小，所以直接采用浮点数的高斯消元即可，需要注意的是，浮点数消元的时候为了避免精度误差，每次找最大的行，乘上一个小于1的系数进行消元，这样可以把误差降到最小。

本来很快就搞定的，但是脑残了，总是pE，后来又wa，一定要把b[][]变成浮点数，才行，不然会wa 有误差。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#define LL __int64
#define eps 1e-8
const int maxn = 100+10;
using namespace std;
int equ, var;
double a[maxn][maxn], x[maxn];

int Gauss()
{
    int i, j, k, max_r, col;
    double tmp;
    col = 0;

    for(k = 0; k<equ && col<var; k++, col++)
    {
        max_r = k;
        for(i = k+1; i < equ; i++)
            if(fabs(a[i][col])-fabs(a[max_r][col]) > eps)
            max_r = i;

        if(max_r != k)
            for(j = k; j < var+1; j++)
            swap(a[k][j], a[max_r][j]);

        if(fabs(a[k][col]) < eps)
        {
            k--;
            continue;
        }
        for(i = k+1; i < equ; i++)
        {
            if(fabs(a[i][col]) > eps)
            {
                double t = a[i][col]/a[k][col];
                a[i][col] = 0.0;

                for(j = col; j < var+1; j++)
                a[i][j] -= a[k][j]*t;
            }
        }
    }
    for(i = var-1; i >= 0; i--)
    {
        if(fabs(a[i][i]) < eps) continue;
        tmp = a[i][var];
        for(j = i+1; j < var; j++)
        if(a[i][j] != 0)
        tmp -= a[i][j]*x[j];

        //if(tmp%a[i][i] != 0) return -2;
        x[i] = tmp/a[i][i];
    }
    return 0;
}
int mht(int x1, int y1, int x2, int y2)
{
    return abs(x1-x2)+abs(y1-y2);
}

int main()
{
    int n, m, d;
    int f = 0, i, j, k, l;
    double b[maxn][maxn];
    while(cin>>m>>n>>d)
    {
        if(m==0&&n==0&&d==0) break;
        if(f)
        cout<<endl;
        f = 1;
        equ = n*m;
        var = n*m;
        memset(a, 0, sizeof(a));
        memset(x, 0, sizeof(x));
        for(i = 0; i < n; i++)
        for(j = 0; j < m; j++)
        cin>>b[i][j];
        for(i = 0; i < n; i++)
        for(j = 0; j < m; j++)
        for(k = 0; k < n; k++)
        for(l = 0; l < m; l++)
        if(mht(i, j, k, l) <= d)
        {
            a[i*m+j][k*m+l] = 1;
            a[i*m+j][n*m] += b[i][j];
        }
        Gauss();
        for(i = 0; i < n; i++)
        {
            for(j = 0; j < m; j++)
            printf("%8.2lf", x[i*m+j]);
            cout<<endl;
        }
    }
    return 0;
}