「LOJ6538」烷基计数 加强版 加强版

题目

点这里看题目。

省流版：求结点个数为 \(n\) 的结点儿子数不超过 3 的无标号有根树个数。

对于 \(100\%\) 的数据，满足 \(1\le n\le 10^5\)。

分析

首先，遇到这种问题，不难猜想使用生成函数；并且，其余一大堆树计数的问题，都可以用生成函数解决。

设 \(G(x)\) 为问题的生成函数，也就是 \([x^n]G(x)\) 表示“结点个数为 \(n\) 且结点儿子数不超过 3 的无标号有根树个数”。我们认为 \([x^0]G(x)=1\)，也就是包含空树。

尝试列出一个关于 \(G\) 的方程：先确定一个根，而后考虑根的儿子。由于 \(G\) 包含空树，我们将所有的情况都看作是根有三个儿子。

我们不妨认为儿子之间有顺序，这样很容易计算出方案数为 \(G^3(x)\)；现在，如果两个方案是相同的，就存在一种儿子之间的映射，使得映射后子树对应相同。不难发现这就是在求 \(S_3\) 划分下等价类个数。所以可以联想到 Burnside 引理。最终方程为：

\[G(x)=1+\frac x 6(G^3(x)+3G(x^2)G(x)+2G(x^3)) \]

稍微移一下项，就可以得到：

\[\frac x 6G^3(x)+\left(\frac x 2G(x^2)-1\right)G(x)+\frac x 3G(x^3)+1=0 \]

这个方程有两种处理方法：一者，分治乘法，比较好想但实现难度未知；二者，我们用牛顿迭代解出这个方程。

这里有一个值得注意的细节，牛顿迭代实际上解的是如下的方程：

\[f(G,x)\equiv 0\pmod {x^n} \]

其中，\(f(G,x)\) 同时是关于 \(G,x\) 的多项式方程。其中用到的泰勒级数实际上也是在对 \(G\) 这个变量求偏导。所以最终的迭代过程为：

若：

\[f(G_{t},x)\equiv 0\pmod {x^{\lceil\frac n 2\rceil}}\\ \]

则：

\[G_{t+1}\equiv G_t-\frac{f(G_{t},x)}{\frac{\partial f}{\partial G}(G_t,x)}\pmod {x^n} \]

在这个情境下，有 \(f(G,x)=\frac x 6G^3(x)+\left(\frac x 2G(x^2)-1\right)G(x)+\frac x 3G(x^3)+1\)。但是存在 \(G(x^2),G(x^3)\) 这样的项，怎么处理呢？

注意到，在 \(G_t\) 的时候，实际上我们已经得到了 \(G_{t+1}(x^3)\equiv G_{t}(x^3),G_{t+1}(x^2)\equiv G_t(x^2)\pmod {x^n}\)（迭代过程中尤其注意次数）。所以我们可以认为 \(G_t(x^2),G_t(x^3)\) 为常量，求偏导得到的是 \(\frac{\partial f}{\partial G}(G,x)=\frac{x}{2}G^2(x)+\frac{x}{2}G(x^2)-1\)。

之后套用牛顿迭代即可，复杂度为 \(O(n\log n)\)。

小结：

1. 注意生成函数在树计数、图计数问题中的应用。
2. 注意有无标号的处理方法；一般来说，有标号会倾向于使用 EGF，而无标号会倾向于使用先标号，再去重或者排序等处理手段。
3. 学习一下完整的牛顿迭代法。其中的小细节不能忘了！
4. 牛顿迭代的倍增过程尤其注意次数；或者说，注意 \(\pmod {x^n}\)，这是一个很好的辅助条件。

代码

#include <cstdio>
#include <iostream>
#include <algorithm>

#define rep( i, a, b ) for( int i = (a) ; i <= (b) ; i ++ )
#define per( i, a, b ) for( int i = (a) ; i >= (b) ; i -- )

const int mod = 998244353, g = 3, phi = mod - 1;
const int inv2 = ( mod + 1 ) / 2, inv3 = ( mod + 1 ) / 3, inv6 = ( mod + 1 ) / 6;
const int MAXN = ( 1 << 18 ) + 5;

template<typename _T>
void read( _T &x ) {
    x = 0; char s = getchar(); bool f = false;
    while( s < '0' || '9' < s ) { f = s == '-', s = getchar(); }
    while( '0' <= s && s <= '9' ) { x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar(); }
    if( f ) x = -x;
}

template<typename _T>
void write( _T x ) {
    if( x < 0 ) putchar( '-' ), x = -x;
    if( 9 < x ) write( x / 10 );
    putchar( x % 10 + '0' );
}

int ans[MAXN];

int N;

inline int Qkpow( int, int );
inline int Inv( const int a ) { return Qkpow( a, mod - 2 ); }
inline int Mul( int x, const int v ) { return 1ll * x * v % mod; }
inline int Sub( int x, const int v ) { return ( x -= v ) < 0 ? x + mod : x; }
inline int Add( int x, const int v ) { return ( x += v ) >= mod ? x - mod : x; }

inline int Qkpow( int base, int indx ) {
    int ret = 1;
    while( indx ) {
        if( indx & 1 ) ret = Mul( ret, base );
        base = Mul( base, base ), indx >>= 1;
    }
    return ret;
}

