烷基计数加强版加强版

Posted 2020-11-29 chasedeath

烷基计数加强版加强版

定义(F(x))表示答案多项式

枚举一个点为根，那么我们可以列出(F(x))的转移方程

(F(x)=xfrac{F(x)^3+3F(x^2)F(x)+2F(x^3)}{6}+1)

其中(F(x^2))表示选了两个相同大小，(F(x^3))表示选了三个大小相同，1表示空树

这个东西可以用(burnside)定理得出

burnside 定理

对于原来选择方案的一个等价变换(F ightarrow G)

如这个题中，三个子树编号为(1,2,3)，那么我们把三个子树任意变换到别的位置都是等价的

如((1,2,3) ightarrow (3,1,2))为一个等价的变化

(方案数=frac{对于每一个变换，变换之后方案完全相同的数量}{总的变换数量})

这个题所有的变换就是排列，一共有(6)个

对于这个题来说，排列之后不变，意味着排列之后这两个点完全相同

((1,2,3) ightarrow (1,2,3):F(x)^3)

((1,2,3) ightarrow (1,3,2):F(x)F(x^2)(3=2))

((1,2,3) ightarrow (2,1,3):F(x^2)F(x)(1=2))

((1,2,3) ightarrow (2,3,1):F(x^3)(1=2=3))

((1,2,3) ightarrow (3,1,2):F(x^3)(1=2=3))

((1,2,3) ightarrow (3,2,1):F(x^2)F(x)(1=3))

带入就得到了上面的转移方程式

[ ]

我们考虑解这个方程牛顿迭代法

[f(F(x))=F(x)-xfrac{F(x)^3+3F(x^2)F(x)+2F(x^3)}{6}-1 ]

[f‘(F(x))=1-frac{x}{6}(3F(x)^2+3F(x^2)) ]

预先求出(G(x) = F(x) (mod x^frac{n}{2}))

因为(F(x^2),F(x^3))可以通过(G(x))的系数直接得到,所以我们认为(F(x^2),F(x^3))是常数

设(F(x^2)=A,F(x^3)=B)

(f(F)=F-xfrac{F^3+3Acdot F+2B}{6}-1=0)

[f‘(F)=1-frac{x}{6}(3F^2+3A)=0 ]

求出多项式(f(F))在(G)上的泰勒展开

[f(F)=sum _0^infty frac{f^{(i)}(G)}{i!}(F-G)^i ]

(i=0,f(G))

(i=1,f‘(G)(F-G))

(i>1,(F-G)^iequiv 0 (mod x^n))

(0=f(G)+f‘(G)(F-G))

(F=G-frac{f(G)}{f‘(G)})

[F=G-frac{G-frac{x}{6}(G^3+3Acdot G+2B)-1}{1-frac{x}{6}(3G^2+3A)} ]

[F=frac{G-frac{x}{6}(3G^3+3Acdot G)-G+frac{x}{6}(G^3+3Acdot G+2B)+1}{1-frac{x}{6}(3G^2+3A)} ]

[F=frac{frac{x}{6}(2B-2 G^3)+1}{1-frac{x}{6}(3G^2+3A)} ]

代码，上面全部都是多项式的模板

#include<cstdio>
#include<algorithm>
#include<iostream>
#include<cctype>
#include<vector>
using namespace std;

#define reg register
typedef long long ll;
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)

#define pb push_back
template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); }
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); }

