小C的岛屿 题解

IOI2018 集训队作业 by 任轩笛

statement

小C有(n)个岛屿。她为了管理岛上的人，决定让这些岛屿连通。

每个岛屿i有一个概率值(p_i)，和一个友好列表(A_i)。

小C首先来到了1号岛屿，她会依次执行下面操作：

1. 设现在在岛屿(x)，有(p_x)的概率产生一条无向边连接两个随机岛屿，这两个岛屿不会相同，也不会已经有边相连。（即在尚不存在的无向边中随机一条加入图中，不会加自环）

2. 如果此时所有岛屿连通，她就会心满意足地离开。

3. 否则她会在当前点的友好列表中，等概率选择一个，走到那座岛屿上。她的不满意度会+1，并重复第1步。

求小C的期望不满意度，答案模1000000007。

tutorial

设 (f[x][m]) 表示现在在 (x) 节点, 已经加入了 (m) 条边, 图还未联通, 不满意度的期望; 显然 (f[1][0]) 就是答案.

设 (g[m]) 表示当前有 (m) 条边且图没有联通, 并且加入原图中不存在的一条边后联通的概率.

[f[x][m]=frac{(1-g[m])p_x}{|S_x|}sum_{kin S_x}(f[k][m+1]+1)+frac{1-p_x}{|S_x|}sum_{kin S_x}(f[k][m]+1) ]

发现这个式子可以分层高斯消元, 复杂度 (O(N^5)).

现在看如何求 (g[m]), 设 (A_{n,m}) 表示 (n) 个点 (m) 条边的无向图方案数, (B_{n,m}) 表示 (n) 个点 (m) 条边联通无向图方案数:

[g[m]=frac{sum_{i=1}^{n-1}inom{n-1}{i-1}sum_{j=0}^{m}B_{i,j}B_{n-i,m-j}i(n-i)}{(A_{n,m}-B_{n,m}) imes (inom{n}{2}-m)} ]

显然, (A_{n,m}=inom{inom{n}{2}}{m}). 现在考虑如何求 (B_{n,m}): 这是一个经典问题, 可以做到 (O(N^4)). 但是由于这道题要对所有的 ((n,m)) 都要求, 所以还是 (O(N^5)) 的. 但是问题不大

考虑将整个图分成若干子集, 使得不同子集间不能有边, 同一子集中不做要求. 设 (h(k)) 表示将整个图分成 (k) 个子集的方案数, 那么含有 (x) 个联通块的方案会被恰好算到 (S(x,k)) 次, 其中 (S) 是第二类斯特林数.

考虑到 (sum_{i=1}^{n}S(n,i)(-1)^{i-1}(i-1)!=[n=1]). 那么我们将所有的 (h(k)) 乘上 ((-1)^{k-1}(k-1)!) 再求和, 得到的就是只有一个联通块的方案数.

于是考虑设 (h[i][j]) 表示用了 (i) 个点, (j) 条边的方案数. 由于需要凑出 ((-1)^{k-1}(k-1)!) 的系数, 可以在转移的时候枚举任意一个点的联通块(而不是1号点所在的联通块)进行转移, 这样就可以带上 (k!) 的系数, 最后再用一次枚举 1 号点联通块的转移即可. 最终求 (B_{n,m}) 的时候还要乘上 (inom{j}{m}) 的系数, 因为转移的时候对于每个子集内部是作完全图来考虑的.

代码:

#include <bits/stdc++.h>
#define N 55
#define M 2555
using namespace std;
const int mod = 1e9 + 7;

inline int fsp(int x, int k = mod - 2) {
    int res = 1;
    while (k) {
        if (k & 1) res = 1ll * res * x % mod;
        x = 1ll * x * x % mod, k >>= 1;
    }
    return res;
}

inline void inc(int &x, long long y) {
    x = (x + y) % mod;
}

int n, m, p[N], a[N][N], B[N][M], C[M][M];
int f[N][M], g[M], h[N][M];
vector<int> S[N];

inline int A(int x, int y) {
    return C[C[x][2]][y];
}

void Gauss() {
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n + 1; ++j)
            a[i][j] = (a[i][j] % mod + mod) % mod;
    for (int i = 1; i <= n; ++i) {
        int inv = fsp(a[i][i]);
        for (int j = i + 1; j <= n; ++j) {
            int tmp = 1ll * inv * a[j][i] % mod;
            for (int k = i; k <= n + 1; ++k)
                inc(a[j][k], mod - 1ll * a[i][k] * tmp % mod);
        }
    }
    for (int i = n; i >= 1; --i) {
        for (int j = i + 1; j <= n; ++j)
            inc(a[i][n + 1], mod - 1ll * a[j][j] * a[i][j] % mod);
        a[i][i] = 1ll * a[i][n + 1] * fsp(a[i][i]) % mod;
    }
}

int main() {
    cin >> n, m = n * (n - 1) / 2;
    for (int i = 1; i <= n; ++i)
        scanf("%d", p + i);
    for (int i = 1; i <= n; ++i) {
        int k; scanf("%d", &k);
        S[i].resize(k);
        for (int j = 0; j < k; ++j)
            scanf("%d", &S[i][j]);
    }
    for (int i = 0; i <= m; ++i)
        for (int j = C[i][0] = 1; j <= i; ++j)
            C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % mod;
    h[0][0] = 1;
    for (int i = 1; i <= n; ++i)
        for (int j = 0; j <= C[i][2]; ++j)
            for (int k = 1; k <= i && C[k][2] <= j; ++k)
                inc(h[i][j], mod - 1ll * h[i - k][j - C[k][2]] * C[i][k] % mod);
    for (int i = n; i >= 1; --i)
        for (int j = C[i][2]; j >= 0; --j) {
            h[i][j] = 0;
            for (int k = 1; k <= i && C[k][2] <= j; ++k)
                inc(h[i][j], 1ll * h[i - k][j - C[k][2]] * C[i - 1][k - 1]);
        }
    for (int i = 1; i <= n; ++i)
        for (int j = 0; j <= m; ++j)
            for (int k = j; k <= m; ++k)
                inc(B[i][j], 1ll * h[i][k] * C[k][j]);
    for (int i = 0, u, v; i <= m; ++i) {
        u = 0, v = 1ll * (A(n, i) - B[n][i] + mod) % mod * (m - i) % mod;
        for (int x = 1; x < n; ++x)
            for (int e = 0; e <= i; ++e)
                inc(u, 1ll * C[n - 1][x - 1] * B[x][e] % mod * B[n - x][i - e] % mod * x * (n - x));
        g[i] = 1ll * u * fsp(v) % mod;
    }
    for (int j = m; j >= 0; --j) {
        memset(a, 0, sizeof a);
        for (int i = 1; i <= n; ++i) {
            int inv = fsp(S[i].size()), x;
            x = 1ll * (mod + 1 - g[j]) * p[i] % mod * inv % mod;
            a[i][i] = 1;
            a[i][n + 1] = mod + 1 - p[i];
            for (int k : S[i]) {
                a[i][k] = 1ll * (p[i] - 1) * inv % mod;
                inc(a[i][n + 1], 1ll * x * (f[k][j + 1] + 1));
            }
        }

        Gauss();
        for (int k = 1; k <= n; ++k)
            f[k][j] = a[k][k];
    }
    return cout << f[1][0] << endl, 0;
}