1. 程式人生 > >【數學】稀疏圖的隨機遊走問題

【數學】稀疏圖的隨機遊走問題

.... 稀疏 tor 存在 lin 需要 ans 生成 %d

Description

給出一張n個點,m條邊的平面圖,從1號點開始隨機遊走,抵達n號點則結束,問期望步數?
n<=5000

Solution

這題在wxh的IOI2019國家候選隊論文中也提到了

首先考慮平面圖有什麽好性質,它的邊數不會很多!實際上(根據論文),大於等於3個點的平面圖邊數不會超過3n-6,也就是說邊數和點數是同階的。

我們可以將概率寫成數列的形式,實際上它是一個線性遞推
具體來說,我們設遊走在恰好第i步結束的概率是\(f_i\),那麽存在一個長為n,下標從1開始的數列\(a\),使得對於任意\(j\geq n\),有\(f_j=\sum\limits_{k=1}^{n}f_{j-k}a_k\)

這是由於如果我們將轉移寫成矩陣,實際上就是矩陣反復相乘,而對於任意n*n的矩陣M,\(I,M,M^2,M^3,....\)構成n階線性遞推,且系數就是M的特征多項式的系數(符號什麽的移下項就好)。
然後矩陣的某一個元素列出來也構成了線性遞推。

具體可以參考zzq的IOI2019國家候選隊論文。

既然這樣,我們不妨BFS求出第\(1\)到第\(2*n\)步時結束的概率,然後采用Berlekamp-Massey算法求出\(f\)的遞推式,用生成函數寫出來,一定有\(F(x)=F(x)A(x)+R(x)\),其中\(F(x)\)\(f\)的生成函數,\(A(x)\)是遞推系數的生成函數,\(R(x)\)

是一個小於n次的多項式,表示前n項可能不滿足遞推式的補足部分。

我們求出了\(F(x)\)的前\(2n\)項,求出了\(A(x)\),自然也可以求出\(R(x)\)
那麽\(F(x)={R(x)\over 1-A(x)}\)

然而我們要求的是\(\sum f_i*i\)
實際上把\(F(x)=\sum\limits f_ix^i\)求導,\(F'(x)=\sum f_i*ix^{i-1}\),代入\(F'(1)\)就是答案。

實際上也可以不用求導,我們設\(g_i\)為走到第i步還未結束的概率
那麽有\(\sum\limits f_i*i=\sum\limits g_i\)
這樣前面直接暴力算g,然後直接BM算出g的遞推式,然後一樣算答案即可,省去了求導,就不需要f了。

分析復雜度,由於線性遞推是n階的,所以這裏BM算法最壞是\(O(n^2)\)的,由於m與n同階,因此前面的BFS也可以看做是\(O(n^2)\)
因此我們就用\(O(n^2)\)的時間復雜度解決了問題,比\(O(n^3)\)的暴力高斯消元要優秀很多。

Code

#include <bits/stdc++.h>
#define fo(i,a,b) for(int i=a;i<=b;++i)
#define fod(i,a,b) for(int i=a;i>=b;--i)
#define N 5005
#define LL long long
#define mo 998244353
using namespace std;
LL ny[N],f[N],g[N];
int m1,n,m,t,rd[N],n1;
LL ksm(LL k,LL n)
{
    LL s=1;
    for(;n;n>>=1,k=k*k%mo) if(n&1) s=s*k%mo;
    return s;
}
LL ap[2*N],fp[2][N],gp[N];
int d[2][N];
vector<int> a1[N];
/*
inline void inc(LL &x,const LL &v)
{
    x=(x+v)%mo;
}*/
#define inc(x,v) (x=(x+v)%mo)
void bfs()
{
    memset(fp,0,sizeof(fp));
    LL sum=1;
    fp[0][1]=1,ap[0]=1;
    int p,r;LL v;
    vector<int>::iterator it;
    fo(i,0,2*n-1)
    {
        int i1=i&1;
        fo(j,1,n) fp[1^i1][j]=gp[j]=0;
        fo(k,1,n-1)
        {
            if(fp[i1][k]&&k!=n)
            {
                v=fp[i1][k]*ny[rd[k]]%mo*ny[2]%mo;
                for(it=a1[k].begin();it!=a1[k].end();it++)
                {
                    p=*it;
                    inc(fp[1^i1][p],v);
                }
            }
        }
        fo(k,1,n) gp[k]=fp[1^i1][k];
        fo(k,1,n) 
            if(gp[k]) 
            {
                v=gp[k]*ny[rd[k]]%mo;
                for(it=a1[k].begin();it!=a1[k].end();it++)
                {
                    p=*it;
                    inc(fp[1^i1][p],v);
                }
            }
        sum=(sum-fp[1^i1][n]+mo)%mo;
        ap[i+1]=sum;
        fp[1^i1][n]=0;
    }
}

LL rc[4*N],rp[4*N],le,le1,rw[4*N];
void BM()
{
    le=le1=0;
    memset(rc,0,sizeof(rc));
    memset(rp,0,sizeof(rp));
    int lf=0;LL lv=0;
    fo(i,0,n1)
    {
        LL v=0;
        fo(j,1,le) inc(v,rc[j]*ap[i-j]%mo);
        if(v==ap[i]) continue;
        if(le==0) 
        {
            le=i+1;
            fo(j,1,le) rc[j]=rp[j]=0;
            le1=0,lf=i,lv=(ap[i]-v)%mo;
            continue;
        }
        v=(ap[i]-v+mo)%mo;
        LL mul=v*ksm(lv,mo-2)%mo;
        
        fo(j,1,le) rw[j]=rc[j];
        
        inc(rc[i-lf],mul);
        fo(j,i-lf+1,i-lf+le1) inc(rc[j],(mo-mul*rp[j-(i-lf)]%mo)%mo);
        if(le<i-lf+le1)
        {
            swap(le1,le);
            le=i-lf+le,lf=i,lv=v;
            fo(j,1,le1) rp[j]=rw[j];
        }
        
        v=0;
        fo(j,1,le) inc(v,rc[j]*ap[i-j]%mo);
    }
}

int main()
{
    ny[1]=1;
    fo(i,2,5000) ny[i]=(-ny[mo%i]*(LL)(mo/i)%mo+mo)%mo;
    cin>>t;
    while(t--)
    {
        scanf("%d%d",&n,&m);
        fo(i,1,n)
        {
            int x,y;
            scanf("%d%d",&x,&y);
            a1[i].clear();
        }
        memset(rd,0,sizeof(rd));
        m1=0;
        fo(i,1,m)
        {
            int x,y;
            scanf("%d%d",&x,&y);
            a1[x].push_back(y),a1[y].push_back(x); 
            rd[x]++,rd[y]++;
        }
        n1=2*n;
        bfs();
        BM();
        rc[0]=1,rp[0]=0;
        fo(i,1,le) rc[i]=(mo-rc[i])%mo,rp[i]=0;
        fo(i,0,le)
            fo(j,0,n1) if(i+j<=le) inc(rp[i+j],rc[i]*ap[j]%mo);
        LL ans=0,sv=0;
        fo(i,0,le+n1) inc(ans,rp[i]);
        fo(i,0,le) inc(sv,rc[i]);
        printf("%lld\n",ans*ksm(sv,mo-2)%mo);
    }
}

【數學】稀疏圖的隨機遊走問題