LOJ

分析

要求的就是每个节点的儿子数都不超过 $3$ 的无标号有根树数。

设其 OGF 为 $F(x)$,考虑列出一个关于 $F(x)$ 的方程。

不妨规定 $[x^0]F(x) = 1$,那么可以看成每个节点的儿子数都恰好为 $3$。

直接得出 $F(x) = 1 + F(x) ^ 3$ 显然是错误的,因为我们同构的情况被重复计算了。

因此可以考虑 Burnside 引理。对于 $(1, 2, 3)$,所有情况都是不动点;对于 $(1, 3, 2)$、$(3, 2, 1)$、$(2, 1, 3)$,在置换中成环的两棵子树同构的情况是不动点;对于 $(3, 1, 2)$、$(2, 3, 1)$,三棵子树都同构的情况是不动点。

这样子我们就可以列出方程

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

考虑使用牛顿迭代解出 $F(x)$。设 $G(F(x)) = 1 + x \frac{F(x) ^ 3 + 3 F(x ^ 2) F(x) + 2 F(x ^ 3)}{6} - F(x) = 0$。

假设当前已知 $F_0(x) = F(x) \bmod x ^ {\frac{n}{2}}$,要求 $F(x) \bmod x ^ n$。那么事实上 $F(x ^ 2) \bmod{x ^ n}$ 和 $F(x ^ 3) \bmod x ^ n$ 已知,因此可以将其看作常量 $A(x)$ 和 $B(x)$。这样子有

$$ G'(F(x)) = x \frac{F(x) ^ 2 + 3 A(x)}{6} - 1 $$

套牛顿迭代的式子

