高斯消元

高斯消元的步骤

基本思路

和正常的解方程一样，我们对于每一个元和每一个列中只有一个元的系数为1 即形成 1000000=B1 0100000=B2 0010000=B3 0001000=B4 0000100=B5 0000010=B6 0000001=B7 那么就可以很方便的得到解了

基本的步骤

<1>选择一个选有绝对值最大的为主元的主列P <2>把主列的系数改变使第i项为1 <3>于是再我们把其他的列都减去P的k倍（k即为这一列中第i项的系数）

void Gauss(int n){ for(int i=0;i<n;i++)A[i][n]=B[i]; for(int i=0;i<n;i++){ int mr=i;//为了避免精度误差我们选绝对值最大的为主元 for(int j=i;j<n;j++) if(abs(A[mr][i])<abs(A[j][i]))mr=j; for(int j=0;j<=n;j++)swap(A[mr][j],A[i][j]); for(int j=i+1;j<=n;j++)A[i][j]/=A[i][i]; //把主列的系数改变使A[i][i]=1; A[i][i]=1; for(int j=0;j<n;j++){ if(i==j)continue; for(int k=i+1;k<=n;k++) A[j][k]-=A[j][i]*A[i][k]; //把其他列中的含第i项的系数减为0即A[j][i]=0; } 最后我们得到了一个矩阵 /* 1000000=B1 0100000=B2 0010000=B3 0001000=B4 0000100=B5 0000010=B6 0000001=B7 于是x[i]=B[i]=A[i][n] */ } }

应用

在概率dp的过程中对于一些不满足拓扑序的情况（例如图） 我们可以得到很多的方程，于是我们可以使用高斯消元来完成

例题 Where is the canteen

有一个地图，一个人从某个点出发，问走到花园的期望步数为多少

void init(){ for(int i=0;i<n;i++) for(int j=0;j<m;j++){ //我们用A[f(i,j)]=B[f(i,j)]来表示 一个有关f(i,j)的方程 if(str[i][j]=='$'||!mark[i][j]){ A[f(i,j)][f(i,j)]=1;B[f(i,j)]=0; continue; } int has=0; for(int k=0;k<4;k++){ int x=i+rx[k],y=j+ry[k]; if(x>=0&&y>=0&&x<n&&y<m&&mark[x][y]){ A[f(i,j)][f(x,y)]=-1;has++; } } //即 e[f(i,j)]=(e[id1]+e[id2]+e[id3]+e[id4]+4)/4 //即 4*e[f(i,j)]-e[id1]-e[id2]-e[id3]-e[id4]=4 B[f(i,j)]=1.0*has; A[f(i,j)][f(i,j)]=1.0*has; } } #include<bits/stdc++.h> #define f(a,b) a*m+b using namespace std; const int M=20; double A[M<<4][M<<4],B[M<<4]; char str[M][M]; int rx[]={0,1,0,-1}; int ry[]={1,0,-1,0}; int n,m,cnt; bool mark[M][M]; void dfs(int x,int y){ mark[x][y]=1; if(str[x][y]=='$')cnt++; for(int i=0;i<4;i++){ int x1=x+rx[i],y1=y+ry[i]; if(x1>=0&&x1<n&&y1>=0&&y1<m&&str[x1][y1]!='#'&&!mark[x1][y1])dfs(x1,y1); } } void Gauss(int n){ for(int i=0;i<n;i++)A[i][n]=B[i]; for(int i=0;i<n;i++){ int mr=i; for(int j=i;j<n;j++) if(abs(A[mr][i])<abs(A[j][i]))mr=j; for(int j=0;j<=n;j++)swap(A[mr][j],A[i][j]); for(int j=i+1;j<=n;j++)A[i][j]/=A[i][i]; for(int j=0;j<n;j++){ if(i==j)continue; for(int k=i+1;k<=n;k++) A[j][k]-=A[j][i]*A[i][k]; } } } void init(){ for(int i=0;i<n;i++) for(int j=0;j<m;j++){ //我们用A[f(i,j)]=B[f(i,j)]来表示 一个有关f(i,j)的方程 if(str[i][j]=='$'||!mark[i][j]){ A[f(i,j)][f(i,j)]=1;B[f(i,j)]=0; continue; } int has=0; for(int k=0;k<4;k++){ int x=i+rx[k],y=j+ry[k]; if(x>=0&&y>=0&&x<n&&y<m&&mark[x][y]){ A[f(i,j)][f(x,y)]=-1;has++; } } //即 e[f(i,j)]=(e[id1]+e[id2]+e[id3]+e[id4]+4)/4 //即 4*e[f(i,j)]-e[id1]-e[id2]-e[id3]-e[id4]=4 B[f(i,j)]=1.0*has; A[f(i,j)][f(i,j)]=1.0*has; } } void solve(){ int s1,s2; cnt=0; memset(A,0,sizeof(A)); memset(B,0,sizeof(B)); memset(mark,0,sizeof(mark)); for(int i=0;i<n;i++){ scanf("%s",str[i]); for(int j=0;j<m;j++) if(str[i][j]=='@')s1=i,s2=j; } dfs(s1,s2); if(!cnt){puts("-1");return;} init(); Gauss(n*m); printf("%.6f\n",A[f(s1,s2)][n*m]); } int main(){ while(~scanf("%d %d",&n,&m))solve(); }