Where is the canteen高斯消元

http://acm.hdu.edu.cn/showproblem.php?pid=2262

主要是递推迭代有点多

#include<bits/stdc++.h>
#define ls o<<1
#define rs o<<1|1
#define MS(x,y) memset(x,y,sizeof(x))
#define MC(x,y) memcpy(x,y,sizeof(x))
typedef long long LL; 
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;}
const int N=20,M=0,Z=1e9+7,W=13,L=2e6;
const double eps=1e-10;
const int dy[4]={-1,0,0,1},dx[4]={0,-1,1,0};
using namespace std;
int casenum,casei;
int id;
int n,m,i,j,h,t;
char a[N][N];
int b[N][N]; //how many connects
int vis[N][N]; //the i st reached,
int qy[L],qx[L];// queue
double A[255+5][255+5]; //matrix
int sty,stx;// start  point
bool flag; //reach or not 

//~(-1)=0
void inq(int y,int x){//y i
	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;//minus later
		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);//starts with 0,since id starts with 0
	for(int i=1;i<=n;i++){
		for(int j=1;j<=m;j++)if(~vis[i][j]){//not 0
			int u=vis[i][j];//every point
			A[u][u]=1;// coeffi
			if(a[i][j]=='$')continue;//means 0;
			A[u][id]=1;//every equ same
			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;
				}
			}
		}
	}
}
void gauss_jordan(int equ,int val){
	int i,j,r,c,maxr;
	//0,0 1,1 2,2
	for(r=c=0;r<equ&&c<val;c++){//down
		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]);
		}
		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];
			//delete get the value for every x
			//you need to get the val by divide
			A[r][c]=0;//upper row
		}
	}
}

int main(){
//	freopen("in.txt","r",stdin);
	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("%.06lf\n",A[0][id]/A[0][0]);
	}
	return 0;
}

 https://blog.csdn.net/snowy_smile/article/details/49452699

注意几点

1. fabs()<eps一定要continue

2. A[][]表示系数矩阵,从0开始存

3.从第0列开始,每次取当前列最大值所在的行,与当前行进行交换

4.列方程每一行就是所有的参数,表示一种当前列系数为1时的状态,所有状态加起来概率是1

5.当在‘$’点的时候期望是0

6.最后用增广列除掉系数就是解

7.感觉队列开的有点大

8.解的个数就是参数的个数

9.每次用方阵的感觉算

10.这一列系数无限接近0的时候,这里对应的解有可能是0,就是增广那部分也可能很小

发布了800 篇原创文章 · 获赞 19 · 访问量 6万+
展开阅读全文

没有更多推荐了,返回首页

©️2019 CSDN 皮肤主题: 大白 设计师: CSDN官方博客

分享到微信朋友圈

×

扫一扫,手机浏览