BZOJ 3270 博物馆 ——概率DP 高斯消元

用$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]]);
}