[国家集训队作业]城市规划 题解

传送门：洛谷P4841

题意：求$n$个点简单无向连通图数量。

首先我们设$g(n)$表示$n$个点的简单无向图（不要求连通）数量，对于每一对点都可以连或者不连边，所以$g(n)=2^{C_n^2}$

设$f(n)$表示$n$个点简单无向连通图数量，那么枚举$1$号点所在的连通块至少是多大，可得

$$g(n)=\sum _{i=1}^n{C}_{n-1}^{i-1}f(i)g(n-i)$$

考虑怎么求出$f(n)$。将组合数展开

$$2^{C_n^2}=\sum _{i=1}^n\frac{(n-1)!f(i)2^{C_{n-i}^2}}{(i-1)!(n-i)!}$$

移项得到

$$\frac{2^{C_n^2}}{(n-1)!}=\sum _{i=1}^n\frac{f(i)}{(i-1)!}\frac{2^{C_{n-i}^2}}{(n-i)!}$$

如果设$a_i=\{\begin{array}{ll}\frac{f(i)}{(i-1)!},& i\ge 1\\ 0,& i=0\end{array}$，

$b_i=\frac{2^{C_i^2}}{i!}$，

$c_i=\{\begin{array}{ll}\frac{2^{C_i^2}}{(i-1)!},& i\ge 1\\ 0,& i=0\end{array}$，

并设它们的生成函数$A(x)=\sum _{i\ge 0}{a}_i{x}^i$，$B(x)=\sum _{i\ge 0}{b}_i{x}^i$，$C(x)=\sum _{i\ge 0}{c}_i{x}^i$，那么不难发现

$$C(x)=A(x)B(x)$$

现在已知$C(x)$和$B(x)$，对$B(x)$求逆再与$C(x)$卷积即可知道$A(x)$，最终答案就是$[x^n]A(x)$。

#include <cctype>
#include <cstdio>
#include <climits>
#include <algorithm>

template <typename T> inline void read(T& x) {
    int f = 0, c = getchar(); x = 0;
    while (!isdigit(c)) f |= c == '-', c = getchar();
    while (isdigit(c)) x = x * 10 + c - 48, c = getchar();
    if (f) x = -x;
}
template <typename T, typename... Args>
inline void read(T& x, Args&... args) {
    read(x); read(args...); 
}
template <typename T> void write(T x) {
    if (x < 0) x = -x, putchar('-');
    if (x > 9) write(x / 10);
    putchar(x % 10 + 48);
}
template <typename T> inline void writeln(T x) { write(x); puts(""); }
template <typename T> inline bool chkmin(T& x, const T& y) { return y < x ? (x = y, true) : false; }
template <typename T> inline bool chkmax(T& x, const T& y) { return x < y ? (x = y, true) : false; }

typedef long long LL;

const int maxn = 4e5 + 7;
const LL P = 1004535809, G = 3, Gi = 334845270;

inline LL qpow(LL x, LL k) {
    LL s = 1;
    for (; k; x = x * x % P, k >>= 1)
        if (k & 1) s = s * x % P;
    return s;
}

inline void ntt(LL *A, int *r, int lim, int tp) {
    for (int i = 0; i < lim; ++i)
        if (i < r[i]) std::swap(A[i], A[r[i]]);
    for (int mid = 1; mid < lim; mid <<= 1) {
        LL wn = qpow(tp == 1 ? G : Gi, (P - 1) / (mid << 1));
        for (int j = 0; j < lim; j += mid << 1) {
            LL w = 1;
            for (int k = 0; k < mid; ++k, w = w * wn % P) {
                LL x = A[j + k], y = w * A[j + k + mid] % P;
                A[j + k] = (x + y) % P;
                A[j + k + mid] = (x - y + P) % P;
            }
        }
    }
    if (tp == -1) {
        LL inv = qpow(lim, P - 2);
        for (int i = 0; i < lim; ++i)
            A[i] = A[i] * inv % P;
    }
}

struct Inv {
    LL a[maxn], b[maxn], c[maxn];
    int r[maxn];
    void work(int deg) {
        if (deg == 1) { b[0] = qpow(a[0], P - 2); return; }
        work((deg + 1) >> 1);
        int l = 0, lim = 1;
        while (lim < (deg << 1)) lim <<= 1, ++l;
        for (int i = 0; i < lim; ++i)
            r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
        std::copy(a, a + deg, c);
        std::fill(c + deg, c + lim, 0);
        ntt(c, r, lim, 1);
        ntt(b, r, lim, 1);
        for (int i = 0; i < lim; ++i)
            b[i] = (2ll - c[i] * b[i] % P + P) % P * b[i] % P;
        ntt(b, r, lim, -1);
        std::fill(b + deg, b + lim, 0);
    }
};
Inv inv;

LL fac[maxn], ifac[maxn];
LL b[maxn], c[maxn];
int r[maxn];
int n;

int main() {
    read(n);
    fac[0] = 1;
    for (int i = 1; i <= n; ++i)
        fac[i] = fac[i - 1] * i % P;
    ifac[n] = qpow(fac[n], P - 2);
    for (int i = n - 1; i; --i)
        ifac[i] = ifac[i + 1] * (i + 1) % P;
    ifac[0] = 1;
    b[0] = b[1] = b[2] = 1;
    for (int i = 3; i <= n; ++i)
        b[i] = qpow(2, 1ll * i * (i - 1) / 2) * ifac[i] % P;
    c[0] = 0;
    for (int i = 1; i <= n; ++i)
        c[i] = 1ll * i * b[i];
    std::copy(b, b + n + 1, inv.a);
    inv.work(n + 1);
    std::copy(inv.b, inv.b + n + 1, b);
    int lim = 1, l = 0;
    while (lim <= (n << 1)) lim <<= 1, ++l;
    for (int i = 0; i < lim; ++i)
        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
    ntt(b, r, lim, 1);
    ntt(c, r, lim, 1);
    for (int i = 0; i < lim; ++i)
        b[i] = b[i] * c[i] % P;
    ntt(b, r, lim, -1);
    writeln(b[n] * fac[n - 1] % P);
    return 0;
}