BZOJ 3456: 城市规划

有标号集合计数
对于 \(A\) 中元素形成的集合的指数型生成函数 EGF 为 $$B(x)=\sum_{i\geq 0}\frac{A(x)^i}{i!}=e$$
现在设 \(G\) 为所有简单无向图，则 \(G\) 的 EGF 为 $$G(x)=\sum_{i\geq 0}2^{{\binom{n}{2}}\frac{x}i}{i!}$$
而 \(C\) 为简单无向连通图，连通图可以看成若干个连通分量组成的集合，也就是 $$G(x)=e^{C(x)}$$
那么 $$C(x) = \ln G(x)$$
多项式取对数即可

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
const int MOD = 1004535809, inv2 = (MOD + 1) / 2;
const int N = 1 << 20;
int inv[N];

inline void M(int &x) { if (x >= MOD) x -= MOD; if (x < 0) x += MOD; }
inline int mul_mod(int x, int y) {return (ll)x * y % MOD;}
inline int qp(int a, int n = MOD - 2) {
  int ret = 1;
  for (a %= MOD; n; n >>= 1) {
    if (n & 1) ret = mul_mod(ret, a);
    a = mul_mod(a, a);
  }
  return ret;
}

namespace FFT {
  #define toll(val) ((long long)((val)+(val<0?-0.5:0.5))%MOD)
  const db PI = acos(-1.0L);
  const int M0 = 1 << 15;
  struct Complex {
    db re, im;
    Complex(db _re = 0, db _im = 0) : re(_re), im(_im) {}
    Complex operator!() const {return Complex(re, -im);}
    Complex operator+(const Complex &b) const {return Complex(re + b.re, im + b.im);}
    Complex operator-(const Complex &b) const {return Complex(re - b.re, im - b.im);}
    Complex operator*(const Complex &b) const {return Complex(re * b.re - im * b.im, re * b.im + im * b.re);}
    Complex& operator+=(const Complex &b) {re += b.re, im += b.im; return *this;}
    Complex& operator*=(const Complex &b) {db x = re * b.re - im * b.im, y = re * b.im + im * b.re; re = x, im = y; return *this;}
    Complex& operator/=(db d) {re /= d, im /= d; return *this;}
  };
  int base = 1;
  int rev[N] = {0, 1};
  Complex roots[N] = {{0, 0}, {1, 0}};
  inline void ensure_base(int nbase) {
    if (base > nbase) return;
    //rev.resize(1 << nbase), roots.resize(1 << nbase);
    for (int i = 0; i < (1 << nbase); i++)
      rev[i] = (rev[i >> 1] >> 1) + ((i & 1) << (nbase - 1));
    for (; base < nbase; ++base) {
      db angle = 2 * PI / (1 << (base + 1));
      for (int i = 1 << (base - 1); i < (1 << base); i++) {
        db angle_i = angle * (2 * i + 1 - (1 << base));
        roots[i << 1] = roots[i];
        roots[(i << 1) + 1] = Complex(cos(angle_i), sin(angle_i));
      }
    }
  }
  inline void _fft(Complex a[], int n, bool isInv = false) {
    int zeros = __builtin_ctz(n);
    ensure_base(zeros);
    int shift = base - zeros;
    for (int i = 0; i < n; ++i)
      if (i < (rev[i] >> shift)) swap(a[i], a[rev[i] >> shift]);
    for (int k = 1; k < n; k <<= 1) {
      for (int i = 0; i < n; i += 2 * k) {
        for (int j = 0; j < k; ++j) {
          Complex z(a[j + i + k] * (isInv ? !roots[j + k] : roots[j + k]));
          a[j + i + k] = a[j + i] - z, a[j + i] += z;
        }
      }
    }
    if (isInv) for (int i = 0; i < n; ++i) a[i] /= n;
  }
  inline int get_n(int x) {
    int n = 0x80000000u >> __builtin_clz(x);
    return n < x ? n << 1 : n;
  }
  Complex pa[N], pb[N], tem[N];
  inline void mul_force(int a[], int len_a, int b[], int len_b, int res[]) {
    static int tempSum[N];
    int *sum = res, tot = len_a + len_b - 1;
    if (res == a || res == b) sum = tempSum;
    memset(sum, 0, sizeof(int)*tot);
    for (int i = 0; i < len_a; ++i) {
      for (int j = 0; j < len_b; ++j) {
        sum[i + j] += mul_mod(a[i], b[j]);
        if (sum[i + j] >= MOD) sum[i + j] -= MOD;
      }
    }
    if (sum == tempSum) memcpy(res, sum, sizeof(int)*tot);
  }
  inline void convo(int a[], int len_a, int b[], int len_b, int res[]) {
    if (ll(len_a) * len_b <= 1 << 18) {
      mul_force(a, len_a, b, len_b, res); return;
    }
    int tot = len_a + len_b - 1, n = get_n(tot);
    for (int i = 0; i < len_a; ++i)
      pa[i] = Complex(a[i] / M0, a[i] % M0);
    memset(pa + len_a, 0, sizeof(Complex) * (n - len_a));
    _fft(pa, n);
    if (a == b && len_a == len_b) {
      memcpy(pb, pa, sizeof(Complex)*n);
    } else {
      for (int i = 0; i < len_b; ++i)
        pb[i] = Complex(b[i] / M0, b[i] % M0);
      memset(pb + len_b, 0, sizeof(Complex) * (n - len_b));
      _fft(pb, n);
    }
    tem[0] = Complex(pa[0].im * pb[0].im, 0);
    for (int i = 1; i < n; ++i)
      (tem[i] = (pa[i] - !pa[n - i]) * (pb[i] - !pb[n - i])) /= -4; // tem = Fra*Frb
    for (int i = 0; i < n; ++i)
      (pa[i] *= pb[i]) += tem[i];
    _fft(pa, n, true), _fft(tem, n, true);// IDFT(Fka*Fkb + i(Fra*Fkb + Fka*Frb))、IDFT(Fra*Frb)
    for (int i = 0; i < tot; ++i)
      res[i] = (toll(pa[i].re) * M0 * M0 + toll(pa[i].im) * M0 + toll(tem[i].re)) % MOD;
  }
  void poly_inv(int a[], int b[], int n) {
    if (n == 1) {b[0] = qp(a[0], MOD - 2); b[1] = 0; return;}
    static int res[N];
    int hN = (n + 1) >> 1;
    poly_inv(a, b, hN);
    convo(b, hN, b, hN, res);
    convo(a, n, res, hN * 2 - 1, res);
    for (int i = hN; i < n; ++i) b[i] = res[i] == 0 ? 0 : MOD - res[i];
  }
  void poly_sqr(int *a, int *b, int n) {
    static int A[N], B[N];
    b[0] = 1; b[1] = 0;
    int len, lim;
    for (len = 1; len < (n << 1); len <<= 1) {
      lim = len << 1;
      for (int i = 0; i < len; i++) A[i] = a[i];
      poly_inv(b, B, lim >> 1);
      convo(A, len, B, len, A);
      for (int i = 0; i < len; i++) b[i] = 1ll * (b[i] + A[i]) * inv2 % MOD;
      for (int i = len; i < lim; i++) b[i] = 0;
    }
    for (int i = 0; i < len; i++) A[i] = B[i] = 0;
    for (int i = n; i < len; i++) b[i] = 0;
  }
  void Int(int *a, int *b, int n) {
    for (int i = n; i >= 1; i--)
      b[i] = 1ll * a[i - 1] * inv[i] % MOD;
    b[0] = 0;
  }
  void Diff(int *a, int *b, int n) {
    for (int i = 0; i < n - 1; i++)
      b[i] = 1ll * a[i + 1] * (i + 1) % MOD;
    b[n - 1] = 0;
  }
  void poly_ln(int *a, int *b, int n) {
    static int c[N], d[N];
    poly_inv(a, c, n);
    Diff(a, d, n);
    convo(d, n - 1, c, n, d);
    Int(d, b, n);
  }
  void poly_exp(int *a, int *b, int n) {
    static int c[N], d[N];
    d[0] = 1;
    for (int len = 1; len < (n << 1); len <<= 1) {
      //memset(b + len, 0, sizeof(int) * len);
      //memset(c, 0, sizeof(int) * len);
      poly_ln(d, c, len);
      for (int i = 0; i < len; i++)
        M(c[i] = (i == 0) + a[i] - c[i]);
      convo(d, len, c, len, d);
    }
    for (int i = 0; i < n; i++)
      b[i] = d[i];
  }
  void poly_pow(int *a, int *b, int n, int k) {
    static int c[N], d[N];
    poly_ln(a, c, n);
    for (int i = 0; i < n; i++)
      d[i] = 1LL * c[i] * k % MOD;
    poly_exp(d, b, n);
  }
  void poly_Pow(int *a, int *b, int n, int k) {
    int d = n;
    for (int i = 0; i < n; i++)
      if (a[i]) {
        d = i;
        break;
      }
    static int c[N];
    //memset(b, 0, sizeof(int) * n);
    if (1LL * k * d >= n) return;
    int fi = a[d], invv = qp(fi);
    for (int i = 0; i < n - d; i++)
      b[i] = 1LL * a[i + d] * invv % MOD;
    poly_pow(b, c, n, k);
    for (int i = 0; i < k * d; i++)
      b[i] = 0;
    int p = qp(fi, k);
    for (int i = k * d; i < n; i++)
      b[i] = 1LL * c[i - k * d] * p % MOD;
  }
}

int C(int n) {
  return 1LL * n * (n - 1) / 2 % (MOD - 1);
}

int fac[N], fac_inv[N], a[N], b[N];

int main() {
  int n;
  scanf("%d", &n);
  inv[1] = fac[0] = fac[1] = fac_inv[0] = fac_inv[1] = 1;
  a[0] = a[1] = 1;
  for (int i = 2; i <= 2 * n; i++) {
    inv[i] = 1ll * (MOD - MOD / i) * inv[MOD % i] % MOD;
    fac[i] = 1LL * i * fac[i - 1] % MOD;
    fac_inv[i] = 1LL * fac_inv[i - 1] * inv[i] % MOD;
    a[i] = 1LL * qp(2, C(i)) * fac_inv[i] % MOD;
  }
  FFT::poly_ln(a, b, n + 1);
  printf("%lld\n", 1LL * b[n] * fac[n] % MOD);
  return 0;
}