12/26
23:29
OI

烷基计数 [Burnside引理][生成函数][FFT]

烷基计数 加强版 加强版

因为我太菜没有约,只能去学高中数学/化学(

题意就是求无标号,有根,儿子不超过3个的树的数量。

做法还是很神奇的。

做一个简单DP:令f_x表示x个点的答案:

f_x = \sum_{i+j+k+1 = x} f_i f_j f_k

\uparrow明显是错的。要去除重复的答案。

其实暴力容斥就行,练一下数数,用Burnside引理。方案数是不动点的平均值。

推一推可以找出六种置换中不动点分别有多少。

f_x = \frac{ \sum_{i+j+k+1 = x} f_i f_j f_k + 3\sum_{2i + j + 1 = x} f_i f_j + 2 \sum_{3i + 1 = x} f_i }{6}

写成生成函数:

F(x) = x \frac{ F(x)^3 + 3 F(x^2)F(x) + 2F(x^3) }{6} + 1

考虑牛顿迭代。在\mod \lfloor \frac{x}{2} \rfloor意义下求出F_0(x)后,实际上当前的x,F(x^2),F(x^3)都知道了。它们可以看做常数项。

A(F(x)) \equiv x \frac{ F(x)^3 + 3 F(x^2)F(x) + 2F(x^3) }{6} + 1 – F(x) \equiv 0

F(x) = F_0(x) – \frac{A(F_0(x))}{A(F_0(x))’}

A(F(x))’ = x \frac{3F(x)^2 + 3F(x^2)}{6} – 1

带进去就可以做了

#include <bits/stdc++.h>

typedef long long ll;

const int N = 8e5 + 233, P = 998244353;

inline ll fpow(ll x, int y) {
  ll ret = 1;
  for ( ; y; y >>= 1, x = x * x % P)
    if (y & 1) ret = ret * x % P;
  return ret;
}

const int Gp = 3, Gpi = fpow(Gp, P - 2);

int rev[N], L, lim;

inline void init(int n) {
  L = 0, lim = 1;
  while (lim < n) lim <<= 1, ++L;
  for (int i = 1; i < lim; ++i)
    rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (L - 1));
}

inline void NTT(ll f[], int o) {
  for (int i = 1; i < lim; ++i)
    if (i < rev[i])
      std::swap(f[i], f[rev[i]]);
  for (int i = 1; i < lim; i <<= 1) {
    ll T = fpow(o == 1 ? Gp : Gpi, (P - 1) / (i << 1));
    for (int j = 0; j < lim; j += i << 1) {
      ll w = 1;
      for (int k = 0; k < i; ++k, w = w * T % P) {
        ll nx = f[j + k], ny = f[i + j + k] * w % P;
        f[j + k] = nx + ny >= P ? nx + ny - P : nx + ny;
        f[i + j + k] = nx - ny >= 0 ? nx - ny : nx - ny + P;
      }
    }
  }
  if (o == -1) {
    ll inv = fpow(lim, P - 2);
    for (int i = 0; i < lim; ++i)
      f[i] = f[i] * inv % P;
  }
}

void get_inv(ll F[], ll G[], int n) {
  if (n == 1) {
    G[0] = fpow(F[0], P - 2);
    return;
  }
  get_inv(F, G, (n + 1) >> 1);
  static ll T[N];
  init(n << 1);
  for (int i = 0; i < lim; ++i)
    T[i] = i < n ? F[i] : 0;
  NTT(G, 1), NTT(T, 1);
  for (int i = 0; i < lim; ++i)
    G[i] = ((2 - G[i] * T[i]) % P + P) * G[i] % P;
  NTT(G, -1);
  for (int i = n; i < lim; ++i)
    G[i] = 0;
}

int n; ll A[N], B[N], C[N], D[N], E[N], F[N];

void solve(int n) {
  if (n == 1) return A[0] = 1, void();
  solve((n + 1) >> 1); init(n << 1);
  for (int i = 0; i < lim; ++i)
    B[i] = C[i] = D[i] = E[i] = F[i] = 0;
  for (int i = 0; i < n; ++i) {
    if (i * 2 < n) B[i * 2] = A[i];
    if (i * 3 < n) C[i * 3] = A[i];
  }
  NTT(A, 1), NTT(B, 1), NTT(C, 1);
  for (int i = 0; i < lim; ++i) {
    D[i] = (A[i] * A[i] % P * A[i] % P + 3 * B[i] * A[i] % P + 2 * C[i]) % P;
    E[i] = (3 * A[i] * A[i] % P + 3 * B[i]) % P;
  }
  NTT(A, -1), NTT(D, -1), NTT(E, -1);
  for (int i = lim - 1; i; --i)
    D[i] = D[i - 1], E[i] = E[i - 1];
  D[0] = 6; E[0] = P - 6;
  for (int i = 0; i < lim; ++i)
    D[i] = (D[i] - 6 * A[i] % P + P) % P;
  get_inv(E, F, lim);
  NTT(D, 1), NTT(F, 1);
  for (int i = 0; i < lim; ++i)
    F[i] = (P - D[i] * F[i] % P) % P;
  NTT(F, -1);
  for (int i = 0; i < lim; ++i)
    A[i] = (A[i] + F[i]) % P;
  for (int i = n; i < lim; ++i)
    A[i] = 0;
}

int main() {
  std::cin >> n;
  init(n);
  solve(n + 1);
  std::cout << A[n] << std::endl;
  return 0;
}