好像四联通八联通的网格图消元都可以那么干,以下是四联通的,具体地就是每条方程都只和
(x,y),(x−1,y),(x+1,y),(x,y−1),(x,y+1)
(
x
,
y
)
,
(
x
−
1
,
y
)
,
(
x
+
1
,
y
)
,
(
x
,
y
−
1
)
,
(
x
,
y
+
1
)
四个未知量有关。
对于这样的高斯消元,我们可以选择从
(0,0),(1,0),(2,0)...(0,1)...(n,m)
(
0
,
0
)
,
(
1
,
0
)
,
(
2
,
0
)
.
.
.
(
0
,
1
)
.
.
.
(
n
,
m
)
这样的顺序,也就是从下往上一行一行消,那么消到
(x,y)
(
x
,
y
)
这条方程的时候,与其相关的未知量只有
(x,y),(x+1,y),(x+2,y)...(0,y+1)...(x,y+1)
(
x
,
y
)
,
(
x
+
1
,
y
)
,
(
x
+
2
,
y
)
.
.
.
(
0
,
y
+
1
)
.
.
.
(
x
,
y
+
1
)
这
O(n)
O
(
n
)
个,而且与
(x,y)
(
x
,
y
)
相关的方程也只有上述这些位置的方程。那么我们每次消就只要那
(x,y)
(
x
,
y
)
消
O(n)
O
(
n
)
条方程,而且消每条方程的时候也只需要消
O(n)
O
(
n
)
个位置,那么复杂度就可以从
O(n6)
O
(
n
6
)
优化到
O(n4)
O
(
n
4
)
。
然后圆里的散点也是网格图,这题就做完了。
代码:
#include<iostream>
#include<cstdio>
#include<cstring>
#define N 55
#define ll long long
#define inv(x) (ksm(x,mod-2))
#define up(x,y) (x=(x+(y))%mod)
using namespace std;
const int mod=1000000007;
int R,gx[4]={-1,0,1,0},gy[4]={0,-1,0,1},id[N<<1][N<<1],tot,pl[N*N<<2][2],A[7850][7850];
ll P[4],d[N*N<<2];
ll ksm(ll a,ll b){ll r=1;for(;b;b>>=1){if(b&1)r=r*a%mod;a=a*a%mod;}return r;}
void gauss()
{
int m=0;
for(int i=1;i<=tot;i++)
{
m=max(m,id[pl[i][0]][pl[i][1]+1]);
for(int j=i+1;j<=m;j++)
{
ll w=mod-inv(A[i][i])*A[j][i]%mod;
for(int k=i;k<=m;k++)
up(A[j][k],w*A[i][k]);
up(A[j][tot+1],w*A[i][tot+1]);
}
}
}
void getans()
{
d[tot+1]=1;
for(int i=tot;i;i--)
{
for(int j=i+1;j<=tot+1;j++)
up(d[i],mod-d[j]*A[i][j]%mod);
d[i]=d[i]*inv(A[i][i])%mod;
}
}
int main()
{
scanf("%d%d%d%d%d",&R,&P[0],&P[1],&P[2],&P[3]);
ll s=inv((P[0]+P[1]+P[2]+P[3])%mod);
for(int t=0;t<4;t++)
P[t]=P[t]*s%mod;
for(int y=-R;y<=R;y++)
for(int x=-R;x<=R;x++)
if(x*x+y*y<=R*R) id[x+N][y+N]=++tot,pl[tot][0]=x+N,pl[tot][1]=y+N;
for(int i=1;i<=tot;i++)
{
int x=pl[i][0],y=pl[i][1];
A[i][i]=mod-1;A[i][tot+1]=1;
for(int t=0;t<4;t++)
if(id[x+gx[t]][y+gy[t]]) A[i][id[x+gx[t]][y+gy[t]]]=P[t];
}
gauss();
getans();
printf("%lld",d[id[N][N]]);
return 0;
}