BZOJ4079 [WF2014]Pachinko

Link
首先有一个显然的(O((nm)^3))的高斯消元。
观察发现每次消元最多影响(m)个方程，因此只保留矩阵中的这一部分项（第二维可以用作差法记录），同时只对这些项消元即可。
这样复杂度就降到了(O(nm^3))。

#include<cmath>
#include<cstdio>
#include<cstring>
typedef double db;
const int N=10007,M=23,dx[4]={-1,1,0,0},dy[4]={0,0,-1,1};const db eps=1e-10;
int n,m,cnt,pos[N*M],vis[N*M];
db p[4],a[N*M][M*2],g[N*M];
char s[N][M];
#define f(x,y) a[(x)][(x)-(y)+(M)]
int min(int a,int b){return a<b? a:b;}
int max(int a,int b){return a>b? a:b;}
int main()
{
    while(~scanf("%d%d",&m,&n))
    {
	cnt=0,memset(a,0,sizeof a),memset(g,0,sizeof g),memset(vis,0,sizeof vis);
	for(int i=0;i<4;++i) scanf("%lf",p+i),p[i]/=100;
	for(int i=0;i<n;++i) scanf("%s",s[i]);
	for(int i=0;i<m;++i) if(s[0][i]=='.') ++cnt;
	for(int i=0;i<n;++i)
	    for(int j=0;j<m;++j)
		if(s[i][j]^'X')
		{
		    int P=i*m+j;
		    f(P,P)=1;
		    if(!i) g[P]=1.0/cnt;
		    for(int k=0;k<4;++k)
		    {
			int x=i+dx[k],y=j+dy[k];
			if(x<0||x>=n||y<0||y>=m||s[x][y]=='X') f(P,P)-=p[k];
			else if(s[x][y]=='.') f(P,x*m+y)-=p[k^1];
		    }
		    if(s[i][j]=='T') f(P,P)=1;
		}
	for(int i=0;i<n*m;++i)
	{
	    int l=max(0,i-m),r=min(n*m-1,i+m),k=-1;
	    for(int j=l;j<=r;++j) if(!vis[j]) if(fabs(f(j,i))>eps) {k=j;break;}
	    if(!~k) {pos[i]=-1;continue;}
	    vis[pos[i]=k]=1;
	    for(int j=k+1;j<=r;++j)
	    {
		db t=f(j,i)/f(k,i);
		for(int x=i;x<n*m&&x<=k+m;++x) f(j,x)-=f(k,x)*t;
		g[j]-=g[k]*t;
	    }
	}
	for(int i=n*m-1;~i;--i)
	    if(~pos[i])
	    {
		int k=pos[i],l=max(0,i-m),r=min(n*m-1,i+m);
		for(int j=l;j<=r;++j) if(j^k) g[k]-=f(k,j)*g[j];
		g[k]/=f(k,i);
	    }
	for(int i=0;i<n;++i) for(int j=0;j<m;++j) if(s[i][j]=='T') printf("%.9f
",g[pos[i*m+j]]);
	puts("");
    }
}