高斯消元模板[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

;
int b

;
int vis

;
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;
}