[Codeforces963] E-Circles of Waiting 网格图高斯消元

好像四联通八联通的网格图消元都可以那么干,以下是四联通的,具体地就是每条方程都只和 (x,y),(x1,y),(x+1,y),(x,y1),(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;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值