「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 引理。最终方程为:
稍微移一下项,就可以得到:
这个方程有两种处理方法:一者,分治乘法,比较好想但实现难度未知;二者,我们用牛顿迭代解出这个方程。
这里有一个值得注意的细节,牛顿迭代实际上解的是如下的方程:
其中,\(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)\)。
小结:
- 注意生成函数在树计数、图计数问题中的应用。
- 注意有无标号的处理方法;一般来说,有标号会倾向于使用 EGF,而无标号会倾向于使用先标号,再去重或者排序等处理手段。
- 学习一下
完整的牛顿迭代法。其中的小细节不能忘了! - 牛顿迭代的倍增过程尤其注意次数;或者说,注意 \(\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;
}