高斯消元模板[HDU2262]

#include<stdio.h>
#include<iostream>
#include<string.h>
#include<ctype.h>
#include<math.h>
#include<map>
#include<set>
#include<vector>
#include<queue>
#include<functional>
#include<string>
#include<algorithm>
#include<time.h>
void fre(){freopen("c://test//input.in","r",stdin);freopen("c://test//output.out","w",stdout);}
template <class T> inline void scand(T &x){char c;x=0;while((c=getchar())<'0');while(c>='0'&&c<='9')x=x*10+(c-48),c=getchar();}
#define MS(x,y) memset(x,y,sizeof(x))
#define MC(x,y) memcpy(x,y,sizeof(x))
#define ls o<<1
#define rs o<<1|1
typedef long long LL;
typedef unsigned int UI;
typedef int Int;
template <class T> inline void gmax(T &a,T b){if(b>a)a=b;}
template <class T> inline void gmin(T &a,T b){if(b<a)a=b;}
using namespace std;
const int N=20,M=0,L=2e6,Z=1e9+7,maxint=2147483647,ms31=522133279,ms63=1061109567,ms127=2139062143;
const double eps=1e-10,PI=acos(-1.0);//.0
const int dy[4]={-1,0,0,1},dx[4]={0,-1,1,0};
map<int,int>mop;
int casenum,casei;
int id;
int n,m,i,j,h,t;
char a[N][N];
int b[N][N];
int vis[N][N];
int qy[L],qx[L];
double A[225+5][225+5];
int sty,stx;
bool flag;
/*
本质上对线性方程组的方程有三种变换
1，交换任意两个方程的位置
2，某一方程两边乘不为零的常数
3，把某一方程的倍数加到另一方程上去

对应着初等行变换——
1，交换矩阵任意两行的位置
2，用不为0的数k乘上矩阵某一行的所有元素
3，将矩阵某一行(乘数k)加到另一张上去

高斯-若尔当(Gauss-Jordan)消元法
===================
概念一——"行阶梯形"
行阶梯形矩阵的特点是：
非零行的第一个非零元素的列标号随着行标号的增加而严格增加
===================
概念二——"行最简形"
行最简形矩阵是特殊的行阶梯形矩阵，它的特点是：
非零行的第一个非零元素为1，而该元素所在列的其他元素全为0
行最简形对应着方程组的解
===================
步骤是
(1)写出线性方程组的增广矩阵
(2)将增广矩阵用初等行变换化成行阶梯形矩阵（等价于消元过程）
(3)判断线性方程组是否有解
    如果行阶梯形矩阵的最后一个非零行代表矛盾方程0=d≠0，则方程组无解
    否则线性方程组有解，并进行下一步
(4)将行阶梯形矩阵用初等行变换化成行最简形矩阵（等价于代入过程）
(5)由行最简形矩阵得线性方程组的解

时间复杂度:
高斯消元代码的时间复杂度是稳定在O(n^3)的
差别一部分体现在矩阵的构建上，另一部分则体现在消元的优化上
（比如利用eps可以使得些没有值的行不需要计算了）

以下模板以HDU2262为例
*/
void inq(int y,int x)
{
    if(y<1||y>n||x<1||x>m||a[y][x]=='#'||~vis[y][x])return;
    vis[y][x]=id++;
    if(a[y][x]=='$'){flag=1;return;}
    qy[t]=y;
    qx[t++]=x;
}
void bfs()
{
    MS(vis,-1);flag=0;
    h=t=0;inq(sty,stx);
    while(h<t)
    {
        int y=qy[h];
        int x=qx[h++];
        b[y][x]=4;
        for(int i=0;i<4;i++)
        {
            if(a[y+dy[i]][x+dx[i]]=='#')b[y][x]--;
            else inq(y+dy[i],x+dx[i]);
        }
    }
}
void build_matrix()
{
    MS(A,0);
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=m;j++)if(~vis[i][j])//有意义的点，从起点可达，$点可以不参与矩阵消元
        {
            int u=vis[i][j];
            A[u][u]=1;//对角线都初始化为1
            if(a[i][j]=='$')continue;
            A[u][id]=1;//否则它的值会传递变成1
            double p=1.0/b[i][j];
            for(int k=0;k<4;k++)//查找所有与其有关系的状态
            {
                int y=i+dy[k];
                int x=j+dx[k];
                if(~vis[y][x])
                {
                    int v=vis[y][x];
                    A[u][v]=-p;
                }
            }
            /*
            如何列式呢？我们设每个点到达终点的期望步数为e[]。
            然后有：e[x]=(e[y]+1)*pxy+(e[y]+1)*pxz
            即得到：e[x]-pxy*e[y]-pxz*e[z]=pxy+pxz
            */
        }
    }
}

