题意:一张 nnn 个点的无向连通图,两个人开始时分别在 a,ba,ba,b。每次在 uuu 时会以 ppp 的概率原地不动,1−p1-p1−p 的概率等概率随机选择到一个相邻的点,当两人在同一点时停止。分别求在每个点相遇的概率。
n≤22n\leq 22n≤22
网上一堆 “从起点走到 (i,j)(i,j)(i,j) 的概率”,看得我一脸懵逼……
很容易分析出转移矩阵 MMM,然后相当于求这个东西:
limt→+∞MtV\lim_{t\to +\infin}M^{t}Vt→+∞limMtV
但这个并没有通用的求法,因为矩阵的特征向量有无数多个。
算法一
比较直观的解法。
设 f(i,j)f(i,j)f(i,j) 表示 (i,j)(i,j)(i,j) 这个状态到达次数的期望,即这个状态在所有世界线中出现次数的平均值。
由于最终点只有可能出现 000 次或 111 次,所以它的期望次数就是概率。
而这个期望随便消一下就可以算出来。
算法二
比较本质的解法。
考虑枚举一个最终状态 (s,s)(s,s)(s,s),在此条件下求 f(i,j)f(i,j)f(i,j) 表示最终到这个状态的概率,令 f(i,i)=[i=s]f(i,i)=[i=s]f(i,i)=[i=s],就可以消元了。但这样是 O(n7)O(n^7)O(n7) 的,无法通过。
考虑一次性把每个点作为终点的 nnn 个答案算出来,即构建出 (n2−n)×(n2)(n^2-n)\times (n^2)(n2−n)×(n2) 的矩阵。这样有 n2n^2n2 个未知数但只有 n2−nn^2-nn2−n 个方程,无法解出,但可以求出 f(a,b)f(a,b)f(a,b) 关于 f(1,1),f(2,2),…,f(n,n)f(1,1),f(2,2),\dots,f(n,n)f(1,1),f(2,2),…,f(n,n) 的线性表达,就可以求出答案了。
复杂度 O(n6)O(n^6)O(n6)
所以算法一算法二以及上面那个假算法写出来都一样的……
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#include <cmath>
#include <vector>
using namespace std;
vector<int> e[25];
double p[25],a[505][505];
int n,m,sa,sb;
inline int id(int x,int y)
{if (x==y) return n*n-n+x;return (x-1)*(n-1)+y-(y>x);
}
void gauss(int n,int m)
{for (int i=1;i<=n;i++){int pos=i;for (int j=i+1;j<=n;j++) if (fabs(a[j][i])>fabs(a[pos][i])) pos=j;if (pos>i) swap(a[i],a[pos]);for (int j=1;j<=n;j++)if (j!=i){double t=a[j][i]/a[i][i];for (int k=i;k<=m;k++)a[j][k]-=t*a[i][k];}}
}
int main()
{scanf("%d%d%d%d",&n,&m,&sa,&sb);if (sa==sb){for (int i=1;i<=n;i++) if (i==sa) printf("%.10f ",1.0);else printf("%.10f ",0.0);return 0;}for (int i=1;i<=m;i++){int u,v;scanf("%d%d",&u,&v);e[u].push_back(v),e[v].push_back(u);}for (int i=1;i<=n;i++) scanf("%lf",&p[i]);for (int u=1;u<=n;u++)for (int v=1;v<=n;v++)if (u!=v){int s=id(u,v);double pu=1.0/e[u].size(),pv=1.0/e[v].size();for (int i=0;i<=(int)e[u].size();i++)for (int j=0;j<=(int)e[v].size();j++){double t=(i<(int)e[u].size()? (1-p[u])*pu:p[u])*(j<(e[v].size())? (1-p[v])*pv:p[v]);a[s][id(i<(int)e[u].size()? e[u][i]:u,j<(int)e[v].size()? e[v][j]:v)]=t;}a[s][s]-=1;}gauss(n*n-n,n*n);int s=id(sa,sb);for (int i=n*n-n+1;i<=n*n;i++) printf("%.10f ",-a[s][i]/a[s][s]);return 0;
}