[SDOI2012]走迷宫

· · 题解

题目

        点这里看题目。

分析

        首先考虑输出INF。第一个,如果从S走不到T,那么步数肯定是INF。第二个,如果在ST的路上中间出现了走不到T的岔路口,由于随机游走中我们有可能走到这个岔路里面去,走进去却走不到T,所以这样的话也应该输出INF。(下图,从1\rightarrow 2的时候有\frac 1 2的概率会走到1\rightarrow 3中,然后就走不到2了)

        其次考虑60pts的做法。这个其实很容易看出来是一个高斯消元。设E(u)为从u出发走到T的期望步数,转移是:

E(u)=\frac 1 {deg(u)}\left(\sum_{(u,v)\in E}E(v)\right)+1

        其中deg(u)表示u的度数;(u,v)\in EE为边集。
        由于可能存在环形转移,所以我们用高斯消元,O(n^3)过掉前12组数据。
        接着考虑100pts。其实题面里有一个指向做法的明显提示了,就是——

每个强连通分量的大小不超过100。

        于是我们考虑用它来作文章。先用Tarjan对图进行缩点。然后这个图就变成了DAGDAG上,对于每个uE(u),我们就可以按照拓扑序直接做DP。而这也提示我们,我们在之前的做法中之所以选择高斯消元,是为了解决环状转移,而分量与分量之间现在没有了环,我们何必去用高斯消元呢?
        所以我们可以按照强连通分量的拓扑序(缩点之后的拓扑序)倒着做,每次用高斯消元解决一个强连通分量内部的E。对于当前分量S中的一个点u,我们枚举它出发的边。如果这条边指向的v\in S,我们就把它加入到方程未知元里面;否则,v属于的强连通分量的拓扑序一定比S大,E(v)一定已经计算好了,我们直接把E(v)当做常数放到等号右边。之后解这个方程组,我们就可以得到S中每一个点的E了。所以我们就可以用这样“化整为零”的方法将高斯消元的时间降下来。
        设每个强连通分量的大小为s,这样做的时间是O(\sum s_i^3),很显然\sum s_i^3\le (\sum s_i) \times(\max(s))^2\le n\times100^2,所以时间不会超过O(10^4n)

代码

#include <cstdio>
#include <vector>
#include <iostream>
using namespace std;

const double eps = 1e-7;
const int MAXN = 1e4 + 5, MAXM = 1e6 + 5, MAXS = 105;

template<typename _T>
void read( _T &x )
{
    x = 0;char s = getchar();int f = 1;
    while( s > '9' || s < '0' ){if( s == '-' ) f = -1; s = getchar();}
    while( s >= '0' && s <= '9' ){x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar();}
    x *= f;
}

template<typename _T>
void write( _T x )
{
    if( x < 0 ){ putchar( '-' ); x = ( ~ x ) + 1; }
    if( 9 < x ){ write( x / 10 ); }
    putchar( x % 10 + '0' );
}

template<typename _T>
_T MIN( const _T a, const _T b )
{
    return a < b ? a : b;
}

template<typename _T>
_T ABS( const _T a )
{
    return a < 0 ? -a : a;
}

struct edge
{
    int to, nxt;
}Graph[MAXM];

vector<int> SCC[MAXN];

double A[MAXS][MAXS];
double E[MAXS], f[MAXN];
int q[MAXN], sta[MAXN], top;
int pos[MAXN], seq[MAXS];
int deg[MAXN], head[MAXN], DFN[MAXN], LOW[MAXN], bel[MAXN];
int N, M, S, T, cnt, ID, tot, siz;
bool vis[MAXN], arr[MAXS];

bool equal( const double a, const double b = 0 ) { return ABS( a - b ) < eps; }

void addEdge( const int from, const int to )
{
    Graph[++ cnt].to = to, Graph[cnt].nxt = head[from];
    head[from] = cnt, deg[from] ++;
}

bool chk()
{
    int h = 1, t = 0;
    q[++ t] = S;
    while( h <= t )
    {
        int u = q[h ++];
        for( int i = head[u], v ; i ; i = Graph[i].nxt )
            if( ! vis[v = Graph[i].to] ) vis[v] = true, q[++ t] = v;
    }
    bool flag = vis[T];
    for( int i = 0 ; i <= N ; i ++ ) vis[i] = false;
    return flag;
}