/*gauss消元处理的是增广矩阵。
该增广矩阵含有n个方程（n行），m个未知数（m+1列）
a11*x1 + a12*x2 + ... + a1m*xm == b1
a21*x1 + a22*x2 + ... + a2m*xm == b2
...
an1*x1 + an2*x2 + ... + anm*xm == bn*/
//有equ个等式，val个未知数，若算上等式的值列，则共有val+1列
void gauss_jordan(int equ,int val)
{
    int i,j,r,c,maxr;
    for(r=c=0;r<equ&&c<val;c++)//现在想对第c个未知数进行消元
    {
        maxr=r;//为了减少精度误差，常常是找未知数系数绝对值最大的一个
        for(i=r+1;i<equ;i++)if(fabs(A[i][c])>fabs(A[maxr][c]))maxr=i;

        //找不到的话就换下一个未知数
        if(fabs(A[maxr][c])<eps)continue;

        //找到最大行但是不是初始行的话交换上去
        if(maxr!=r)
        {
            for(j=c;j<=val;j++)swap(A[maxr][j],A[r][j]);
            //上面一行也可以直接写成swap(A[maxr],A[r]);
        }

        //然后往下面，针对同一列的所有需要消除的行
        for(i=r+1;i<equ;i++)if(fabs(A[i][c])>eps)
        {
            double k=A[i][c]/A[r][c];
            for(j=c;j<=val;j++)A[i][j]-=k*A[r][j];
        }

        r++;//一行消除完之后再消除下一行      
    }

    //现在是消除得到了行阶梯形，完成下面消除后才能够得到行最简形
    for(c=val-1;c>=1;c--)//枚举每一列（即每一个未知数）
    {
        for(r=c-1;r>=0;r--)if(fabs(A[r][c])>eps)//对上面的每一行做消除
        {
            A[r][id]-=A[r][c]/A[c][c]*A[c][id];
            //其实这里省略了一行——
            A[r][c]=0;
        }
    }

    /*高斯消元还经常要考虑无解或无数解的情况。对于n元线性方程组：
    当增广矩阵的行阶梯形的最后一个非零行代表矛盾方程时，方程组无解；否则方程组有解，且——
    当增广矩阵的行阶梯形有n个非零行时，方程组有唯一解
    当增广矩阵的行阶梯形少于n个非零行时，方程组有无穷多解*/
}
//优化:有时候，如果所求未知量对应于最后一列，在得到行阶梯形后就可以提前输出答案

int main()
{
    while(~scanf("%d%d",&n,&m))
    {
        id=0;

        for(i=1;i<=n;i++)
        {
            scanf("%s",a[i]+1);
            for(j=1;j<=m;j++)if(a[i][j]=='@')
            {
                sty=i;
                stx=j;
            }
        }
        for(i=0;i<=n+1;i++)a[i][0]=a[i][m+1]='#';
        for(j=0;j<=m+1;j++)a[0][j]=a[n+1][j]='#';

        bfs();
        if(flag==0){printf("-1\n");continue;}

        build_matrix();
        gauss_jordan(id,id);
        printf("%.6lf\n",A[0][id]/A[0][0]);
    }
    return 0;
}