【HDU】3359

来源:互联网 发布:linux 更改umask 编辑:程序博客网 时间:2024/06/05 06:46

题目链接:这里写链接内容


这里写图片描述


高斯消元模板题。


代码如下:

#include<queue>#include<cmath>#include<stack>#include<cstdio>#include<vector>#include<cstring>#include<iostream>#include<algorithm>using namespace std;typedef long long LL;#define INF 0x3f3f3f3f#define CLR(a,b) memset(a,b,sizeof(a))#define PI acos(-1.0)const double eps = 1e-8;struct Matrix{    int w,h;    double m[105][105];    void init(int th,int tw)    {        h = th;        w = tw;        CLR(m,0);    }    void Guass()    {        for (int i = 0 ; i < h-1 ; i++)        {            //交换行(减小误差)             int t = i;            for (int j = i+1 ; j < h ; j++)            {                if (fabs(m[j][i]) > fabs(m[t][i]))                    t = j;            }            if (t != i)            {                for (int j = 0 ; j < w ; j++)                    swap(m[t][j],m[i][j]);            }            //消元            for (int j = i+1 ; j < h ; j++)            {                if (fabs(m[j][i]) < eps)                    continue;                double mul = m[j][i] / m[i][i];                for (int k = i ; k < w ; k++)                    m[j][k] -= m[i][k] * mul;            }        }        //反向求解        for (int i = h-1 ; i >= 0 ; i--)        {            if (fabs(m[i][i]) < eps)                continue;            for (int j = i+1 ; j < h ; j++)                m[i][w-1] -= m[i][j]*m[j][w-1];            m[i][w-1] /= m[i][i];        }    }};int main(){    int w,h,d;    double a[15][15];    bool isFirst = true;    while (~scanf ("%d %d %d",&w,&h,&d) && h)    {        for (int i = 0 ; i < h ; i++)        {            for (int j = 0 ; j < w ; j++)                scanf ("%lf",&a[i][j]);        }        if (isFirst)            isFirst = false;        else            printf ("\n");        Matrix ans;        ans.init(h*w,h*w+1);        //构造增广矩阵         for (int i = 0 ; i < h ; i++)        {            for (int j = 0 ; j < w ; j++)            {                int cnt = 0;        //计数                 for (int k = i-d ; k <= i+d ; k++)      //枚举行                 {                    if (k < 0 || k >= h)                        continue;                    int last = d - abs(i-k);        //剩余步数                     for (int l = j-last ; l <= j+last ; l++)        //枚举列                     {                        if (l < 0 || l >= w)                            continue;                        cnt++;                        ans.m[i*w+j][k*w+l] = 1;        //对应的是i*w+j个方程                     }                }                ans.m[i*w+j][w*h] = cnt*a[i][j];            }        }        ans.Guass();//      puts("=======");//      for (int i = 0 ; i < w*h ; i++)//      {//          for (int j = 0 ; j < w*h+1 ; j++)//              printf ("%-8.2lf ",ans.m[i][j]);//          printf ("\n");//      }//      puts("=======");        //输出结果        for (int i = 0 ; i < h ; i++)        {            for (int j = 0 ; j < w ; j++)                printf ("%8.2lf",ans.m[i*w+j][ans.w-1]);            printf ("\n");        }    }    return 0;}