BZOJ 3270 && BZOJ 1778 (期望DP && 高斯消元)

BZOJ 3270 :设置状态为Id(x,y)表示一人在x，一人在y这个状态的概率。 所以总共有n^2种状态。

p[i]表示留在该点的概率，Out[i]=(1-p[i])/Degree[i]表示离开该点的概率.

那么对于每一种状态a，b 则有P(a,b)=p[a]∗p[b]∗P(a,b)+Out[u]∗p[b]∗P(u,b)+p[a]∗Out[v]∗P(a,v)+Out[u]∗Out[v]∗P(u,v) 则有n^2个方程

对于起始状态a，b，则有P(a,b)=p[a]∗p[b]∗P(a,b)+Out[u]∗p[b]∗P(u,b)+p[a]∗Out[v]∗P(a,v)+Out[u]∗Out[v]∗P(u,v)+1

 
#include <cstdio>
#include <cstring>
const int Maxn=10010;
double p[Maxn],Out[Maxn],a[510][510];
int Degree[Maxn],n,m,head[Maxn],cnt,s,t,u,v;
struct Edge{int to,next;}edge[Maxn<<1];
inline void Add(int u,int v) 
{edge[cnt].to=v;edge[cnt].next=head[u];head[u]=cnt++;}
inline int Id(int x,int y) {return (x-1)*n+y;}
inline void Swap(double &x,double &y) {double t=x;x=y;y=t;}
void Build(int x,int y)
{
    int S=Id(x,y); 
    a[S][S]-=1;
    for (int i=head[x];i!=-1;i=edge[i].next)
        for (int j=head[y];j!=-1;j=edge[j].next)
        {
            int u=edge[i].to,v=edge[j].to,T=Id(u,v);
            if (u!=v) 
            {
                if (u==x && v==y) a[S][T]+=p[u]*p[v];
                else if (u==x) a[S][T]+=p[u]*Out[v];
                else if (v==y) a[S][T]+=Out[u]*p[v];
                else a[S][T]+=Out[u]*Out[v];
            }
        }
}
inline void Guass(int n)
{
    for (int i=1;i<=n;i++)
    {
        for (int j=i;j<=n;j++) 
            if (a[j][i])
            {
                for (int k=1;k<=n+1;k++) Swap(a[i][k],a[j][k]);
                for (int k=1;k<=n+1;k++) if (k!=i) a[i][k]/=a[i][i];
                a[i][i]=1;
                break;
            }
        if (!a[i][i]) continue;
        for (int j=1;j<=n;j++) 
        {
            if (j==i) continue;
            double t=a[j][i];
            for (int k=1;k<=n+1;k++) a[j][k]-=t*a[i][k];
        }
    }
}
int main()
{
    scanf("%d%d%d%d",&n,&m,&s,&t);
    memset(head,-1,sizeof(head));
    for (int i=1;i<=m;i++)
    {
        scanf("%d%d",&u,&v);
        Add(u,v),Add(v,u);
        Degree[u]++,Degree[v]++;
    }
    for (int i=1;i<=n;i++) scanf("%lf",&p[i]);
    for (int i=1;i<=n;i++) Out[i]=(1-p[i])/Degree[i];
    for (int i=1;i<=n;i++) Add(i,i);
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++) Build(i,j);
    a[Id(s,t)][n*n+1]=-1;
    Guass(n*n);
    for (int i=1;i<=n;i++) printf("%.6lf ",a[Id(i,i)][n*n+1]);
    return 0;
}

 BZOJ 1778:和上面的题目差不多但是要注意的是自己不可能转移到自己因为要么炸要么移动

P(u)=Out(v)*P(v) 特别的 P(1)=Out(v)*P(v)+P/Q然后高斯消元即可.

#include <cstdio>
#include <cmath>
#include <cstring>
#define LD long double
const int Maxn=510;
const int Maxm=50100;
const LD eps=1e-13;
int Degree[Maxn],head[Maxn],n,m,P,Q,cnt,u,v;
LD a[Maxn][Maxn],t,Out[Maxn];
struct  Edge{int to,next;}edge[Maxm<<2];
inline void Swap(LD &x,LD &y) {LD t=x;x=y;y=t;}
inline int dcmp(LD x) {if (fabs(x)<eps) return 0; return (x>0)?1:-1;}
inline void Add(int u,int v)
{edge[cnt].to=v;edge[cnt].next=head[u];head[u]=cnt++;}
inline void Build(int u)
{
    a[u][u]=1.0;
    for (int i=head[u];i!=-1;i=edge[i].next)
        a[u][edge[i].to]=-Out[edge[i].to];
}
inline void Gauss()
{
    for (int i=1;i<=n;i++)
    {
        for (int j=i;j<=n;j++) 
            if (dcmp(a[i][j])!=0)
            {
                for (int k=1;k<=n+1;k++) Swap(a[i][k],a[j][k]);
                for (int k=1;k<=n+1;k++) if (k!=i) a[i][k]/=a[i][i];
                a[i][i]=1;
                break;
            }
        if (dcmp(a[i][i])==0) continue;
        for (int j=1;j<=n;j++) 
        {
            if (j==i) continue;
            LD t=a[j][i];
            for (int k=1;k<=n+1;k++) a[j][k]-=t*a[i][k];
        }
    }
}
int main()
{
    scanf("%d%d%d%d",&n,&m,&P,&Q); t=(LD)P/(LD)Q; t=1.0-t;
    memset(head,-1,sizeof(head)); cnt=0;
    for (int i=1;i<=m;i++) 
    {
        scanf("%d%d",&u,&v);
        Add(u,v),Add(v,u),Degree[u]++,Degree[v]++;
    }
      
    for (int i=1;i<=n;i++) Out[i]=t/(LD)Degree[i];
    for (int i=1;i<=n;i++) Build(i);
    a[1][n+1]=1-t;
    Gauss();
    for (int i=1;i<=n;i++) printf("%.9Lf\n",a[i][n+1]);
    return 0;
}