HDU3359 – Kind of a Blur(高斯消元)

Kind of a Blur

【题意】矩阵变换，从矩阵a到b，a中每个点的曼哈顿距离内的所有sum{a[i][j]}/总个数，为矩阵b[][]中对应a的点值，现给出矩阵b求原矩阵a。

【分析】可以算是第一题高斯消元，学习了。首先，sum{a[i][j]} = 每个sum{a[][]*b[][]}；这样就形成了方程：把a矩阵中的每个点当成未知数x，总共h*w个点，然后把每个点和所有点建一个方程，如果在该点曼哈顿距离内则把方程左边项数值变成1，否则为0，等于b矩阵的值之和。则高斯消元即可求出矩阵a。复杂度(h*w)^3。

【AC CODE】93ms

#include <cstdio>
#include <cstring>
#include <cctype>
#include <cmath>
#include <set>
#include <bitset>
//#include <unordered_set>
#include <queue>
#include <stack>
#include <vector>
#include <string>
#include <algorithm>
using namespace std;
typedef long long LL;
#define rep(i,a,n) for(int i = a; i < n; i++)
#define repe(i,a,n) for(int i = a; i <= n; i++)
#define per(i,n,a) for(int i = n; i >= a; i--)
#define clc(a,b) memset(a,b,sizeof(a))
const int INF = 0x3f3f3f3f, MAXN = 10+5;
double b[MAXN][MAXN];
const double eps = 1e-8;
int n,m;//n-方程个数,m-未知数个数
double a[MAXN*MAXN][MAXN*MAXN],x[MAXN*MAXN];//a[][]方程左边的矩阵,x[]-方程右边的值以及最后的未知数答案

inline bool zero(double num){return fabs(num) < eps;}
int gauss()//返回0无解,2有多个解(但没算出),1有一个解并存在x[]
{
	rep(now,0,m)
	{
		int mx = now;//找出绝对值最大的作为当前行
		rep(i,now+1,n) if(fabs(a[i][now]) > fabs(a[mx][now])) mx = i;
		if(mx != now)
		{
			rep(j,now,m) swap(a[mx][j],a[now][j]);
			swap(x[mx],x[now]);
		}
		if(zero(a[now][now])) return 2;//有多个解
		rep(i,now+1,n)
		{
			double tmp = a[i][now]/a[now][now];
			a[i][now] = 0;
			rep(j,now+1,m) a[i][j] -= tmp*a[now][j];
			x[i] -= tmp*x[now];
		}
	}
	rep(i,0,n)
	{
		bool end = true;
		rep(j,0,m)
		{
			if(!zero(a[i][j]))
			{
				end = false;
				break;
			}
		}
		if(end && !zero(x[i])) return 0;//i行所有a为0，而右边值x不为0，无解
	}
	per(now,m-1,0)//从下往上推,利用已知求未知
	{
		rep(j,now+1,m) x[now] -= a[now][j]*x[j];
		x[now] /= a[now][now];//输出整数解可能需要(int)(x[i]+0.5)
	}
	return 1;
}

int dis(int x1, int y1, int x2, int y2)
{
	return abs(x1-x2)+abs(y1-y2);
}
int main()
{
#ifdef SHY
	freopen("d:\\1.txt", "r", stdin);
#endif
	int w,h,d,count = 0;
	while(~scanf("%d %d %d", &w, &h, &d) && (w+h+d))
	{
		rep(i,0,h)
			rep(j,0,w) scanf("%lf", &b[i][j]);
		clc(a,0);clc(x,0);
		n = m = h*w;
		rep(i,0,h)
		{
			rep(j,0,w)
			{
				rep(k,0,h)
				{
					rep(l,0,w)
					{
						if(dis(i,j,k,l) <= d)
						{
							a[i*w+j][k*w+l] = 1;
							x[i*w+j] += b[i][j];
						}
					}
				}
			}
		}
		int f = gauss();
		if(f != 1) printf("f = %d\n", f);
		if(count++) putchar('\n');
		rep(i,0,h)
		{
			rep(j,0,w) printf("%8.2f", x[i*w+j]);
			putchar('\n');
		}
	}
	return 0;
}