「LOJ6538」烷基計數 加強版 加強版
阿新 • • 發佈：2022-02-22
題目
點這裡看題目。
省流版：求結點個數為 \(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)\)。
小結：
- 注意生成函式在樹計數、圖計數問題中的應用。
- 注意有無標號的處理方法；一般來說，有標號會傾向於使用 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;
}