[BZOJ 3456]城市规划

danihao123 posted @ 2018年4月19日 19:56 in 题解 with tags bzoj 多项式求逆 NTT 图计数 , 671 阅读
转载请注明出处:http://danihao123.is-programmer.com/

aji多项式求逆毁我青春,,,

设\(f_i\)表示\(i\)个点的有标号无向联通图,考虑所有可能的图(记\(F_i\)为\(i\)个点的有标号无向图,显然\(F_i = 2^{\binom{i}{2}}\))和\(f\)的关系(使用图计数的经典套路:枚举1所在的联通块大小):

\[F_n = \sum_{i = 1}^n \binom{n-1}{i-1} f_i F_{n - i}\]

看起来事卷积?但是这个卷积没有办法直接用FFT/NTT求(当然分离一下项啊,移下项就可以分治NTT力)。

考虑进一步化简柿子。完全展开后会发现右边有个非常碍眼的\((n-1)!\),所以两边除一下:

\[\frac{F_n}{(n-1)!} =\sum_{i = 1}^n \frac{f_i}{(i-1)!}\cdot \frac{F_{n - i}}{(n - i)!}\]

然后这个问题就很毒瘤了:我们要求的答案多项式(除上那个阶乘)和一个已知多项式做卷积,可以得到另一个已知多项式……这样就需要多项式除法了,于是乎多项式逆元派上了用场。

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
typedef long long ll;
const int maxn = 1020010;
const ll ha = 1004535809LL;
const ll bs = 3LL;
ll pow_mod(ll a, ll b) {
  ll ans = 1LL, res = a % ha;
  while(b) {
    if(1LL & b) ans = (ans * res) % ha;
    res = (res * res) % ha;
    b >>= 1;
  }
  return ans;
}
ll inv(ll x) {
  return pow_mod(x,  ha - 2LL);
}
 
int flip(int bi, int x) {
  int ans = 0;
  for(int i = 0; i < bi; i ++) {
    if((1 << i) & x) {
      ans += (1 << (bi - i - 1));
    }
  }
  return ans;
}
void NTT(ll *A, int bi, bool flag = false) {
  int n = 1 << bi;
  for(int i = 0; i < n; i ++) {
    int v = flip(bi, i);
    if(v < i) std::swap(A[v], A[i]);
  }
  for(int L = 1; L < n; L <<= 1) {
    ll xi = pow_mod(3LL, (ha - 1LL) / (ll(L << 1)));
    if(flag) xi = inv(xi);
    for(int i = 0; i < n; i += (L << 1)) {
      ll w = 1LL;
      for(int j = i; j < i + L; j ++) {
        ll v1 = A[j], v2 = A[j + L];
        A[j] = (v1 + (w * v2) % ha) % ha;
        A[j + L] = (v1 - (w * v2) % ha + ha) % ha;
        w = (w * xi) % ha;
      }
    }
  }
}
void poly_mul(ll *A, ll *B, int bi, ll *C) {
  static ll T1[maxn], T2[maxn];
  int n = (1 << bi);
  std::copy(A, A + n, T1);
  std::copy(B, B + n, T2);
  NTT(T1, bi); NTT(T2, bi);
#ifdef LOCAL
  puts("poly_mul :");
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", A[i]);
  }
  puts("");
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", B[i]);
  }
  puts("");
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", T1[i]);
  }
  puts("");
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", T2[i]);
  }
  puts("");
#endif 
  for(int i = 0; i < n; i ++) {
    T1[i] = (T1[i] * T2[i]) % ha;
  }
#ifdef LOCAL
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", T1[i]);
  }
  puts("");
#endif 
  NTT(T1, bi, true);
  ll inv_n = inv(n);
  for(int i = 0; i < n; i ++) {
    T1[i] = (T1[i] * inv_n) % ha;
  }
  std::copy(T1, T1 + n, C);
#ifdef LOCAL
  for(int i = 0; i < (n); i ++) {
    printf("%lld ", C[i]);
  }
  puts("");
#endif 
}
 
void poly_inv(int mod, ll *B, ll *BB) {
  if(mod == 1) {
    BB[0] = inv(B[0]);
  } else {
    poly_inv((mod + 1) >> 1, B, BB);
    int bi = 0, sz = 1;
    while(sz <= ((mod * 2) + 1)) {
      bi ++; sz <<= 1;
    }
    ll inv_sz = inv(sz);
    static ll tmp[maxn];
    std::copy(B, B + mod, tmp);
    std::fill(tmp + mod, tmp + sz, 0LL);
    NTT(tmp, bi); NTT(BB, bi);
    for(int i = 0; i < sz; i ++) {
      tmp[i] = (tmp[i] * BB[i]) % ha;
      tmp[i] = (tmp[i] * (ha - 1LL)) % ha;
      tmp[i] = (tmp[i] + 2LL) % ha;
      tmp[i] = (tmp[i] * BB[i]) % ha;
    }
    NTT(tmp, bi, true);
    for(int i = 0; i < sz; i ++) {
      tmp[i] = (tmp[i] * inv_sz) % ha;
    }
    std::copy(tmp, tmp + mod, BB);
    std::fill(BB + mod, BB + sz, 0LL);
  }
}
 
int main() {
  static ll fac[maxn], A[maxn], B[maxn], BB[maxn];
  int n; scanf("%d", &n);
  int bi = 0, sz = 1;
  while(sz <= n + 1) {
    bi ++; sz <<= 1;
  }
  fac[0] = 1LL;
  for(int i = 1; i <= n; i ++) {
    fac[i] = (fac[i - 1] * (ll(i))) % ha;
  }
  B[0] = 1LL;
  for(int i = 1; i <= n; i ++) {
    B[i] = pow_mod(2LL, (ll(i)) * (ll(i - 1)) / 2LL);
    B[i] = (B[i] * inv(fac[i])) % ha;
  }
  for(int i = 1; i <= n; i ++) {
    A[i] = pow_mod(2LL, (ll(i)) * (ll(i - 1)) / 2LL);
    A[i] = (A[i] * inv(fac[i - 1])) % ha;
  }
  poly_inv(n + 1, B, BB);
  poly_mul(A, BB, bi, A);
  printf("%lld\n", (A[n] * fac[n - 1]) % ha);
  return 0;
}

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter