小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_{k\in S_x}(f[k][m+1]+1)+\frac{1-p_x}{|S_x|}\sum_{k\in 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}\binom{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})\times (\binom{n}{2}-m)}$$

显然, $A_{n,m}=\binom{\binom{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}$ 的时候还要乘上 $\binom{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;
}