$$ F(x) = F_0(x) - \frac{G(F_0(x))}{G'(F_0(x))} $$

这样子就可以求出 $F(x)$ 了。时间复杂度 $\mathcal{O}(n\log n)$。

代码

// ====================================
//   author: M_sea
//   website: https://m-sea-blog.com/
// ====================================
#include <bits/stdc++.h>
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
#define debug(...) fprintf(stderr, __VA_ARGS__)
using namespace std;
typedef long long ll;

int read() {
    int X = 0, w = 1; char c = getchar();
    while (c < '0' || c > '9') { if (c == '-') w = -1; c = getchar(); }
    while (c >= '0' && c <= '9') X = X * 10 + c - '0', c = getchar();
    return X * w;
}

const int N = 262144 + 10;
const int mod = 998244353, inv2 = 499122177;
int qpow(int a, int b) {
    int c = 1;
    for (; b; b >>= 1, a = 1ll * a * a % mod)
        if (b & 1) c = 1ll * c * a % mod;
    return c;
}

namespace Poly {
    int r[N];
    void NTT(int *A, int n, int op) {
        for (int i = 0; i < n; ++i)
            if (i < r[i]) swap(A[i], A[r[i]]);
        for (int i = 1; i < n; i <<= 1) {
            int rot = qpow(op == 1 ? 3 : 332748118, (mod - 1) / (i << 1));
            for (int j = 0; j < n; j += i << 1)
                for (int k = 0, w = 1; k < i; ++k, w = 1ll * w * rot % mod) {
                    int x = A[j + k], y = 1ll * w * A[j + k + i] % mod;
                    A[j + k] = (x + y) % mod, A[j + k + i] = (x - y + mod) % mod;
                }
        }
        if (op == -1) {
            int inv = qpow(n, mod - 2);
            for (int i = 0; i < n; ++i) A[i] = 1ll * A[i] * inv % mod;
        }
    }
    int NTT_init(int n) {
        int lim = 1, l = 0;
        for (; lim < n; lim <<= 1, ++l);
        for (int i = 0; i < lim; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
        return lim;
    }
    void PolyInv(int *F, int *G, int n) {
        static int A[N], B[N];
        if (n == 1) { G[0] = qpow(F[0], mod - 2); return; }
        PolyInv(F, G, n >> 1);
        for (int i = 0; i < n; ++i) A[i] = F[i], B[i] = G[i];
        int lim = NTT_init(n << 1);
        NTT(A, lim, 1), NTT(B, lim, 1);
        for (int i = 0; i < lim; ++i) A[i] = 1ll * A[i] * B[i] % mod * B[i] % mod;
        NTT(A, lim, -1);
        for (int i = 0; i < n; ++i) G[i] = (2ll * G[i] - A[i] + mod) % mod;
        for (int i = 0; i < lim; ++i) A[i] = B[i] = 0;
    }
}

void Newton(int *F, int n) {
    static int F0[N], A[N], B[N], X[N], Y[N], iY[N];
    F[0] = 1;
    for (int lim = 2; lim <= n; lim <<= 1) {
        for (int i = 0; i < lim; ++i) F0[i] = F[i];
        for (int i = 0; i < lim; i += 2) A[i] = F[i / 2];
        for (int i = 0; i < lim; i += 3) B[i] = F[i / 3];
        
        Poly::NTT_init(lim << 1);
        Poly::NTT(F0, lim << 1, 1), Poly::NTT(A, lim << 1, 1);
        for (int i = 0; i < lim << 1; ++i)
            X[i] = (1ll * F0[i] * F0[i] % mod * F0[i] + 3ll * A[i] * F0[i]) % mod;
        Poly::NTT(X, lim << 1, -1), Poly::NTT(F0, lim << 1, -1), Poly::NTT(A, lim << 1, -1);
        for (int i = 0; i < lim; ++i) X[i] = (X[i] + 2ll * B[i]) % mod;
        for (int i = lim - 1; i; --i) X[i] = X[i - 1];
        X[0] = 6;
        for (int i = 0; i < lim; ++i) X[i] = (X[i] + 1ll * (mod - 6) * F0[i]) % mod;
        Poly::NTT(F0, lim << 1, 1);
        for (int i = 0; i < lim << 1; ++i) Y[i] = 3ll * F0[i] * F0[i] % mod;
        Poly::NTT(Y, lim << 1, -1), Poly::NTT(F0, lim << 1, -1);
        for (int i = 0; i < lim; ++i) Y[i] = (Y[i] + 3ll * A[i]) % mod;
        for (int i = lim - 1; i ; --i) Y[i] = Y[i - 1];
        Y[0] = mod - 6;

        for (int i = lim; i < lim << 1; ++i) X[i] = Y[i] = 0;
        Poly::PolyInv(Y, iY, lim);
        Poly::NTT_init(lim << 1);
        Poly::NTT(X, lim << 1, 1), Poly::NTT(iY, lim << 1, 1);
        for (int i = 0; i < lim << 1; ++i) X[i] = 1ll * X[i] * iY[i] % mod;
        Poly::NTT(X, lim << 1, -1);
        for (int i = 0; i < lim; ++i) F[i] = (F0[i] - X[i] + mod) % mod;
        
        for (int i = 0; i < lim << 1; ++i) F0[i] = A[i] = B[i] = X[i] = Y[i] = iY[i] = 0;
    }
}

int F[N], G[N], H[N];

int main() {
    int n = read(), lim = 1;
    for (; lim <= n; lim <<= 1);
    Newton(F, lim);

    Poly::NTT_init(lim << 1);
    Poly::NTT(F, lim << 1, 1);
    for (int i = 0; i < lim << 1; ++i) G[i] = 1ll * F[i] * F[i] % mod;
    Poly::NTT(G, lim << 1, -1), Poly::NTT(F, lim << 1, -1);
    for (int i = 0; i < lim; i += 2) G[i] = (G[i] + F[i / 2]) % mod;
    for (int i = 0; i < lim; ++i) G[i] = 1ll * G[i] * inv2 % mod;
    for (int i = lim - 1; i; --i) G[i] = G[i - 1];
    G[0] = 0;
    for (int i = lim; i < lim << 1; ++i) G[i] = 0;

    Poly::NTT(G, lim << 1, 1);
    for (int i = 0; i < lim << 1; ++i) H[i] = 1ll * G[i] * G[i] % mod;
    Poly::NTT(H, lim << 1, -1), Poly::NTT(G, lim << 1, -1);
    for (int i = 0; i < lim; i += 2) H[i] = (H[i] + G[i / 2]) % mod;
    for (int i = 0; i < lim; ++i) H[i] = 1ll * H[i] * inv2 % mod;

    for (int i = 2; i <= n; ++i) printf("%d\n", H[i]);
    return 0;
}
最后修改:2021 年 02 月 15 日 08 : 42 PM