Circles of Waiting
求一个整点四连通随机游⾛，离原点距离超过R期望步数。R≤50。
带状矩阵法
本质上就是网格图的随机游走。
[E_x=sum_y P_{x,y}E_y+1
]
相关联的变量较少，所以使用Band Matrix即可。时间复杂度 (O(R^4))。
https://blog.csdn.net/lycheng1215/article/details/80180178
貌似这题不写主元系数为0时的操作也是对的。
CO int dx[4]={-1,0,1,0},dy[4]={0,-1,0,1};
int p[4];
CO int N=103,O=51;
int idx[N][N],tot;
pair<int,int> pt[N*N];
CO int M=7845+5;
int a[M][M];
int main(){
int R=read<int>();
int all=0;
for(int i=0;i<=3;++i) all+=read(p[i]);
all=fpow(all,mod-2);
for(int i=0;i<=3;++i) p[i]=mul(p[i],all);
for(int y=-R;y<=R;++y)for(int x=-R;x<=R;++x)
if(x*x+y*y<=R*R) idx[x+O][y+O]=++tot,pt[tot]=make_pair(x+O,y+O);
for(int i=1;i<=tot;++i){
a[i][i]=1,a[i][tot+1]=1;
int x=pt[i].first,y=pt[i].second;
for(int j=0;j<=3;++j)if(idx[x+dx[j]][y+dy[j]])
a[i][idx[x+dx[j]][y+dy[j]]]=mod-p[j];
}
int lim=0;
for(int i=1;i<=tot;++i){
lim=max(lim,idx[pt[i].first][pt[i].second+1]);
int inv=fpow(a[i][i],mod-2);
for(int j=i+1;j<=lim;++j)if(a[j][i]){
int coef=mul(mod-a[j][i],inv);
for(int k=i;k<=lim;++k) a[j][k]=add(a[j][k],mul(coef,a[i][k]));
a[j][tot+1]=add(a[j][tot+1],mul(coef,a[i][tot+1]));
}
}
for(int i=tot;i>=1;--i){
for(int j=i+1;j<=tot;++j)
a[i][tot+1]=add(a[i][tot+1],mod-mul(a[i][j],a[j][tot+1])),a[i][j]=0;
a[i][tot+1]=mul(a[i][tot+1],fpow(a[i][i],mod-2)),a[i][i]=1;
}
printf("%d
",a[idx[O][O]][tot+1]);
return 0;
}
主元法
可以对每一行最左边的元素设置未知变量。
每个点可以向上下左右走，那么如何转移系数呢？
[E_{x-1,y}=p_1E_{x-2,y}+p_2E_{x-1,y-1}+p_3E_{x,y}+p_4E_{x-1,y+1}+1\
E_{x,y}=frac{1}{p_3}(E_{x-1,y}-p_1E_{x-2,y}-p_2E_{x-1,y-1}-p_4E_{x-1,y+1}-1)
]
所以说系数的转移还可以间接来整。
时间复杂度 (O(R^3))。
CO int N=100+10;
int p[4];
int c[N][N][N];
int a[N][N];
int main(){
int R=read<int>();
int all=0;
for(int i=0;i<=3;++i) all+=read(p[i]);
all=fpow(all,mod-2);
for(int i=0;i<=3;++i) p[i]=mul(p[i],all);
int n=2*R+1;
function<bool(int,int)> check=[R](int x,int y)->bool{
return (x-R-1)*(x-R-1)+(y-R-1)*(y-R-1)<=R*R;
};
for(int x=1;x<=n;++x)for(int y=1;y<=n;++y){
if(!check(x,y)) continue;
if(!check(x-1,y)){
c[x][y][y]=1;
continue;
}
int inv=fpow(p[2],mod-2);
for(int i=1;i<=n+1;++i){
int sum=c[x-1][y][i];
sum=add(sum,mod-mul(p[0],c[x-2][y][i]));
sum=add(sum,mod-mul(p[1],c[x-1][y-1][i]));
sum=add(sum,mod-mul(p[3],c[x-1][y+1][i]));
if(i==n+1) sum=add(sum,mod-1);
c[x][y][i]=mul(sum,inv);
}
}
for(int y=1;y<=n;++y){
int x=sqrt(R*R-(y-R-1)*(y-R-1))+R+1;
for(int i=1;i<=n+1;++i){
int sum=c[x][y][i];
sum=add(sum,mod-mul(p[0],c[x-1][y][i]));
sum=add(sum,mod-mul(p[1],c[x][y-1][i]));
sum=add(sum,mod-mul(p[3],c[x][y+1][i]));
if(i==n+1) sum=add(sum,mod-1);
a[y][i]=sum;
}
a[y][n+1]=mod-a[y][n+1];
}
for(int i=1;i<=n;++i){
int p=i;
for(int j=i;j<=n;++j)
if(a[j][i]) {p=j;break;}
if(p!=i) swap(a[p],a[i]);
int inv=fpow(a[i][i],mod-2);
for(int j=1;j<=n;++j)if(j!=i){
int coef=mul(mod-a[j][i],inv);
for(int k=i;k<=n+1;++k) a[j][k]=add(a[j][k],mul(coef,a[i][k]));
}
}
for(int i=1;i<=n;++i)
a[i][n+1]=mul(a[i][n+1],fpow(a[i][i],mod-2)),a[i][i]=1;
int ans=c[R+1][R+1][n+1];
for(int i=1;i<=n;++i) ans=add(ans,mul(c[R+1][R+1][i],a[i][n+1]));
printf("%d
",ans);
return 0;
}