高斯消元模板[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|1typedef 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);//.0const 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 == b1a21*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;}