[BZOJ]2707: [SDOI2012]走迷宫 期望+高斯消元

Description

Morenan被困在了一个迷宫里。迷宫可以视为N个点M条边的有向图，其中Morenan处于起点S，迷宫的终点设为T。可惜的是，Morenan非常的脑小，他只会从一个点出发随机沿着一条从该点出发的有向边，到达另一个点。这样，Morenan走的步数可能很长，也可能是无限，更可能到不了终点。若到不了终点，则步数视为无穷大。但你必须想方设法求出Morenan所走步数的期望值。

题解：

这道题似乎是目前做的这种题目中最难的一道了……也很难写……几乎是全程%zyf2000，说一下大概的做法：一个显然的状态表示：$f[i]$表示到i点的期望步数。首先跑一次Tarjan，对于同一个连通分量里面的点，它们的$f$是会互相影响的，所以对于这些点需要用高斯消元来解，就是根据$f[y]=\sum_{x可以到y}\frac{f[x]+1}{degree[y]}$（$degree[y]$表示y的入度）这个式子，来列方程求解。对于不在一个连通分量里面的点，若$x$有连向$y$的边，那就直接加到$f[y]$上就好了。因为到了$t$点就停止，所以我们把原来的边反过来，从$t$开始推，$f[t]=0$就方便推了。然后说一下INF的判断方法：1、到不了$t$。2、存在一个可以到达的不包含$t$的连通分量，并且它的出度为0，因为这样就有几率到达不了终点，步数为INF，期望为步数*概率，还是INF。还有为什么我边数开小会TLE？

代码：

#include<bits/stdc++.h>
using namespace std;
#define LL long long
const int Maxn=10010;
const double eps=1e-10;
int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
    return x*f;
}
int n,m,s,t;
int ex[1000010],ey[1000010];
struct Edge{int y,next;}e[1000010],_e[1000010];
int last[Maxn],len=0;
int _last[Maxn],_len=0;
void ins(int x,int y)
{
    int t=++len;
    e[t].y=y;e[t].next=last[x];last[x]=t;
    t=++_len;swap(x,y);
    _e[t].y=y;_e[t].next=_last[x];_last[x]=t;
}
int b[Maxn][110];
int dfn[Maxn],low[Maxn],sta[Maxn],top=0,id=0,bel[Maxn],cnt=0,po[Maxn];
bool in[Maxn];
void Tarjan(int x)
{
    sta[++top]=x;in[x]=true;
    low[x]=dfn[x]=++id;
    for(int i=last[x];i;i=e[i].next)
    {
        int y=e[i].y;
        if(!dfn[y])Tarjan(y),low[x]=min(low[x],low[y]);
        else if(in[y])low[x]=min(low[x],dfn[y]);
    }
    if(low[x]==dfn[x])
    {
        int i;cnt++;po[cnt]=0;
        do
        {
            i=sta[top--];
            in[i]=false;
            b[cnt][++po[cnt]]=i;
            bel[i]=cnt;
        }while(i!=x);
    }
}
double a[110][110],f[Maxn];
void gauss(int n)
{
    for(int i=1;i<=n;i++)
    {
        if(abs(a[i][i])<=eps)
        {
            for(int j=i+1;j<=n;j++)
            if(abs(a[j][i])>eps)
            {
                for(int k=i;k<=n+1;k++)swap(a[j][k],a[i][k]);
                break;
            }
        }
        for(int j=i+1;j<=n;j++)
        if(abs(a[j][i])>eps)
        {
            double t=a[j][i]/a[i][i];
            for(int k=i;k<=n+1;k++)a[j][k]-=t*a[i][k];
        }
    }
    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];
    }
}
bool can1[Maxn],can2[Maxn];
void work(int x)
{
    can1[x]=can2[bel[x]]=true;
    for(int i=last[x];i;i=e[i].next)
    {
        int y=e[i].y;
        if(can1[y])continue;
        work(y);
    }
}
int outd[Maxn],degree1[Maxn],num[Maxn];
queue<int>q;
int main()
{
    n=read();m=read();s=read();t=read();
    for(int i=1;i<=m;i++)
    {
        ex[i]=read(),ey[i]=read();
        ins(ex[i],ey[i]);degree1[ex[i]]++;
    }
    for(int i=1;i<=n;i++)if(!dfn[i])Tarjan(i);
    work(s);if(!can1[t]){puts("INF");return 0;}
    for(int i=1;i<=m;i++)
    if(bel[ex[i]]!=bel[ey[i]])outd[bel[ex[i]]]++;
    for(int i=1;i<=cnt;i++)
    if(bel[t]!=i&&!outd[i]&&can2[i]){puts("INF");return 0;}
    q.push(bel[t]);
    while(!q.empty())
    {
        int o=q.front();q.pop();
        memset(a,0,sizeof(a));
        for(int i=1;i<=po[o];i++)num[b[o][i]]=i;
        for(int j=1;j<=po[o];j++)
        {
            int x=b[o][j];a[j][j]-=degree1[x];
            if(abs(f[x])>eps)a[j][po[o]+1]=-f[x]*degree1[x];
            for(int i=last[x];i;i=e[i].next)
            {
                int y=e[i].y;
                if(bel[y]!=o)continue;
                a[j][num[y]]+=1.0;a[j][po[o]+1]-=1.0;
            }
            if(x==t)
            {
                for(int i=1;i<=po[o]+1;i++)a[j][i]=0.0;
                a[j][j]=1.0;
            }
        }
        gauss(po[o]);
        for(int j=1;j<=po[o];j++)
        {
            int x=b[o][j];
            f[x]=a[j][po[o]+1];
            for(int i=_last[x];i;i=_e[i].next)
            {
                int y=_e[i].y,by=bel[y];
                if(by==o)continue;
                outd[by]--;if(!outd[by])q.push(by);
                f[y]+=(f[x]+1.0)/degree1[y];
            }
        }
    }
    printf("%.3lf",f[s]);
}