图上随机游走

Posted 2021-06-13 __Archer

这两天做了几道图上随机游走的题，虽然这部分题大多是一个固定的套路，但是还是值得记录一下。
值得一提的是，发生随机游走的局面往往是在一个平等的关系下，即可以反复横跳，这一点在第三题有所体现，也可能是随机游走的本质。
第一题:$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;
}

View Code

 

第二题:$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;
}

View Code

 

 

第三题:$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;
}

View Code

 

 

第四题:$2019icpc沈阳B$