LOJ 6538 烷基计数 （生成函数）（Burnside）（牛顿迭代）

传送门
题意：求儿子数不超过 3 的有根树个数
$a_x$ 表示有 x 个碳的方案树
直接上生成函数，$A(x)=\sum a_i{x}^i$

我们可以先强行钦定有 3 个儿子
如果只有一个或两个的话相当于那个儿子的 $size$ 为 0

题目要求不同构，考虑用 $Burnside$ 来计数
考虑同构的总方案数为每种置换下的不动点个数的平均值
对于置换 $(1,2,3)$，方案数是 3 个子树的拼接 $A(x{)}^3$
对于置换$(2,1,3),(1,3,2),(3,2,1)$，有两个树是同构的，方案数是 $3A(x)\ast A(x^2)$
对于置换$(3,1,2)(2,1,3)$，三个子树是同构的，方案数是 $2\ast A(x^3)$

于是有方程
$A(x)=1+x\frac{A(x{)}^3+3A(x)A(x^2)+2A(x^3)}{6}$

考虑牛顿迭代
求 $F(A(x))=1-A(x)+x\frac{A(x{)}^3+3A(x)A(x^2)+2A(x^3)}{6}$ 的 0 点
有 $A(x)=A_0(x)-\frac{F(A_0(x))}{F^{\prime }(A_0(x))}$

现在问题就是如何求导
发现牛顿迭代求出前 $x^{\frac{n}{2}}$ 的系数的时候， $x,A(x^2),A(x^3)$ 是已知的
将它们看成常数
那么 $F^{\prime }(A(x))=x\frac{3\ast A(x{)}^2+3\ast A(x^2)}{6}-1$

#include<bits/stdc++.h>
#define cs const
using namespace std;
cs int Mod = 998244353;
int add(int a, int b){ return a + b >= Mod ? a + b - Mod : a + b; }
int mul(int a, int b){ return 1ll * a * b % Mod; }
int dec(int a, int b){ return a - b < 0 ? a - b + Mod : a - b; }
int ksm(int a, int b){ int ans = 1; for(;b;b>>=1,a=mul(a,a)) if(b&1) ans = mul(ans, a); return ans; }
typedef vector<int> poly;
cs int inv6 = ksm(6, Mod-2);
int n, up, bit; poly rev;
void init(int deg){
	up = 1; bit = 0; while(up < deg) up <<= 1, ++bit; rev.resize(up);
	for(int i = 0; i < up; i++) rev[i] = (rev[i>>1]>>1)|((i&1)<<(bit-1));
} 
poly operator + (poly a, poly b){
	int deg = max(a.size(), b.size()); a.resize(deg); b.resize(deg);
	for(int i = 0; i < deg; i++) a[i] = add(a[i], b[i]); return a;
}
poly operator - (poly a, poly b){
	int deg = max(a.size(), b.size()); a.resize(deg); b.resize(deg);
	for(int i = 0; i < deg; i++) a[i] = dec(a[i], b[i]); return a;
}
poly operator * (poly a, int b){ for(int i = 0; i < (int)a.size(); i++) a[i] = mul(a[i], b); return a; }
void NTT(poly &a, int typ){
	for(int i = 0; i < up; i++) if(i < rev[i]) swap(a[i], a[rev[i]]);
	for(int i = 1; i < up; i <<= 1){
		int wn = ksm(3, (Mod-1)/(i<<1));
		for(int j = 0; j < up; j += (i<<1))
		for(int k = 0, w = 1; k < i; k++, w = mul(w, wn)){
			int x = a[k + j], y = mul(w, a[k + j + i]);
			a[k + j] = add(x, y); a[k + j + i] = dec(x, y);
		}
	} if(typ == -1){
		reverse(a.begin() + 1, a.end());
		for(int i = 0, inv = ksm(up, Mod-2); i < up; i++) a[i] = mul(a[i], inv);
	}
}
poly operator * (poly a, poly b){
	int deg = a.size() + b.size() - 1; init(deg);
	a.resize(up); b.resize(up); NTT(a, 1); NTT(b, 1);
	for(int i = 0; i < up; i++) a[i] = mul(a[i], b[i]);
	NTT(a, -1); a.resize(deg); return a;
}
poly inv(poly a){
	int deg = a.size(); poly b(1, ksm(a[0], Mod - 2)), c;
	for(int len = 4; len <= (deg << 2); len <<= 1){
		c = a; c.resize(len >> 1);
		init(len); c.resize(up); b.resize(up); NTT(b, 1); NTT(c, 1);
		for(int i = 0; i < up; i++) b[i] = mul(b[i], dec(2, mul(b[i], c[i])));
		NTT(b, -1); b.resize(len >> 1); 
	} b.resize(deg); return b;
}
poly Get(poly a, int coef, int len){
	poly b; b.resize(len);
	for(int i = 0; i < len; i++){
		if(i * coef < len) b[i * coef] = a[i];
	} return b;
}
poly Newton(int deg){
	poly b; b.push_back(1);
	for(int len = 2; len <= (deg << 1); len <<= 1){
		poly c = b * b * b + Get(b, 2, len) * b * 3 + Get(b, 3, len) * 2;
		c = c * inv6; c.insert(c.begin(), 0);
		c.resize(len); ++c[0]; c = c - b;
		poly d = b * b * 3 + Get(b, 2, len) * 3;
		d = d * inv6; d.insert(d.begin(), 0);
		d.resize(len); d[0] = dec(d[0], 1);
		b = b - c * inv(d); b.resize(len);
	} return b;
}
int main(){
	cin >> n; poly A = Newton(n); cout << A[n];
	return 0;
}