[BZOJ3270]博物馆 (高斯消元)

发布于 2018-05-23  53 次阅读


本文章由SYCstudio或本站其它作者所有,严禁转载,转载必究

本文链接地址:[BZOJ3270]博物馆 (高斯消元)

Description

有一天Petya和他的朋友Vasya在进行他们众多旅行中的一次旅行,他们决定去参观一座城堡博物馆。这座博物馆有着特别的样式。它包含由m条走廊连接的n间房间,并且满足可以从任何一间房间到任何一间别的房间。
两个人在博物馆里逛了一会儿后两人决定分头行动,去看各自感兴趣的艺术品。他们约定在下午六点到一间房间会合。然而他们忘记了一件重要的事:他们并没有选好在哪儿碰面。等时间到六点,他们开始在博物馆里到处乱跑来找到对方(他们没法给对方打电话因为电话漫游费是很贵的)
不过,尽管他们到处乱跑,但他们还没有看完足够的艺术品,因此他们每个人采取如下的行动方法:每一分钟做决定往哪里走,有Pi 的概率在这分钟内不去其他地方(即呆在房间不动),有1-Pi 的概率他会在相邻的房间中等可能的选择一间并沿着走廊过去。这里的i指的是当期所在房间的序号。在古代建造是一件花费非常大的事,因此每条走廊会连接两个不同的房间,并且任意两个房间至多被一条走廊连接。
两个男孩同时行动。由于走廊很暗,两人不可能在走廊碰面,不过他们可以从走廊的两个方向通行。(此外,两个男孩可以同时地穿过同一条走廊却不会相遇)两个男孩按照上述方法行动直到他们碰面为止。更进一步地说,当两个人在某个时刻选择前往同一间房间,那么他们就会在那个房间相遇。
两个男孩现在分别处在a,b两个房间,求两人在每间房间相遇的概率。

Tag

高斯消元

解决思路

设$$F[i][j]$$表示$A$在$i$号点,$B$在$j$号点的概率,那么$i$和$j$分别有不动和从旁边转移过来两种方式。列出$O(n^2)$个方程,消元求解。
需要注意的是,不能从两个相同的点转移过来,因为这样表示已经达到目标状态了。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;

#define ll long long
#define ld long double
#define mem(Arr,x) memset(Arr,x,sizeof(Arr))

const int maxN=23;
const int maxM=maxN*maxN;
const ld eps=1e-13;
const int inf=2147483647;

int n,m;
int S,T;
ld P[maxN];
int edgecnt=0,Head[maxN],Next[maxM],V[maxM];
int Id[maxN][maxN];
ld Degree[maxN],To[maxN];
ld Mat[maxN*maxN][maxN*maxN];

int main()
{
    mem(Head,-1);
    scanf("%d%d%d%d",&n,&m,&S,&T);
    for (int i=1;i<=m;i++)
    {
        int u,v;scanf("%d%d",&u,&v);
        edgecnt++;Next[edgecnt]=Head[u];Head[u]=edgecnt;V[edgecnt]=v;
        edgecnt++;Next[edgecnt]=Head[v];Head[v]=edgecnt;V[edgecnt]=u;
        Degree[u]+=1.0;Degree[v]+=1.0;
    }
    for (int i=1;i<=n;i++) scanf("%LF",&P[i]);
    for (int i=1;i<=n;i++) To[i]=(1.0-P[i])/Degree[i];

    int idcnt=0;
    for (int i=1;i<=n;i++) for (int j=1;j<=n;j++) Id[i][j]=++idcnt;

    for (int u=1;u<=n;u++)
        for (int v=1;v<=n;v++)
        {
            if (u!=v) Mat[Id[u][v]][Id[u][v]]=(P[u]*P[v]-1.0);
            else Mat[Id[u][v]][Id[u][v]]=1.0;
            for (int e1=Head[u];e1!=-1;e1=Next[e1])
                for (int e2=Head[v];e2!=-1;e2=Next[e2])
                {
                    if ((u!=V[e2])&&(v!=V[e2])) Mat[Id[u][v]][Id[u][V[e2]]]=P[u]*(1.0-P[V[e2]])/(ld)Degree[V[e2]];
                    if ((u!=V[e1])&&(v!=V[e1])) Mat[Id[u][v]][Id[V[e1]][v]]=(1.0-P[V[e1]])*P[v]/(ld)Degree[V[e1]];
                    if (V[e1]!=V[e2]) Mat[Id[u][v]][Id[V[e1]][V[e2]]]=(1.0-P[V[e1]])*(1.0-P[V[e2]])/(ld)Degree[V[e1]]/(ld)Degree[V[e2]];
                }
        }
    Mat[Id[S][T]][idcnt+1]=1;
    for (int i=1;i<=idcnt;i++)
    {
        int j=i;
        while (fabs(Mat[j][i])<eps) j++;
        for (int k=1;k<=idcnt+1;k++) swap(Mat[j][k],Mat[i][k]);
        ld d=(ld)1.0/Mat[i][i];
        for (int k=1;k<=idcnt+1;k++) Mat[i][k]*=d;
        for (int k=1;k<=idcnt;k++)
            if ((k!=i)&&(fabs(Mat[k][i])>eps))
            {
                d=Mat[k][i]/Mat[i][i];
                for (int l=1;l<=idcnt+1;l++) Mat[k][l]=Mat[k][l]-Mat[i][l]*d;
            }
    }
    ld sum=0;
    for (int i=1;i<=n;i++) sum=sum+Mat[Id[i][i]][idcnt+1];
    for (int i=1;i<=n;i++) printf("%.6LF ",Mat[Id[i][i]][idcnt+1]/sum);
    return 0;
}

本文章由SYCstudio或本站其它作者所有,严禁转载,转载必究

本文链接地址:[BZOJ3270]博物馆 (高斯消元)


HNCJ OIer 一枚