[生成函数][NTT][多项式求逆]BZOJ 3456: 城市规划

Description

求n个有标号点的联通图的方案数。

Solution

设fn为n个有标号点的联通图的方案数。
考虑容斥。n个有标号点的一般图的方案数为2(i2)。
考虑图中的一个点所在联通块大小，设其为i。那么就有(n−1i−1)种选法，剩下的乱选，为2(n−i2)。
所以就有了这样的递推式：

fn=2(i2)−∑i=1n−1(n−1i−1)2(n−i2)fi

移项整理就得到了：

2(n2)(n−1)!=∑i=1nfi(i−1)!2(n−12)(n−i)!

考虑生成函数

F(x)G(x)H(x)FG=====∑n≥02(n2)(n−1)!xn∑n≥0fn(n−1)!xn∑n≥02(n2)n!xnG∗HF∗H−1

多项式求逆就好啦。。
好像现在才知道FFT NTT都是循环卷积的，所以要开大一倍。。。

UPD：
后来又推了一发。。
设

G(x)C(x)==∑n≥02(n2)n!xn∑n≥0fnn!xn

就有

G′(x)G(x)==C′(x)⋅G(x)eC(x)

好像就可以牛顿迭代做了呢。。

#include <bits/stdc++.h>using namespace std;const int N = 404040;const int MOD = 1004535809;const int G = 3;typedef long long ll;inline char get(void) {    static char buf[100000], *S = buf, *T = buf;    if (S == T) {        T = (S = buf) + fread(buf, 1, 100000, stdin);        if (S == T) return EOF;    }    return *S++;}inline void read(int &x) {    static char c; x = 0;    for (c = get(); c < '0' || c > '9'; c = get());    for (; c >= '0' && c <= '9'; c = get()) x = x * 10 + c - '0';}int w[2][N];int g, ig, num;int n, m;int R[N];int fac[N], inv[N];int F[N], H[N], iH[N];inline int Pow(int a, int b) {    int c = 1;    while (b) {        if (b & 1) c = (ll)c * a % MOD;        b >>= 1; a = (ll)a * a % MOD;    }    return c;}inline int Inv(int x) {    return Pow(x, MOD - 2);}inline int Mod(int x) {    while (x >= MOD) x -= MOD; return x;}void Prep(int n) {    g = Pow(G, (MOD - 1) / n);    ig = Inv(g); num = n;    w[0][0] = w[1][0] = 1;    for (int i = 1; i <= n; i++) {        w[0][i] = (ll)w[0][i - 1] * ig % MOD;        w[1][i] = (ll)w[1][i - 1] * g % MOD;    }}inline void FFT(int *a, int n, int r) {    static int x, y, INV;    for (int i = 0; i < n; i++)        if (R[i] > i) swap(a[i], a[R[i]]);    for (int i = 1; i < n; i <<= 1)        for (int j = 0; j < n; j += (i << 1))            for (int k = 0; k < i; k++) {                x = a[j + k]; y = (ll)a[j + k + i] * w[r][num / (i << 1) * k] % MOD;                a[j + k] = Mod(x + y); a[j + k + i] = Mod(x - y + MOD);            }    if (!r) {        INV = Inv(n);        for (int i = 0; i < n; i++)            a[i] = (ll)a[i] * INV % MOD;    }}void GetInv(int *a, int *b, int n) {    static int tmp[N];    if (n == 1) return (void)(b[0] = Inv(a[0]));    GetInv(a, b, n >> 1);    for (int i = 0; i < n; i++) {        tmp[i] = a[i]; tmp[i + n] = 0;    }    int L = 0; while (!(n >> L & 1)) L++;    for (int i = 1; i < (n << 1); i++)        R[i] = (R[i >> 1] >> 1) | ((i & 1) << L);    FFT(tmp, n << 1, 1); FFT(b, n << 1, 1);    for (int i = 0; i < (n << 1); i++)        tmp[i] = (ll)b[i] * (2 + MOD - (ll)tmp[i] * b[i] % MOD) % MOD;    FFT(tmp, n << 1, 0);    for (int i = 0; i < n; i++) {        b[i] = tmp[i]; b[n + i] = 0;    } }inline int Calc(int x) {    return (ll)x * (x - 1) / 2 % (MOD - 1);}int main(void) {    freopen("1.in", "r", stdin);    read(n); inv[1] = 1;    for (m = 1; m <= n; m <<= 1);    for (int i = 2; i < m; i++)        inv[i] = (ll)(MOD - MOD / i) * inv[MOD % i] % MOD;    fac[0] = inv[0] = 1;    for (int i = 1; i < m; i++) {        fac[i] = (ll)fac[i - 1] * i % MOD;        inv[i] = (ll)inv[i - 1] * inv[i] % MOD;    }    for (int i = 0; i < m; i++) {        F[i] = (ll)Pow(2, Calc(i)) * inv[i - 1] % MOD;        H[i] = (ll)Pow(2, Calc(i)) * inv[i] % MOD;    }    Prep(m << 1); GetInv(H, iH, m);    FFT(F, m <<= 1, 1); FFT(iH, m, 1);    for (int i = 0; i < m; i++) F[i] = (ll)F[i] * iH[i] % MOD;    FFT(F, m, 0);    printf("%d\n", (ll)F[n] * fac[n - 1] % MOD);    return 0;}