char IO;
template<class T=int> T rd(){
	T s=0;
	int f=0;
	while(!isdigit(IO=getchar())) if(IO==‘-‘) f=1;
	do s=(s<<1)+(s<<3)+(IO^‘0‘);
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

const int N=1<<21|10,P=998244353;

int n;

ll qpow(ll x,ll k) {
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res;
}

namespace Polynomial{
	typedef vector <int> Poly;
	#define Mod1(x) ((x>=P)&&(x-=P))
	#define Mod2(x) ((x<0)&&(x+=P))
	void Show(Poly a){ for(int i:a) printf("%d ",i); puts(""); }

	int rev[N];
	int PreMake(int n){
		int R=1,cc=-1;
		while(R<n) R<<=1,cc++;
		rep(i,1,R-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<cc);
		return R;
	}

	void NTT(int n,Poly &a,int f){
		rep(i,0,n-1) if(rev[i]<i) swap(a[i],a[rev[i]]);
		for(int i=1;i<n;i<<=1) {
			int w=qpow(f==1?3:(P+1)/3,(P-1)/i/2);
			for(int l=0;l<n;l+=i*2) {
				ll e=1;
				for(int j=l;j<l+i;++j,e=e*w%P) {
					int t=a[j+i]*e%P;
					a[j+i]=a[j]-t,Mod2(a[j+i]);
					a[j]+=t,Mod1(a[j]);
				}
			}
		}
		if(f==-1) {
			ll base=qpow(n,P-2);
			rep(i,0,n-1) a[i]=a[i]*base%P;
		}
	}

	Poly operator * (Poly a,Poly b){
		int n=a.size()+b.size()-1,R=PreMake(n);
		a.resize(R),b.resize(R);
		NTT(R,a,1),NTT(R,b,1);
		rep(i,0,R-1) a[i]=1ll*a[i]*b[i]%P;
		NTT(R,a,-1);
		a.resize(n);
		return a;
	}

	Poly operator + (Poly a,Poly b) { 
		int n=max(a.size(),b.size()); 
		a.resize(n),b.resize(n);
		rep(i,0,a.size()-1) a[i]+=b[i],Mod1(a[i]); 
		return a; 
	}
	Poly operator - (Poly a,Poly b) { 
		int n=max(a.size(),b.size()); 
		a.resize(n),b.resize(n);
		rep(i,0,a.size()-1) a[i]-=b[i],Mod2(a[i]); 
		return a; 
	} 

	Poly Inv(Poly a) {
		int n=a.size();
		if(n==1) { Poly tmp; tmp.pb(qpow(a[0],P-2)); return tmp; }
		Poly b=a; b.resize((n+1)/2); b=Inv(b);
		int R=PreMake(n<<1);
		a.resize(R),b.resize(R);
		NTT(R,a,1),NTT(R,b,1);
		rep(i,0,R-1) a[i]=(2-1ll*a[i]*b[i]%P+P)*b[i]%P;
		NTT(R,a,-1);
		a.resize(n);
		return a;
	}

	Poly operator / (vector <int> a,vector <int> b){
		reverse(a.begin(),a.end()),reverse(b.begin(),b.end());
		int n=a.size(),m=b.size();
		b.resize(n-m+1),b=Inv(b);
		a=a*b,a.resize(n-m+1);
		reverse(a.begin(),a.end());
		return a;
	}

	Poly operator % (vector <int> a,vector <int> b) { 
		a=a-a/b*b;
		a.resize(b.size()-1);
		return a;
	}

	Poly Sqrt(vector <int> a){
		int n=a.size();
		if(n==1) { Poly tmp; tmp.pb(1); return tmp; }
		Poly b=a; b.resize((n+1)/2),b=Sqrt(b),b.resize(n);
		Poly c=Inv(b);
		int R=PreMake(n*2);
		a.resize(R),c.resize(R);
		NTT(R,a,1),NTT(R,c,1);
		rep(i,0,R-1) a[i]=1ll*a[i]*c[i]%P;
		NTT(R,a,-1);
		a.resize(n);
		rep(i,0,n-1) a[i]=1ll*(P+1)/2*(a[i]+b[i])%P;
		return a;
	}

	Poly Deri(Poly a){
		rep(i,1,a.size()-1) a[i-1]=1ll*i*a[i]%P;
		a.pop_back();
		return a;
	}

	int Mod_Inv[N];
	Poly IDeri(Poly a) {
		Mod_Inv[0]=Mod_Inv[1]=1;
		rep(i,2,a.size()+1) Mod_Inv[i]=1ll*(P-P/i)*Mod_Inv[P%i]%P;
		a.pb(0);
		drep(i,a.size()-1,1) a[i]=1ll*a[i-1]*Mod_Inv[i]%P;
		a[0]=0;
		return a;
	}

	Poly Ln(Poly a){
		int n=a.size();
		a=Inv(a)*Deri(a);
		a.resize(n-1);
		return IDeri(a);
	}

	Poly Exp(Poly a){
		int n=a.size();
		if(n==1) { Poly tmp; tmp.pb(1); return tmp; }
		Poly b=a; b.resize((n+1)/2),b=Exp(b),b.resize(n);
		Poly c=Ln(b);
		int R=PreMake(n<<1);
		a.resize(R),b.resize(R),c.resize(R);
		NTT(R,b,1),NTT(R,a,1),NTT(R,c,1);
		rep(i,0,R-1) a[i]=1ll*b[i]*(1-c[i]+a[i]+P)%P;
		NTT(R,a,-1),a.resize(n);
		return a;
	}

}
using namespace Polynomial;

const int Inv6=qpow(6,P-2);
const int Inv3=(P+1)/3;
const int Inv2=(P+1)/2;

Poly Calc(int n){
	if(n==1) return Poly{1};
	int m=(n+1)>>1;
	Poly G=Calc(m),A,B;
	A.resize(n),B.resize(n);
	rep(i,0,(n-1)/2) A[i*2]=G[i];
	rep(i,0,(n-1)/3) B[i*3]=G[i];

	Poly x=Poly{1}+Poly{0,Inv3}*(B-G*G*G);
	Poly y=Poly{1}-Poly{0,Inv2}*(G*G+A);
	x=x*Inv(y);
	x.resize(n);
	return x;
}


int main(){
	int n=rd();
	Poly Res=Calc(n+1);
	printf("%d
",Res[n]);
}