烷基计数 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$$ 带进去就可以做了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
#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;
}