用$F(i,j)$表示A在i，B在j的概率。
然后很容易列出转移方程。
然后可以高斯消元了！
被一个问题困扰了很久，为什么起始点的概率要加上1。
（因为其他博客上都是直接写成-1，雾）
考虑初始状态是由什么转移过来的，发现可以由其他点走过来，也可以由初始定义转移。
而初始的定义就决定了它有一个1的概率。
加上之后就可以高斯消元了
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define F(i,j,k) for (int i=j;i<=k;++i)
#define D(i,j,k) for (int i=j;i>=k;--i)
double a[405][405],p[21];
int id[21][21],n,m,x,y,du[21],tot,st;
int h[405],to[405],ne[405],en=0;
void add(int a,int b)
{
to[en]=b;ne[en]=h[a];h[a]=en++;
}
void solve(int x,int y)
{
int num=id[x][y];
a[num][num]=1.0;
for (int i=h[x];i>=0;i=ne[i])
for (int j=h[y];j>=0;j=ne[j])
{
int tx=to[i],ty=to[j]; if (tx==ty) continue;
int tnum=id[tx][ty];
if (tx==x&&ty==y) a[num][tnum]-=p[tx]*p[ty];
else if (tx==x) a[num][tnum]-=p[tx]*(1-p[ty])/(du[ty]*1.0);
else if (ty==y) a[num][tnum]-=p[ty]*(1-p[tx])/(du[tx]*1.0);
else a[num][tnum]-=(1-p[tx])*(1-p[ty])/(du[tx]*1.0)/(du[ty]*1.0);
}
}
void gauss()
{
F(i,1,tot)
{
int tmp=i;
while (!a[tmp][i]&&tmp<=tot) tmp++;
if (tmp>tot) continue;//自由元
F(j,i,tot+1) swap(a[i][j],a[tmp][j]);
F(j,1,tot) if (j!=i)
{
double t=a[j][i]/a[i][i];
F(k,1,tot+1)
a[j][k]-=t*a[i][k];
}
}
}
int main()
{
memset(h,-1,sizeof h);
scanf("%d%d%d%d",&n,&m,&x,&y);
F(i,1,m)
{
int a,b;
scanf("%d%d",&a,&b);
add(a,b);add(b,a);
du[a]++;du[b]++;
}
F(i,1,n) add(i,i);
F(i,1,n) F(j,1,n) id[i][j]=++tot;
a[id[x][y]][tot+1]=1;
F(i,1,n) scanf("%lf",&p[i]);
F(i,1,n) F(j,1,n) solve(i,j);
gauss();
F(i,1,n)
{
int t=id[i][i];
printf("%.6f ",a[t][tot+1]/a[t][t]);
}
printf("
");
}
我怎么还是不会写高斯消元？(⊙_⊙?)
方法二：
考虑初始向量S，转移矩阵A
那么答案就是$ans=SI+SA+SA^2...$
所以$ans=S(I-A)^{-1}$
然后就有$ans[I-A]^{T}=S$
高斯消元即可，或者求逆之后矩阵乘法
#include <map>
#include <ctime>
#include <cmath>
#include <queue>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define F(i,j,k) for (int i=j;i<=k;++i)
#define D(i,j,k) for (int i=j;i>=k;--i)
#define maxn 405
double f[402][402],p[402],d[402][402];
double ans[402],S[402];
int ed,n,m,sa,sb,cnt;
int id[21][21],du[maxn],h[maxn],to[maxn],ne[maxn],en=0;
void add(int a,int b)
{to[en]=b;ne[en]=h[a];h[a]=en++;}
void Build()
{
ed=cnt;
F(i,1,n) F(j,1,n)
if (i!=j) F(k,1,n) F(l,1,n)
f[id[i][j]][id[k][l]]-=d[i][k]*d[j][l];
F(i,1,ed) F(j,1,i-1) swap(f[i][j],f[j][i]);
F(i,1,ed) f[i][i]+=1;
f[id[sa][sb]][ed+1]=1;
}
void Gauss()
{
F(i,1,ed)
{
int tmp=i;
F(j,i+1,ed) if (f[j][i]>f[tmp][i]) tmp=j;
if (tmp!=i) F(j,i,ed+1) swap(f[tmp][j],f[i][j]);
F(j,i+1,ed)
{
double t=f[j][i]/f[i][i];
F(k,i,ed+1) f[j][k]-=t*f[i][k];
}
}
D(i,ed,1)
{
F(j,i+1,ed) f[i][ed+1]-=f[i][j]*ans[j],f[i][j]=0;
ans[i]=f[i][ed+1]/f[i][i];
}
}
int main()
{
memset(h,-1,sizeof h);
scanf("%d%d%d%d",&n,&m,&sa,&sb);
F(i,1,m)
{
int a,b;
scanf("%d%d",&a,&b);
add(a,b);add(b,a);
du[a]++;du[b]++;
}
F(i,1,n) scanf("%lf",&p[i]);
F(i,1,n)
{
d[i][i]=p[i];
for (int j=h[i];j>=0;j=ne[j])
d[i][to[j]]+=(1-p[i])/du[i];
}
F(i,1,n) F(j,1,n) id[i][j]=++cnt;
Build();
Gauss();
F(i,1,n) printf("%.6f ",ans[id[i][i]]);
}