[UOJ 50][UR #3]链式反应

danihao123 posted @ 2018年9月02日 10:05 in 题解 with tags FFT uoj NTT 多项式牛顿迭代 , 47 阅读
转载请注明出处:http://danihao123.is-programmer.com/

给你一个集合\(S\)(保证其中元素都为小于\(n\)的自然数),定义一棵合法的树为一棵编号满足堆的性质,且非叶子节点都一定有两个可能非叶子的儿子,同时有\(c(c\in S)\)个一定为叶子的儿子。对于所有\(i = 1\ldots n\),求有多少大小为\(i\)的形态不同的合法的树。

\(n\leq 200000\)。

毒瘤,毒瘤.jpg

先考虑DP,定义状态\(f(i)\)表示有\(i\)个点的合法树的数量,那么有:

\[
\begin{aligned}
f(n) &= \frac{1}{2}\sum_{j}\sum_{k}[(n - 1 - j - k)\in S]\binom{n - 1}{j}\binom{n - 1 - j}{k}f(j)f(k)\\
&= \frac{1}{2}\sum_{j}\sum_{k}[(n - 1 - j - k)\in S]\frac{(n - 1)!}{(n - 1 - j - k)!}\frac{f(j)}{j!}\frac{f(k)}{k!}
\end{aligned}
\]

那么定义\(g(n) = \sum_{i + j = n} \frac{f(i)f(j)}{i!j!}\),那么我们可以直接枚举\(j + k\),所以转移方程变成了:

\[
\begin{aligned}
f(n) &= \frac{1}{2}\sum_{s}[(n - 1 - s)\in S]\frac{(n - 1)!}{(n - 1 - s)!}g(s)
\end{aligned}
\]

然后就能\(O(n^2)\)做啦!虽然会T飞……

然后这里面出来了阶乘和组合数,考虑上指数型生成函数(下称EGF)!定义\(F(x)\)为\(f(i)\)的EGF,那么\(g(i)\)的EGF就是\(F^2(x)\),我们再定义\(C(x) = \sum_{i = 0}^{+\infty} [i\in S]\frac{1}{i!}\),然后我们得到了这样一个关系:

\[\frac{1}{2}C(x)F^2(x) + 1 = \frac{\mathrm{d}F(x)}{\mathrm{d}x}\]

然后你会发现这个微分方程搞不动了……(就是怎么查解都查不出来

那么上多项式牛迭?我们还是考虑从膜低次的情况(假设为\(F_0(x)\))推到高次情况(假设为\(F(x)\),此时模式为\(x^p\)),顺便设\(h(t) = \frac{1}{2}C(x)t + 1\),那么很显然有\(\frac{\mathrm{d}F(x)}{\mathrm{d}x}\equiv h(F(x))\pmod{x^p}\),那这怎么搞?

我们还是使用传统艺能,将\(h\)在\(F_0(x)\)处泰勒展开,那么得到:

\[
\begin{aligned}
h(F(x))&\equiv h(F_0(x)) + \frac{h'(F_0(x))}{1!}(F(x) - F_0(x))^1 +\ldots\pmod{x^p}\\
&\equiv h(F_0(x)) + h'(F_0(x))F(x) - h'(F_0(x))F_0(x)\pmod{x^p}\\
\frac{\mathrm{d}F(x)}{\mathrm{d}x}&\equiv h(F_0(x)) + h'(F_0(x))F(x) - h'(F_0(x))F_0(x)\pmod{x^p}
\end{aligned}
\]

然后这泥萌还是没法做啊……

我们再挣扎一下,设\(v(x) = e^{-\int h'(F_0(x))\mathrm{d}x}\),那么我们将上面的等式两边乘\(v(x)\),得:

\[
\begin{aligned}
\frac{\mathrm{d}F(x)}{\mathrm{d}x}v(x)&\equiv v(x)(h(F_0(x)) + h'(F_0(x))F(x) - h'(F_0(x))F_0(x))\pmod{x^p}\\
\frac{\mathrm{d}F(x)}{\mathrm{d}x}v(x)-h'(F_0(x))F(x)v(x)&\equiv v(x)(h(F_0(x)) - h'(F_0(x))F_0(x))\pmod{x^p}\\
\frac{\mathrm{d}F(x)}{\mathrm{d}x}v(x)+F(x)\frac{\mathrm{d}v(x)}{\mathrm{d}x}&\equiv v(x)(h(F_0(x)) - h'(F_0(x))F_0(x))\pmod{x^p}\\
\frac{\mathrm{d}(F(x)v(x))}{\mathrm{d}x}&\equiv v(x)(h(F_0(x)) - h'(F_0(x))F_0(x))\pmod{x^p}\\
F(x)v(x)&\equiv \int v(x)(h(F_0(x)) - h'(F_0(x))F_0(x))\mathrm{d}x\pmod{x^p}\\
F(x)&\equiv \frac{1}{v(x)}\int v(x)(h(F_0(x)) - h'(F_0(x))F_0(x))\mathrm{d}x\pmod{x^p}
\end{aligned}
\]

然后这TM就可以倍增了……至于多项式exp、积分、求逆都是传统艺能了……顺便注意下常数……

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cassert>
#include <algorithm>
#include <functional>
#include <utility>
using ll = long long;
constexpr ll ha = 998244353LL;
inline ll pow_mod(ll a, ll b) {
  ll ans = 1, res = a;
  while(b) {
    if(1LL & b) ans = ans * res % ha;
    res = res * res % ha; b >>= 1;
  }
  return ans;
}
inline ll inv(ll x) {
  return pow_mod(x, ha - 2LL);
}
 
const int maxn = 524300;
ll fac[maxn], ifac[maxn], invp[maxn];
inline void process() {
  invp[1] = 1;
  for(int i = 2; i < maxn; i ++) {
    invp[i] = (ha - (ha / (ll)i));
    invp[i] = invp[i] * invp[ha % (ll)i] % ha;
  }
  fac[0] = 1;
  for(int i = 1; i < maxn; i ++) {
    fac[i] = fac[i - 1] * (ll)i % ha;
  }
  ifac[maxn - 1] = inv(fac[maxn - 1]);
  for(int i = maxn - 2; i >= 0; i --) {
    ifac[i] = (ll(i + 1)) * ifac[i + 1] % ha;
  }
}
 
inline void NTT(ll *A, int bi, bool flag = false) {
  int n = 1 << bi;
  for(int i = 1, j = 0; i < n; i ++) {
    int k = n;
    do j ^= (k >>= 1); while ((j & k) == 0);
    if(i < j) std::swap(A[j], A[i]);
  }
  for(int L = 1; L < n; L <<= 1) {
    const bool bb = (L < 4);
    ll c = (ha - 1LL) / (ll(L << 1));
    if(flag) c = ha - 1LL - c;
    ll xi_n = pow_mod(3LL, c);
    ll *B = A + L;
    for(int i = 0; i < n; i += (L << 1)) {
      ll xi = 1;
      if(true) {
        for(int j = i; j < i + L; j ++) {
          ll x = A[j], y = B[j] * xi % ha;
          A[j] = (x + y) % ha;
          B[j] = (x - y + ha) % ha;
          xi = xi * xi_n % ha;
        }
      } else {
        for(int j = i; j < i + L; j += 4) {
          ll x = A[j], y = B[j] * xi % ha;
          A[j] = (x + y); if(A[j] >= ha) A[j] -= ha;
          B[j] = (x - y + ha); if(B[j] >= ha) B[j] -= ha;
          xi = xi * xi_n % ha;
          x = A[j + 1], y = B[j + 1] * xi % ha;
          A[j + 1] = (x + y); if(A[j + 1] >= ha) A[j + 1] -= ha;
          B[j + 1] = (x - y + ha); if(B[j + 1] >= ha) B[j + 1] -= ha;
          xi = xi * xi_n % ha;
          x = A[j + 2], y = B[j + 2] * xi % ha;
          A[j + 2] = (x + y); if(A[j + 2] >= ha) A[j + 2] -= ha;
          B[j + 2] = (x - y + ha); if(B[j + 2] >= ha) B[j + 2] -= ha;
          xi = xi * xi_n % ha;
          x = A[j + 3], y = B[j + 3] * xi % ha;
          A[j + 3] = (x + y); if(A[j + 3] >= ha) A[j + 3] -= ha;
          B[j + 3] = (x - y + ha); if(B[j + 3] >= ha) B[j + 3] -= ha;
          xi = xi * xi_n % ha;
        }
      }
    }
  }
  if(flag) {
    ll inv_n = invp[n];
    for(int i = 0; i < n; i ++) {
      A[i] = A[i] * inv_n % ha;
    }
  }
}
 
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;
    }
    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);
    std::copy(tmp, tmp + mod, BB);
    std::fill(BB + mod, BB + sz, 0LL);
  }
}
inline void poly_ln(int mod, ll *B, ll *BB) {
  int bi = 0, sz = 1;
  while(sz <= ((mod * 2) - 1)) {
    bi ++; sz <<= 1;
  }
  static ll tmp[maxn]; memset(tmp, 0, sizeof(ll) * sz);
  poly_inv(mod, B, BB);
  for(int i = 0; i < mod - 1; i ++) tmp[i] = B[i + 1] * (ll(i + 1)) % ha;
  
  NTT(tmp, bi); NTT(BB, bi);
  for(int i = 0; i < sz; i ++) BB[i] = BB[i] * tmp[i] % ha;
  NTT(BB, bi, true);
  for(int i = mod - 1; i > 0; i --) BB[i] = BB[i - 1] * invp[i] % ha;
  BB[0] = 0; std::fill(BB + mod, BB + sz, 0LL);
  // free(tmp);
}
void poly_exp(int mod, ll *B, ll *BB) {
  if(mod == 1) {
    BB[0] = 1;
  } else {
    poly_exp((mod + 1) >> 1, B, BB);
    int bi = 0, sz = 1;
    while(sz <= ((mod * 2) - 1)) {
      bi ++; sz <<= 1;
    }
     
    std::fill(BB + mod, BB + sz, 0LL);
    static ll tmp[maxn]; memset(tmp, 0, sizeof(ll) * sz);
    poly_ln(mod, BB, tmp);
    for(int i = 0; i < mod; i ++) {
      tmp[i] = (B[i] - tmp[i] + ha);
      if(tmp[i] >= ha) tmp[i] -= ha;
    }
    tmp[0] ++;
    NTT(tmp, bi); NTT(BB, bi);
    for(int i = 0; i < sz; i ++) BB[i] = tmp[i] * BB[i] % ha;
    NTT(BB, bi, true);
     
    std::fill(BB + mod, BB + sz, 0LL);
    // free(tmp);
  }
}
void poly_equation(int mod, ll *B, ll *BB) {
  if(mod == 1) {
    BB[0] = 0;
  } else {
    poly_equation((mod + 1) >> 1, B, BB);
    int bi = 0, sz = 1;
    while(sz <= (mod * 2 - 1)) {
      bi ++; sz <<= 1;
    }

    static ll t1[maxn], rv[maxn], v[maxn];
    std::copy(B, B + mod, t1); memset(t1 + mod, 0, (sz - mod) * sizeof(ll));
    NTT(t1, bi);
    // std::fill(BB + mod, BB + sz, 0LL);
    NTT(BB, bi);
    for(int i = 0; i < sz; i ++) rv[i] = (t1[i] * BB[i] % ha);
    NTT(rv, bi, true);
    for(int i = mod - 1; i > 0; i --) rv[i] = (ha - rv[i - 1] * invp[i] % ha);
    rv[0] = 0;
    poly_exp(mod, rv, v);
    
    static ll rn[maxn]; memset(rn, 0, sizeof(ll) * sz);
    // for(int i = 0; i < sz; i ++) rn[i] = BB[i] * BB[i] % ha;
    // NTT(rn, bi, true); std::fill(rn + mod, rn + sz, 0LL);
    // NTT(rn, bi);
    static const ll ii_2 = ha - inv(2);
    for(int i = 0; i < sz; i ++) {
      rn[i] = (1LL + (ii_2 * (BB[i] * BB[i] % ha) % ha) * t1[i]) % ha;
    }
    NTT(rn, bi, true); std::fill(rn + mod, rn + sz, 0LL);
    memset(rv, 0, sizeof(ll) * sz);
    poly_inv(mod, v, rv); NTT(v, bi); NTT(rn, bi);
    for(int i = 0; i < sz; i ++) rn[i] = rn[i] * v[i] % ha;
    NTT(rn, bi, true); std::fill(rn + mod, rn + sz, 0LL);
    for(int i = mod - 1; i > 0; i --) rn[i] = rn[i - 1] * invp[i] % ha;
    rn[0] = 0;
    NTT(rn, bi); NTT(rv, bi);
    for(int i = 0; i < sz; i ++) rn[i] = rn[i] * rv[i] % ha;
    NTT(rn, bi, true);
    std::copy(rn, rn + mod, BB); std::fill(BB + mod, BB + sz, 0LL);
    // free(rn); free(v); free(rv); free(t1);
  }
}

char S[maxn];
ll A[maxn], C[maxn];
int main() {
#ifdef LOCAL
  freopen("gou.in", "r", stdin);
  freopen("outp", "w", stdout);
  srand(238238);
  for(int i = 0; i < 200000; i ++) {
    if(rand() & 1) S[i] = '1';
  }
#endif
  process();
  int n; scanf("%d%s", &n, S);
  for(int i = 0; i < n; i ++) {
    if(S[i] == '1') {
      C[i] = ifac[i];
    }
  }
  poly_equation(n + 1, C, A);
  for(int i = 1; i <= n; i ++) {
    printf("%lld\n", A[i] * fac[i] % ha);
  }
  return 0;
}

登录 *


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