題目

([SDOI2012]走迷宮) https://ac.nowcoder.com/acm/problem/20575

題目描述

Morenan被困在了一個迷宮裡。迷宮可以視為N個點M條邊的有向圖，其中Morenan處於起點S，迷宮的終點設為T。可惜的是，Morenan非常的腦小，他只會從一個點出發隨機沿著一條從該點出發的有向邊，到達另一個點。這樣，Morenan走的步數可能很長，也可能是無限，更可能到不了終點。若到不了終點，則步數視為無窮大。但你必須想方設法求出Morenan所走步數的期望值。
n<10000,m<1000000，每個SCC大小不超過100
輸入描述:
第1行4個整數，N,M,S,T第[2, M+1]行每行兩個整數o1, o2，表示有一條從o1到o2的邊。
輸出描述:
一個浮點數，保留小數點3位，為步數的期望值。若期望值為無窮大，則輸出"INF"。

輸入

6 6 1 6
1 2
1 3
2 4
3 5
4 6
5 6
輸出
3.000

輸入

9 12 1 9
1 2
2 3
3 1
3 4
3 7
4 5
5 6
6 4
6 7
7 8
8 9
9 7

輸出

9.500

輸入

2 0 1 2

輸出

INF

思路

我們先考慮INF的情況。如果存在一個點X不能到達T，但是S能夠到達X。那麼就是是INF。我們把邊取反。那麼T能夠
到達的點就是能夠到達T的點，再DFS(S)一次。就可以判斷INF。

然後。我們知道期望反推

但是存在一個問題，就是n太大。我們把圖縮點。從T的scc開始。按拓撲反序進行一個一個scc的高斯消元。
對於同一個scc的邊，按上面的公式就可以了，如果\(u->to。scc[u]!=scc[to]\)

F[T]=0, 那麼F[s]就是答案。

坑點：重邊，有自環

#include<bits/stdc++.h>
#define LL long long
using namespace std;

struct Edge {
    int from, to, nxt;
} e[2000005];
int head[100005], cut=0, N=0;
int scc[100005];
vector<int> g[100005];
struct Scc {
    int dfn[100005], low[100005], vis[100005], T=0;
    stack<int> sk;
    void Addedge(int x, int y) {
        e[++cut]= {x, y, head[x]};
        head[x]=cut;
    }
    void Ta(int u) {
        dfn[u]=low[u]=++T;
        sk.push(u);
        vis[u]=1;
        for(int i=head[u]; i; i=e[i].nxt) {
            int to=e[i].to;
            if(!dfn[to]) {
                Ta(to);
                low[u]=min(low[u], low[to]);
            } else if(vis[to]) {
                low[u]=min(low[u], dfn[to]);
            }
        }
        if(low[u]==dfn[u]) {
            N++;
            while(1) {
                int now=sk.top();
                sk.pop();
                vis[now]=0;
                scc[now]=N;
                g[N].push_back(now);
                if(now==u) {
                    break;
                }
            }
        }
    }
} sc;

int d[200005];
vector<int> G[200005], Gf[200005];

const double eps=1e-5;
double a[205][205];
double ans[205];
map<int, int> mp;
void Gauss(int n) {
    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=j;
        if(fabs(a[r][i])<eps) {  //最大為0,無解
            printf("No Solution\n");
            return;
        }
        if(r!=i)
            for(int j=1; j<=n+1; ++j)
                swap(a[i][j],a[r][j]);
        for(int j=1; j<=n; ++j)  //消元
            if(j!=i) {
                double tmp=a[j][i]/a[i][i];
                for(int k=i; k<=n+1; ++k)
                    a[j][k]-=a[i][k]*tmp;
            }
    }
    for(int i=1; i<=n; ++i)
        ans[i]=a[i][n+1]/a[i][i];
}

queue<int> q;
int vis[200005];
vector<int> c;
void dfs(int u){
    vis[u]=1;
    for(auto x: Gf[u]){
        if(!vis[x]){
            dfs(x);
        }
    }
}

int dfsinf(int u, int t){
    if(u==t) return 1;
    if(vis[u]==0) return 0;
    vis[u]++;
    for(int i=head[u]; i; i=e[i].nxt){
        int x=e[i].to;
        if(vis[x]<2){
            if(!dfsinf(x, t)) return 0;
        }
    }
    return 1;
}

double f[100005];
void dfs(int u, int N) {//每個scc構建高斯消元
    int sz=0;
    vis[u]=1;
    c.push_back(u);
    for(int i=head[u]; i; i=e[i].nxt) {
        sz++;
    }
    for(int i=head[u]; i; i=e[i].nxt) {
        int to=e[i].to;
        if(scc[to]==scc[u]) {
            a[mp[u]][mp[to]]+=-1.0/sz;
            if(!vis[to]) {
                dfs(to, N);
            }
        }
    }
    a[mp[u]][mp[u]]+=1;
    a[mp[u]][N+1]=1;
}

int main() {
    int n, m, s, t;
    scanf("%d%d%d%d", &n, &m, &s, &t);
    for(int i=1; i <=m; i++) {
        int x, y;
        scanf("%d%d", &x, &y);
        sc.Addedge(x, y);
        Gf[y].push_back(x);
    }
    for(int i=1; i<=n; i++) {
        if(!sc.dfn[i]) {
            sc.Ta(i);
        }
    }
    for(int i=1; i<=m; i++) {
        int x=scc[e[i].from], y=scc[e[i].to];
        if(x!=y) {
            d[x]++;
            G[y].push_back(x);
        }
    }

    dfs(t);
    if(!dfsinf(s, t)) {
        printf("INF\n");
        return 0;
    }
    memset(vis, 0, sizeof(vis));

    q.push(scc[t]);
    while(!q.empty()) {//反拓撲序dp
        int now=q.front();
        q.pop();
        mp.clear();
        int N=0;
        for(auto x: g[now]) {
            mp[x]=++N;
        }
        for(auto x: c){
            vis[x]=0;
        }
        c.clear();
        memset(a, 0, sizeof(a));
        dfs(g[now][0], N);

        for(auto x: g[now]) {
            int sz=0;
            for(int i=head[x]; i; i=e[i].nxt) {
                sz++;
            }
            for(int i=head[x]; i; i=e[i].nxt) {
                int y=e[i].to;
                if(scc[x]!=scc[y]) {
                    a[mp[x]][N+1]+=(f[y])/sz;
                }
            }
        }
        if(now==scc[t]) {//f[t]=0
            for(int i=1; i<=N+1; i++) {
                a[mp[t]][i]=0;
            }
            a[mp[t]][mp[t]]=1;
        }

        Gauss(N);
        for(auto x: g[now]) {
            f[x]=ans[mp[x]];
        }

        for(auto x: G[now]) {
            d[x]--;
            if(d[x]==0) {
                q.push(x);
            }
        }
    }
    printf("%.3f\n", fabs(f[s]));

    return 0;
}
/*
8 12 2 4
3 6
4 1
6 8
3 1
2 6
5 4
5 2
8 1
7 3
1 4
2 6
8 7
*/