bzoj 1778 [Usaco2010 Hol]Dotp 驱逐猪猡（高斯消元）

【题意】

    炸弹从1开始运动，每次有P/Q的概率爆炸，否则等概率沿边移动，问在每个城市爆炸的概率。

【思路】

    设M表示移动一次后i->j的概率。Mk为移动k次后的概率，则有：

        Mk=M^k

    设S={ 1,0,0,0,… }

    设pi为移动i步后到对应点爆炸的概率矩阵，则有：

　　p0=P/Q * S

　　p1=P/Q * S * M1

       　　 …

　　p+oo=P/Q * S * Mn

    则答案为：sigma{ pi },0<=i<+oo

    即：

        Ans=P/Q * S * sigma{ M^i } ,0<=i<+oo

    根据等差数列的求和公式：

        Ans= P/Q * S * (I-M^(+oo))/(I-M)

           = P/Q * S * I/(I-M)

　　  =>  Ans(I-M) = P/Q * S

    其中I为单位矩阵，I[j][j]=1，其余为0。

    于是可以得到n个线性方程，可以使用高斯消元法求解。

　 注意Ans*(I-M)，所以第i个方程的第j项系数为-(1-P/Q)*Mji。

【代码】

#include<set>
#include<cmath>
#include<queue>
#include<vector>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define trav(u,i) for(int i=front[u];i;i=e[i].nxt)
#define FOR(a,b,c) for(int a=(b);a<=(c);a++)
using namespace std;

typedef long long ll;
const int N = 400; 
const int M = 1e5+10;

ll read() {
    char c=getchar();
    ll f=1,x=0;
    while(!isdigit(c)) {
        if(c=='-') f=-1; c=getchar();
    }
    while(isdigit(c))
        x=x*10+c-'0',c=getchar();
    return x*f;
}

struct Edge {
    int v,nxt;
}e[M];
int en=1,front[N];
void adde(int u,int v) 
{
    e[++en]=(Edge){v,front[u]}; front[u]=en;
}

int n,m,du[N]; 
double P,a[N][N];

void gause()
{
    for(int i=1;i<=n;i++) 
    {
        int r=i;
        for(int j=i+1;j<=n;j++)
            if(fabs(a[j][i])>fabs(a[r][i])) r=i;
        for(int j=1;j<=n+1;j++) swap(a[i][j],a[r][j]);
        for(int j=n+1;j>=i;j--)
            for(int k=i+1;k<=n;k++)
                a[k][j]-=a[k][i]/a[i][i]*a[i][j];
    }
    for(int i=n;i;i--) 
    {
        for(int j=i+1;j<=n;j++)
            a[i][n+1]-=a[i][j]*a[j][n+1];
        a[i][n+1]/=a[i][i];
    }
}

int main()
{
    int x,y;
    n=read(),m=read(),x=read(),y=read();
    P=(double)x/y;
    FOR(i,1,m) {
        x=read(),y=read();
        adde(x,y),adde(y,x);
        du[x]++,du[y]++;
    }
    FOR(u,1,n) {
        a[u][u]=1;
        trav(u,i) {
            int v=e[i].v;
            a[u][v]-=(1.0-P)/du[v];
        }
    }
    a[1][n+1]=P;
    gause();
    FOR(i,1,n)
        printf("%.9lf
",a[i][n+1]);
    return 0;
}