稀疏图的随机游走问题

最新推荐文章于 2023-07-19 17:02:36 发布

BAJim_H

最新推荐文章于 2023-07-19 17:02:36 发布

阅读量869

点赞数 1

分类专栏：题解好题 ---DP ————递推 ————概率与期望 ---线性代数文章标签：概率与期望线性代数

本文链接：https://blog.csdn.net/hzj1054689699/article/details/90046540

版权

题解同时被 3 个专栏收录

353 篇文章 0 订阅

订阅专栏

好题

37 篇文章 0 订阅

订阅专栏

---DP

29 篇文章 0 订阅

订阅专栏

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)$ 的前 $2 n$ 项，求出了 $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);
	}
}