Codeforces-113DMuseum高斯消元求概率

题目链接：https://codeforces.com/problemset/problem/113/D

题目大意：有一个n个点m条边的无向图，两个人分别从a,b出发，每个人每分钟有$p_i$的概率不动, 有$1-p_i$的概率走到随机一个相邻的点。当他们在同一时刻选择前往同一个房间, 他们就会在那个房间相遇并停止。求在每个点相遇的概率。

Examples

Input

Output

0.5000000000 0.5000000000

Input

Output

0.3333333333 0.3333333333 0.1666666667 0.1666666667

做的时候是个sb，写完了感觉是个sb题。。。。。

我们设$P_{i,j}$表示第一个到$i$，第二个人到$j$的概率，那么则有式子：

$P_{x,y}=sum P_{i,j} imes A_{(i,j)(x,y)}$其中$A_{(i,j)(x,y)}$表示从$(i,j)$转移到相邻的$(x,y)$的概率，移项过后可变为：

$-P_{x,y}+sum P_{i,j} imes A_{(i,j)(x,y)}$我们可以将其合并，由于坐标$(i,j)$中不存在$i=j$而$(x,y)$中存在，所以会额外多出n项$-A_{(i,j)(x,x)}$。

其构造代码如下：

int len=n*n-n;
for (int i=1; i<=n; i++) { //i,j到x,y的系数求解
    for (int j=1; j<=n; j++) {
        if (i==j) continue;
        int u=id[i][j];
        a[u][u]=-1;//i,j到i,j的系数,-p_x,y+Σp_i,j*A即i,j->i,j有一个-1的系数
        for (int x=1; x<=n; x++) {
            for (int y=1; y<=n; y++) {
                if (x==y) {//移项后多出的额外的n维
                    a[u][len+x]-=solve(i,j,x,y);
                    continue;
                }
                int v=id[x][y];
                a[u][v]+=solve(i,j,x,y);
            }
        }
    }
}

接下来就是高斯消元了。

以下是AC代码：

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

const double esp=1e-5;

int dis[30][30],du[30],id[30][30];
double p[30],a[500][500];

double solve(int i,int j,int x,int y)
{//i,j直接转移到相邻x,y的概率
    double ans=0;
    if (i==x && j==y) ans=p[i]*p[j];
    else if (i==x && j!=y) ans=p[i]*(1-p[j])/du[j]*dis[j][y];
    else if (i!=x && j==y) ans=p[j]*(1-p[i])/du[i]*dis[i][x];
    else if (i!=x && j!=y) ans=(1-p[i])*(1-p[j])/du[i]/du[j]*dis[i][x]*dis[j][y];
    return ans;
}

int main()
{
    int n,m,st1,st2;
    scanf ("%d%d%d%d",&n,&m,&st1,&st2);
    for (int i=1; i<=m; i++){
        int u,v;
        scanf ("%d%d",&u,&v);
        du[u]++;du[v]++;
        dis[u][v]=dis[v][u]=1;
    }
    for (int i=1; i<=n; i++)
        scanf ("%lf",&p[i]);
    if (st1==st2){
        for (int i=1; i<=n; i++){
            if (i==st1) printf("%.6f%c",1.0,i==n?'
':' ');
            else printf("%.6f%c",0.0,i==n?'
':' ');
        }
        return 0;
    }
    int cnt=0;
    for (int i=1; i<=n; i++)
        for (int j=1; j<=n; j++) 
            if (i!=j)
                id[i][j]=++cnt;

    int len=n*n-n;
    for (int i=1; i<=n; i++){//i,j到x,y的系数求解
        for (int j=1; j<=n; j++){
            if (i==j) continue;
            int u=id[i][j];
            a[u][u]=-1;//i,j到i,j的系数,-p_x,y+Σp_i,j*A即i,j->i,j有一个-1的系数
            for (int x=1; x<=n; x++){
                for (int y=1; y<=n; y++){
                    if (x==y) {//移项后多出的额外的n维
                        a[u][len+x]-=solve(i,j,x,y);
                        continue;
                    }
                    int v=id[x][y];
                    a[u][v]+=solve(i,j,x,y);
                }
            }
        }
    }

    for (int i=1; i<=len; i++){//化为上三角
        if (fabs(a[i][i])<esp){//系数为0
            int pos=-1;
            for (int j=i+1; j<=len; j++)
                if (fabs(a[j][i])>=esp){
                    pos=j;break;
                }
            swap(a[i],a[pos]);
        }

        //消去下面行的该列元素
        for (int j=i+1; j<=len; j++){
            double f=a[j][i]/a[i][i];
            for (int k=i; k<=n*n; k++)
                a[j][k]-=f*a[i][k];
        }
    }

    for (int i=len; i>=1; i--){
        double f=a[i][i];
        for (int j=len+1; j<=n*n; j++)//i,j到k,k在i,j到i,j的系数化为1后的系数
            a[i][j]/=f;

        a[i][i]=1;
        for (int j=i-1; j>=1; j--){//上面的每一行减去该行
            f=a[j][i];//f=a[j][i]/a[i][i]
            for (int k=len+1; k<=n*n; k++)
                a[j][k]-=f*a[i][k];
        }
    }
    for (int i=1; i<=n; i++)
        printf("%.6f%c",a[id[st1][st2]][len+i],i==n?'
':' ');
    return 0;
}

路漫漫兮