void Tarjan( const int u )
{
    DFN[u] = LOW[u] = ++ ID;
    vis[sta[++ top] = u] = true;
    for( int i = head[u], v ; i ; i = Graph[i].nxt )
    {
        v = Graph[i].to;
        if( ! DFN[v] ) Tarjan( v ), LOW[u] = MIN( LOW[u], LOW[v] );
        else if( vis[v] ) LOW[u] = MIN( LOW[u], DFN[v] );
    }
    if( LOW[u] == DFN[u] )
    {
        int v; tot ++;
        do vis[v = sta[top --]] = false, bel[v] = tot, SCC[tot].push_back( v ); while( v ^ u );
    }
}

void Gauss()
{
    double tmp;
    int opt = 1, indx, lst;
    for( int i = 1 ; i <= siz && opt <= siz ; lst = i ++, opt ++ )
    {
        indx = -1;
        for( int j = i ; j <= siz ; j ++ )
            if( ! equal( A[j][opt] ) && ( indx == -1 || ABS( A[indx][opt] ) < ABS( A[j][opt] ) ) )
                indx = j;
        if( ! ( ~ indx ) ) { i --; continue; }
        std :: swap( A[indx], A[i] );
        tmp = A[i][opt];
        for( int j = 1 ; j <= siz + 1 ; j ++ ) A[i][j] /= tmp;
        for( int j = i + 1 ; j <= siz ; j ++ )
            if( ! equal( A[j][opt] ) )
            {
                tmp = A[j][opt];
                for( int k = 1 ; k <= siz + 1 ; k ++ ) A[j][k] -= A[i][k] * tmp;
            }
    }
    for( int i = lst ; i ; i -- )
    {
        indx = -1;
        for( int j = 1 ; j <= siz ; j ++ ) if( equal( A[i][j], 1 ) ) { indx = j; break; }
        if( ! ( ~ indx ) ) continue;
        E[indx] = A[i][siz + 1];
        for( int j = 1 ; j < i ; j ++ ) A[j][siz + 1] -= E[indx] * A[j][indx], A[j][indx] = 0;
    }
}

int main()
{
    read( N ), read( M ), read( S ), read( T );
    for( int i = 1, fr, to ; i <= M ; i ++ ) read( fr ), read( to ), addEdge( fr, to );
    if( ! chk() ) { puts( "INF" ); return 0; }
    //检查连通性,当然也可以Tarjan的DFN检查
    Tarjan( S );
    for( int i = 1, v ; i <= tot ; i ++ )
    //Tarjan的SCC编号正好就是倒着的拓扑序
    {
        siz = 0;
        for( int j = 0 ; j < SCC[i].size() ; j ++ ) 
            arr[i] |= SCC[i][j] == T, seq[pos[SCC[i][j]] = ++ siz] = SCC[i][j];
            //将SCC的点离散到连续正整数上来
        for( int j = 0 ; j < SCC[i].size() ; j ++ )
            for( int k = head[SCC[i][j]] ; k ; k = Graph[k].nxt )
                if( bel[v = Graph[k].to] ^ i ) arr[i] |= arr[bel[v]];
                //用DP检查是否可以走到T
        if( ! arr[i] ) { puts( "INF" ); return 0; }
        for( int i = 1 ; i <= siz ; i ++ )
            for( int j = 1 ; j <= siz + 1 ; j ++ )
                A[i][j] = 0; 
        for( int j = 0, u, p ; j < SCC[i].size() ; j ++ )
        {
            p = pos[u = SCC[i][j]];
            A[p][p] = 1;
            if( u == T ) continue;
            A[p][siz + 1] = 1;
            for( int k = head[u] ; k ; k = Graph[k].nxt )
            {
                if( bel[v = Graph[k].to] ^ i ) A[p][siz + 1] += 1.0 / deg[u] * f[v];
                //如果已知,直接加入到常数项里面
                else A[p][pos[v]] -= 1.0 / deg[u];
            }
        }
        Gauss();
        for( int i = 1 ; i <= siz ; i ++ ) f[seq[i]] = E[i];
    }
    printf( "%.3lf\n", f[S] );
    return 0;
}