HDU - 3359 Kind of a Blur - 高斯-约当消元

来源:互联网 发布:windows平板刷机 编辑:程序博客网 时间:2024/05/14 14:03

题目描述

题目有坑,h,w=0,但d可能不等于0就要结束。
h,w输入的顺序也是不合常理,唉,被坑了。

#include<cstdio>#include<cstring>#include<cmath>#include<algorithm>using namespace std;#define MAXH 10#define MAXN 100const double eps=1e-12;int h,w,d,id[MAXH+10][MAXH+10];double a[MAXN+10][MAXN+10],x[MAXN+10];bool free_x[MAXN+10];void Init(){    memset(a,0,sizeof a);    memset(x,0,sizeof x);    memset(id,0,sizeof id);}int Dist(int p,int q,int x,int y){    return abs(y-q)+abs(x-p);}void read(){    double xx;    for(int i=1;i<=h;i++)        for(int j=1;j<=w;j++)            id[i][j]=(i-1)*w+j;    for(int i=1;i<=h;i++)        for(int j=1;j<=w;j++){            scanf("%lf",&xx);            int cnt=0,t=id[i][j];            for(int p=1;p<=h;p++)                for(int q=1;q<=w;q++)                    if(Dist(p,q,i,j)<=d){                        a[t][id[p][q]]=1.0;                        cnt++;                    }            a[t][h*w+1]=1.0*cnt*xx;        }}void GJ_elimination(int equ,int var,int &row,int &col){    int mx;    for(row=col=1;row<=equ&&col<=var;row++,col++){        mx=row;        for(int i=row+1;i<=equ;i++){            if(fabs(a[i][col])>fabs(a[mx][col]))                mx=i;        }        if(mx!=row){            for(int j=1;j<=var+1;j++)                swap(a[mx][j],a[row][j]);        }        if(fabs(a[row][col])<eps){            row--;            continue;        }        for(int i=1;i<=equ;i++){            if(i==row||fabs(a[i][col])<eps)                continue;            for(int j=var+1;j>=col;j--)                a[i][j]-=a[i][col]/a[row][col]*a[row][j];        }    }}int GJ_Judge(int equ,int var,int row,int col){    //0:唯一解    for(int i=row-1;i>=1;i--){        double tmp=0.0;        for(int j=i+1;j<=var;j++)            if(fabs(a[i][j])>eps)                tmp+=a[i][j]*x[j];        x[i]=(a[i][var+1]-tmp)/a[i][i];    }    return 0;}void Gauss_Jordan(int equ,int var){    int row,col;    GJ_elimination(equ,var,row,col);    GJ_Judge(equ,var,row,col);}void print(){    for(int i=1;i<=h*w;i++){        printf("%8.2lf",x[i]);        if(i%w==0)            printf("\n");    }}int main(){    bool flag=false;    while(scanf("%d%d%d",&w,&h,&d)&&w&&h){        if(flag) puts("");        Init();        read();        Gauss_Jordan(h*w,h*w);        print();        flag=true;    }}
0 0