高斯消元,直接套模板
#include <CSTDIO>
#include <QUEUE>
using namespace std;
// hdu 3359
/*
高斯消元 求期望
也是数字图像处理的一种应用
*/
const int MAXN = 110;
const double _inf = 1e-9;
double a[MAXN][MAXN], x[MAXN]; // 方程左边的矩阵和等式右边的值, x存放最后结果
int equ, val; // 方程数 未知数个数
inline double mabs(double _X){return _X<0?-_X:_X;}
int gauss()
{
int i,j,k,col,max_r;
for(k=0,col=0;k<equ&&col<val;k++,col++)
{
max_r=k;
for(i=k+1;i<equ;i++)
if(mabs(a[i][col])>mabs(a[max_r][col]))
max_r=i;
if(mabs(a[max_r][col])<_inf) return 0;
if(k!=max_r)
{
for(j=col;j<val;j++)
swap(a[k][j],a[max_r][j]);
swap(x[k],x[max_r]);
}
x[k]/=a[k][col];
for(j=col+1;j<val;j++)a[k][j]/=a[k][col];
a[k][col]=1;
for(i=0;i<equ;i++)
if(i!=k)
{
x[i]-=x[k]*a[i][k];
for(j=col+1;j<val;j++)a[i][j]-=a[k][j]*a[i][col];
a[i][col]=0;
}
}
return 1;
}
int w, h, d;
double data[MAXN][MAXN];
int main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt", "r", stdin);
#endif
int i, j, k, f = 1;
while (scanf("%d%d%d", &w, &h, &d) == 3 && (w+h+d))
{
if (!f)
puts("");
else
f = 0;
for (i = 0; i< h; ++i)
{
for (j = 0; j< w; ++j)
{
scanf("%lf", &data[i][j]);
}
}
equ = val = h*w;
memset(a, 0, sizeof a);
memset(x, 0, sizeof x);
for (i = 0; i<h; ++i)
{
for (j = 0; j< w; ++j)
{
k = i*w+j;
x[k] = 0.0;
int cnt = 0, aa, bb;
for ( aa = -d; aa<= 0; ++aa)
{
for ( bb = -aa-d; bb<= aa+d; ++bb)
{
if ((i+aa)>=0 && (i+aa)<h && (j+bb)>=0 && (j+bb)<w)
{
++cnt;
a[k][(i+aa)*w+(j+bb)] = 1.0;
}
}
}
for ( aa = 1; aa<= d; ++aa)
{
for ( bb = aa-d; bb<= -aa+d; ++bb)
{
if ((i+aa)>=0 && (i+aa)<h && (j+bb)>=0 && (j+bb)<w)
{
++cnt;
a[k][(i+aa)*w+(j+bb)] = 1.0;
}
}
}
x[k] = data[i][j]*cnt;
}
}
gauss();
for (i = 0, k = 0; i< h; ++i)
{
for (j = 0; j< w; ++j, ++k)
{
printf("%8.2lf", x[k]);
}
puts("");
}
}
return 0;
}