这两天做了几道图上随机游走的题，虽然这部分题大多是一个固定的套路，但是还是值得记录一下。
值得一提的是，发生随机游走的局面往往是在一个平等的关系下，即可以反复横跳，这一点在第三题有所体现，也可能是随机游走的本质。
第一题:$bzoj\;3143$
题意:一个无向图，从$1$号节点随机游走，走到$n$号节点停止，问每条边被经过的期望次数。
$solution:$一个经典的套路，将边走过的期望转化为两个端点的期望，设经过点$u,v$的期望值为$f(u),f(v)$,那么边$e$的期望是$\frac{f(u)}{deg(u)}+\frac{f(v)}{deg(v)}$。
至于点期望的求法，设$p_{u,v}$表示若$v$与$u$相邻，那么从$v$走到$u$的概率，即$\frac{1}{deg(v)}$,否则是$0$,那么对于点$u$,可以列出等式$f(u)=p_{u,1}*f(1)+p_{u,2}*f(2)+\dots+p_{u,n}*f(n)$,因为题目说明走到$n$停止，所以$f(n)=0$,同理，$f(1)$本身应该$+1$，因为一开始就在$1$号点。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod = 1000000007;
const int maxn = 505;
const double eps = 1e-9;
vector<int>e[maxn*maxn];
int deg[maxn];
double a[maxn][maxn],p[maxn];
pair<int,int>edge[maxn*maxn];
double val[maxn*maxn];
void add(int u,int v) {
e[u].push_back(v);
}
void Gauss(int n,double a[maxn][maxn]) {
for(int i=0;i<n;i++) {
int mx = i;
for(int j=i+1;j<n;j++) {
if(fabs(a[j][i])>fabs(a[mx][i])) {
mx = j;
}
}
for(int j=0;j<=n;j++) swap(a[i][j],a[mx][j]);
assert(fabs(a[i][i]) > eps);
for(int j=0;j<n;j++) {
if(j!=i) {
double tmp = a[j][i]/a[i][i];
for(int k=i+1;k<=n;k++) {
a[j][k]-=a[i][k]*tmp;
}
}
}
}
for(int i=0;i<n;i++) a[i][n]/=a[i][i];
}
int main() {
int n,m;scanf("%d%d",&n,&m);
for(int i=1;i<=m;i++) {
int u,v;scanf("%d%d",&u,&v);
u--;v--;
if(v>u) swap(v,u);
deg[u]++;deg[v]++;
add(u,v);add(v,u);
edge[i] = make_pair(v,u);
}
a[0][n]=-1.0;
for(int i=0;i<n;i++) a[i][i]=-1.0;
for(int u=0;u<n;u++) {
if(u==n-1) continue;
for(auto v : e[u]) {
a[v][u] += 1.0/deg[u];
}
}
Gauss(n,a);
for(int i=0;i<n;i++) p[i] = a[i][n];
for(int i=1;i<=m;i++) {
int v = edge[i].first,u = edge[i].second;
val[i] = 1.0/deg[v] * p[v];
if(u!=n-1) val[i] += 1.0/deg[u] * p[u];
}
sort(val+1,val+m+1);
double ans = 0;
for(int i=1;i<=m;i++) {
ans += val[i] * (m+1-i);
}
printf("%.3f\n",ans);
return 0;
}
第二题:$bzoj 2337$
题意:无向图，有边权，从$1$到$n$，问走过路径的异或和的期望。
$solution:$拆位，对每一位计算对答案的贡献。假设在考虑第$k$位，设$f(u)$表示从$u$到$n$的异或和期望，因为异或的性质，可以列出：$f(u) = \sum\limits_{v\in son(u)}^{}[w_{v,u}==0]\frac{f(v)}{deg(v)}+[w_{v,u}==1]\frac{1-f(v)}{deg(v)}$,需要注意$f(n)=0$
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 115;
const double eps = 1e-8;
struct node {
int nxt,to,w;
};
node e[maxn*maxn*2];
int head[maxn<<1],tot;
double a[maxn][maxn];
int deg[maxn];
void add(int u,int v,int w) {
e[tot].w = w;e[tot].to = v;e[tot].nxt = head[u];head[u]=tot++;
}
void Gauss(int n,double a[maxn][maxn]) {
//1~n+1
for(int i=1;i<=n;i++) {
int mx = i;
for(int j=i+1;j<=n;j++) {
if(fabs(a[j][i])>fabs(a[mx][i])) {
mx=j;
}
}
for(int j=1;j<=n+1;j++) swap(a[i][j],a[mx][j]);
assert(fabs(a[i][i]) > eps);
for(int j=1;j<=n;j++) {
if(j!=i) {
double tmp = a[j][i]/a[i][i];
for(int k=i+1;k<=n+1;k++) {
a[j][k]-=a[i][k]*tmp;
}
}
}
}
for(int i=1;i<=n;i++) a[i][n+1]/=a[i][i];
}
int main() {
memset(head,-1,sizeof(head));
int n,m;scanf("%d%d",&n,&m);
for(int i=1;i<=m;i++) {
int u,v,w;scanf("%d%d%d",&u,&v,&w);
add(u,v,w);
if(u!=v) add(v,u,w),deg[u]++,deg[v]++;
else deg[u]++;
}
double ans = 0.0;
for(int k=0;k<=30;k++) {
memset(a,0,sizeof(a));
for(int u=1;u<=n-1;u++) {
a[u][u] = 1.0*deg[u];
for(int i=head[u];~i;i=e[i].nxt) {
int v = e[i].to,w = e[i].w;
if((w>>k) & 1) {
a[u][v]+=1.0;a[u][n+1]+=1.0;
}
else a[u][v]+=-1.0;
}
}
a[n][n]=1.0;
Gauss(n,a);
ans += 1.0*(1<<k)*a[1][n+1];
//printf("k = %d,delta = %.3f\n",k, 1.0*(1<<k)*a[1][n+1]);
}
printf("%.3f\n",ans);
return 0;
}
第三题:$bzoj 3640$
这题值得思考。
$solution:$设$dp(u,h)$表示当前走到了$u$点，还剩$h$点生命值时的概率，那么转移方程很容易列出:$dp(u,h)=\sum\limits_{v\in son(u)}^{}\frac{dp(v,h+a_u)}{deg(v)}$，边界$dp(1,hp)=1.0$.答案是$\sum\limits_{h=1}^{h\le hp}dp(n,h)$
看起来这似乎就直接遍历图一遍就可以得到答案了，但这是错误的。原因在于$a_i$可以为$0$.
若$a_u=0$那么转移方程变成$dp(u,h)=\sum\limits_{v\in son(u)}^{}\frac{dp(v,h)}{deg(v)}$,会发现$dp(u,h)$与$dp(v,h)$对等了，$v,u$两个点可以来回摩擦了，那直接推就是错误的。如果$a_i \gt 0$，那么每个$dp(u,h)$都会从$dp(v,h+a_u)$转移而来，这是不可逆的，或者说这是一种指向关系，不可来回摩擦。
显然，对于可以来回摩擦的点，这就是随机游走了，对于不可来回摩擦的点，可以直接转移。
但是生命值有$1e4$，直接做的时间复杂度是$O(n^3*hp)$。考虑先把方程组列出来想优化方法。
发现在枚举生命值$h$的过程中，方程组的系数矩阵始终不变，常数矩阵在变。
于是可以:$\boldsymbol{A*X=B} \iff \boldsymbol{X=A^{-1}B}$
其$A$是系数矩阵，$B$是常数矩阵，矩阵求逆预处理$A$,就可以优化到$O(n^2*hp)$了。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 155;
const double eps = 1e-8;
struct node {
int nxt,to;
};
node e[5555*2];
int head[maxn],tot;
int val[maxn],deg[maxn];
double a[maxn][maxn*2];
double dp[maxn][10005];
double delta[maxn];
void add(int u,int v) {
e[tot].to=v;e[tot].nxt=head[u];head[u]=tot++;
}
void Matrix_inv(int n,double a[maxn][maxn*2]) {
for(int i=1,r;i<=n;i++) {
r=i;
for(int j=i+1;j<=n;j++) {
if(a[j][i]>a[r][i]) r=j;
}
if(r!=i) swap(a[i],a[r]);
assert(fabs(a[i][i]) > eps);
double kk = 1.0/a[i][i];
for(int k=1;k<=n;k++) {
if(k==i) continue;
double p=a[k][i]*kk;
for(int j=i;j<=(n<<1);j++) {
a[k][j]-=p*a[i][j];
}
}
for(int j=1;j<=(n<<1);j++) a[i][j] = a[i][j]*kk;
}
}
void init(int n) {
for(int u=1;u<=n;u++) {
a[u][u]=1.0;a[u][n+u]=1.0;
if(val[u]!=0) continue;
for(int i=head[u];~i;i=e[i].nxt) {
int v = e[i].to;
if(v!=n ) {
a[u][v] += -1.0/deg[v];
}
}
}
Matrix_inv(n,a);
}
int main() {
tot=0;
memset(head,-1,sizeof(head));
int n,m,hp;scanf("%d%d%d",&n,&m,&hp);
for(int i=1;i<=n;i++) scanf("%d",&val[i]);
for(int i=1;i<=m;i++) {
int u,v;scanf("%d%d",&u,&v);
if(u!=v) {
add(u,v);add(v,u);
deg[u]++;deg[v]++;
}
else {
add(u,v);deg[u]++;
}
}
init(n);
double ans=0;
for(int h=hp;h>0;h--) {
memset(delta,0,sizeof(delta));
if(h==hp) delta[1]=1.0;
for(int u=1;u<=n;u++) {
if(val[u]==0) continue;
if(val[u]+h>hp) continue;
for(int i=head[u];~i;i=e[i].nxt) {
int v = e[i].to;
if(v!=n) delta[u] += dp[v][h+val[u]]/deg[v];
}
}
for(int i=1;i<=n;i++) {
for(int j=1;j<=n;j++) {
dp[i][h] += delta[j] * a[i][n+j];
}
}
ans += dp[n][h];
}
printf("%.8f\n",ans);
return 0;
}
第四题:$2019icpc沈阳B$