namespace Basic {
    const int L = 18;

    int w[MAXN];

    void NTTInit( const int n = 1 << L ) {
        w[0] = 1, w[1] = Qkpow( g, phi >> L );
        rep( i, 2, n - 1 ) w[i] = Mul( w[i - 1], w[1] );
    }

    void NTT( int *coe, const int n, const int typ ) {
        for( int i = 1, j = 0 ; i < n ; i ++ ) {
            for( int d = n ; j ^= ( d >>= 1 ), ~ j & d ; );
            if( j < i ) std :: swap( coe[j], coe[i] );
        }
        int *wp, p, k;
        for( int s = 1 ; s < n ; s <<= 1 )
            for( int i = 0 ; i < n ; i += s << 1 ) {
                wp = w, p = ( 1 << L ) / ( s << 1 );
                for( int j = 0 ; j < s ; j ++, wp += p )
                    k = Mul( *wp, coe[i + j + s] ),
                    coe[i + j + s] = Sub( coe[i + j], k ),
                    coe[i + j] = Add( coe[i + j], k );
            }
        if( ~ typ ) return ;
        std :: reverse( coe + 1, coe + n );
        int inv = Inv( n ); rep( i, 0, n - 1 ) coe[i] = Mul( coe[i], inv );
    }
}

namespace PolyInv {
    int P[MAXN], U0[MAXN], V[MAXN], U[MAXN];

    void Newton( const int n ) {
        if( n == 1 ) {
            U0[0] = Inv( P[0] );
            return ;
        }
        int h = ( n + 1 ) >> 1, L; Newton( h );
        for( L = 1 ; L < ( n << 1 ) ; L <<= 1 );
        for( int i = 0 ; i < L ; i ++ ) U[i] = V[i] = 0;
        for( int i = 0 ; i < h ; i ++ ) U[i] = U0[i];
        for( int i = 0 ; i < n ; i ++ ) V[i] = P[i];
        Basic :: NTT( U, L, 1 ), Basic :: NTT( V, L, 1 );
        for( int i = 0 ; i < L ; i ++ ) U[i] = Mul( U[i], Sub( 2, Mul( U[i], V[i] ) ) );
        Basic :: NTT( U, L, -1 );
        for( int i = 0 ; i < n ; i ++ ) U0[i] = U[i];
    }

    void PolyInv( int *ret, const int *inp, const int n ) {
        rep( i, 0, n - 1 ) P[i] = U0[i] = V[i] = U[i] = 0;
        rep( i, 0, n - 1 ) P[i] = inp[i];
        Newton( n );
        rep( i, 0, n - 1 ) ret[i] = U0[i];
    }
}

namespace Solve {
    int U0[MAXN], U[MAXN], V[MAXN], W[MAXN];
    int P[MAXN], Q[MAXN];

    void Newton( const int n ) {
        if( n == 1 ) { U0[0] = 1; return ; }
        int h = ( n + 1 ) >> 1, L; Newton( h );
        for( L = 1 ; L < ( n << 1 ) ; L <<= 1 );
        for( int i = 0 ; i < L ; i ++ ) P[i] = Q[i] = U[i] = V[i] = W[i] = 0;
        for( int i = 0 ; i < h ; i ++ ) {
            U[i] = U0[i];
            if( i * 2 < n ) V[i * 2] = U0[i];
            if( i * 3 < n ) W[i * 3] = U0[i];
        }
        Basic :: NTT( U, L, 1 );
        Basic :: NTT( V, L, 1 );
        Basic :: NTT( W, L, 1 );
        int *wp = Basic :: w, p = ( 1 << Basic :: L ) / L;
        for( int i = 0 ; i < L ; i ++, wp += p ) {
            P[i] = Add( 1, Add( Mul( Mul( inv3, *wp ), W[i] ), 
                                Add( Mul( U[i], Sub( Mul( Mul( inv2, *wp ), V[i] ), 1 ) ),
                                     Mul( Mul( inv6, *wp ), Mul( Mul( U[i], U[i] ), U[i] ) ) ) ) );
            Q[i] = Sub( Add( Mul( Mul( inv2, *wp ), Mul( U[i], U[i] ) ),
                             Mul( Mul( inv2, *wp ), V[i] ) ), 1 );
        }
        Basic :: NTT( P, L, -1 );
        Basic :: NTT( Q, L, -1 );
        for( int i = n ; i < L ; i ++ ) P[i] = Q[i] = 0;
        PolyInv :: PolyInv( Q, Q, n );
        Basic :: NTT( P, L, 1 );
        Basic :: NTT( Q, L, 1 );
        for( int i = 0 ; i < L ; i ++ ) P[i] = Mul( P[i], Q[i] );
        Basic :: NTT( P, L, -1 );
        for( int i = 0 ; i < n ; i ++ ) U0[i] = Sub( U0[i], P[i] );
        for( int i = n ; i < L ; i ++ ) U0[i] = 0;
    }

    void Solve( int *ret, const int n ) {
        Newton( n );
        rep( i, 0, n - 1 ) ret[i] = U0[i];
    }
}

int main() {
    Basic :: NTTInit();
    read( N );
    Solve :: Solve( ans, N + 1 );
    write( ans[N] ), putchar( '\n' );
    return 0;
}