无标号有/无根树计数

\(\text{Problem}:\)无标号无根树计数

\(\text{Solution}:\)

引入 \(\text{Euler}\) 变换:将 \(\text{exp}\) 中的有标号元素替换为无标号元素。

考虑 \(\exp\) 的组合意义:将 \(n\) 个有标号的数分为若干非空组。设 \(F(x)\) 表示组内方案数的生成函数,有:

\[\exp(F(x))=\sum\limits_{i=0}^{\infty}\frac{F^{i}(x)}{i!} \]

对于无标号的问题,设大小为 \(i\) 的组有 \(f_{i}\) 种方案,有:

\[\mathcal{E}(F(x))=\prod\limits_{i=1}^{\infty}\frac{1}{(1-x^{i})^{f_{i}}} \]

但是这个定义式中需要利用 \(F(x)\) 的第 \(i\) 项来转移,不便于直接计算,故考虑优化其定义式。

套路的,利用 \(\ln\)\(\prod\) 转化为 \(\sum\) 的形式,有:

\[\ln \mathcal{E}(F(x))=\sum\limits_{i=1}^{\infty}f_{i}\ln(\frac{1}{1-x^{i}})\\ \]

根据付公主的背包(或拆分数)的套路,可以将 \(\ln(\frac{1}{1-x_{i}})\) 展开,故有:

\[\begin{aligned} \ln\mathcal{E}(F(x))&=\sum\limits_{i=1}^{\infty}f_{i}\sum\limits_{j=1}^{\infty}\frac{x^{ij}}{j}\\ &=\sum\limits_{j=1}^{\infty}\frac{1}{j}\sum\limits_{i=1}^{\infty}f_{i}(x^{j})^{i}\\ &=\sum\limits_{j=1}^{\infty}\frac{F(x^{j})}{j}\\ \mathcal{E}(F(x))&=\exp\left(\sum\limits_{j=1}^{\infty}\frac{F(x^{j})}{j}\right) \end{aligned} \]

回到本题。设 \(F(x)\)无标号有根树的生成函数。考虑去掉根结点后,是一个背包问题,故有:

\[F(x)=x\cdot \mathcal{E}(F(x)) \]

方法一:分治 \(\text{NTT}\)

对上式两边取 \(\ln\) 再求导,有:

\[\begin{aligned} \ln F(x)&=\ln x+\ln\mathcal{E}(F(x))\\ \frac{F'(x)}{F(x)}&=\frac{1}{x}+\sum\limits_{i=1}^{\infty}f_{i}\frac{ix^{i-1}}{1-x^{i}} \end{aligned} \]

两边同乘 \(xF(x)\),有:

\[\begin{aligned} xF'(x)&=F(x)+F(x)\sum\limits_{i=1}^{\infty}f_{i}\frac{ix^{i}}{1-x_{i}}\\ &=F(x)+F(x)\sum\limits_{i=1}^{\infty}x^{i}F'(x^{i}) \end{aligned} \]

考虑 \(F'(k-1)=F(k)\cdot k\),故设 \(G(x)=\sum\limits_{i=1}^{\infty}x^{i}F'(x^{i})\),有:

\[g_{n}=\sum\limits_{i\mid n}f_{i}\cdot i \]

同样的,有 \(xF'(x)[x^{n}]=nf_{n}\),故得到:

\[\begin{aligned} nf_{n}&=f_{n}+\sum\limits_{i=1}^{n-1}f_{i}g_{n-i}\\ f_{n}&=\frac{1}{n-1}\sum\limits_{i=1}^{n-1}f_{i}g_{n-i} \end{aligned} \]

至此,可以用分治 \(\text{NTT}\)\(O(n\log^2 n)\) 的时间复杂度内解决。这里与一般的分治 \(\text{FTT}\) 略有不同,具体实现可参考 command_block 半在线卷积小记

方法二:牛顿迭代

设函数 \(G(x)\),使得:

\[G(F(x))=\ln(\frac{F(x)}{x})-\sum\limits_{i=1}^{\infty}\frac{F(x^{i})}{i} \]

考虑递归处理。设已经求出 \(F_{0}(x)\),满足:

\[\begin{aligned} G(F(x))&\equiv 0\pmod {x^{n}}\\ G(F_{0}(x))&\equiv 0\pmod {x^{\lfloor\frac{n}{2}\rfloor}} \end{aligned} \]

发现 \(F(x^{i})\) 的形式决定了当 \(i\geq 2\) 时,\(F(x^{i})\equiv F_{0}(x^{i})\)。而这是一个调和级数,可以直接计算得出 \(i\geq 2\) 的答案,记为 \(w\),有:

\[\begin{aligned} G(F(x))&=\ln(\frac{F(x)}{x})-F(x)-w\\ G'(F(x))&=\frac{1}{F(x)}-1 \end{aligned} \]

由牛顿迭代,有:

\[\begin{aligned} F(x)&=F_{0}(x)-\frac{G(F_{0}(x))}{G'(F_{0}(x))}\\ &=F_{0}(x)-\frac{\ln(\frac{F_{0}(x)}{x})-F_{0}(x)-w}{\frac{1}{F_{0}(x)}-1} \end{aligned} \]

于是可以递归求解 \(F(x)\)。总时间复杂度 \(O(n\log n)\),但常数较大。

现在要求出的是无标号无根树的个数,那么我们钦定树的重心为根。下面进行分类讨论。

\(n\) 为奇数:

根据重心的性质,必然有一棵子树大小大于 \(\lfloor\frac{n}{2}\rfloor\)。设其为 \(i\),有:

\[f_{n}=\sum\limits_{i=\lfloor\frac{n}{2}\rfloor+1}^{n-1}f_{i}f_{n-i} \]

\(n\) 为偶数:

在奇数的基础上,当 \(2i=n\) 时,可能存在两个重心,故减去 \(\binom{f_{n/2}}{2}\)

\(\text{Code}:\) 分治 \(\text{NTT}\)

#include <bits/stdc++.h>
#pragma GCC optimize(3)
//#define int long long
#define ri register
#define mk make_pair
#define fi first
#define se second
#define pb push_back
#define eb emplace_back
#define is insert
#define es erase
#define vi vector<int>
#define vpi vector<pair<int,int>>
using namespace std; const int N=550010, Mod=998244353;
inline int read()
{
	int s=0, w=1; ri char ch=getchar();
	while(ch<'0'||ch>'9') { if(ch=='-') w=-1; ch=getchar(); }
	while(ch>='0'&&ch<='9') s=(s<<3)+(s<<1)+(ch^48), ch=getchar();
	return s*w;
}
int n;
int rev[N],r[24][2],iiv[N];
inline int ksc(int x,int p) { int res=1; for(;p;p>>=1, x=1ll*x*x%Mod) if(p&1) res=1ll*res*x%Mod; return res; }
inline void Get_Rev(int T) { for(ri int i=0;i<T;i++) rev[i]=(rev[i>>1]>>1)|((i&1)?(T>>1):0); }
inline void DFT(int T,vector<int> &s,int type)
{
	for(ri int i=0;i<T;i++) if(rev[i]<i) swap(s[i],s[rev[i]]);
	for(ri int i=2,cnt=1;i<=T;i<<=1,cnt++)
	{
		int wn=r[cnt][type];
		for(ri int j=0,mid=(i>>1);j<T;j+=i)
		{
			for(ri int k=0,w=1;k<mid;k++,w=1ll*w*wn%Mod)
			{
				int x=s[j+k], y=1ll*w*s[j+mid+k]%Mod;
				s[j+k]=x+y;
				if(s[j+k]>=Mod) s[j+k]-=Mod;
				s[j+mid+k]=x-y;
				if(s[j+mid+k]<0) s[j+mid+k]+=Mod;
			}
		}
	}
	if(!type) for(ri int i=0,inv=ksc(T,Mod-2);i<T;i++) s[i]=1ll*s[i]*inv%Mod;
}
inline void NTT(int n,int m,vector<int> &A,vector<int> &B)
{
	int len=n+m;
	int T=1;
	while(T<=len) T<<=1;
	Get_Rev(T);
	A.resize(T), B.resize(T);
	DFT(T,A,1), DFT(T,B,1);
	for(ri int i=0;i<T;i++) A[i]=1ll*A[i]*B[i]%Mod;
	DFT(T,A,0);
}
int f[N],g[N];
void Solve(int l,int r)
{
	if(r-l+1<=30)
	{
		for(ri int i=l;i<=r;i++)
		{
			for(ri int j=l;j<i;j++)
			{
				f[i]=(f[i]+1ll*f[j]*g[i-j]%Mod)%Mod;
				if(l^1) f[i]=(f[i]+1ll*g[j]*f[i-j]%Mod)%Mod;
			}
			if(i^1) f[i]=1ll*f[i]*iiv[i-1]%Mod;
			for(ri int j=i,w=1ll*i*f[i]%Mod;j<=n;j+=i) g[j]=(g[j]+w)%Mod;
		}
		return;
	}
	int mid=(l+r)/2;
	Solve(l,mid);
	if(l==1)
	{
		vector<int> A,B;
		A.resize(mid), B.resize(mid);
		for(ri int i=1;i<=mid;i++) A[i-1]=f[i], B[i-1]=g[i];
		NTT(mid-1,mid-1,A,B);
		for(ri int i=mid+1;i<=r;i++) f[i]=(f[i]+A[i-2])%Mod;
		Solve(mid+1,r);
		return;
	}
	vector<int> A,B;
	A.resize(mid-l+1), B.resize(r-l);
	for(ri int i=l;i<=mid;i++) A[i-l]=f[i];
	for(ri int i=1;i<=r-l;i++) B[i-1]=g[i];
	NTT(mid-l,r-l-1,A,B);
	vector<int> C,D;
	C.resize(mid-l+1), D.resize(r-l);
	for(ri int i=l;i<=mid;i++) C[i-l]=g[i];
	for(ri int i=1;i<=r-l;i++) D[i-1]=f[i];
	NTT(mid-l,r-l-1,C,D);
	for(ri int i=mid+1;i<=r;i++) f[i]=(f[i]+(A[i-l-1]+C[i-l-1])%Mod)%Mod;
	Solve(mid+1,r);
}
signed main()
{
	r[23][1]=ksc(3,119), r[23][0]=ksc(ksc(3,Mod-2),119);
	for(ri int i=22;~i;i--) r[i][0]=1ll*r[i+1][0]*r[i+1][0]%Mod, r[i][1]=1ll*r[i+1][1]*r[i+1][1]%Mod;
	iiv[1]=1;
	for(ri int i=2;i<N;i++) iiv[i]=1ll*(Mod-Mod/i)*iiv[Mod%i]%Mod;
	n=read();
	f[1]=1;
	Solve(1,n);
	int ans=f[n];
	for(ri int i=n/2+1;i<n;i++) ans=(ans-1ll*f[i]*f[n-i]%Mod+Mod)%Mod;
	if(!(n&1)) ans=(ans-1ll*f[n/2]*(f[n/2]-1)/2)%Mod;
	printf("%d\n",(ans+Mod)%Mod);
	return 0;
}

\(\text{Code}:\) 牛顿迭代

暂时咕咕咕
posted @ 2021-05-08 19:19  zkdxl  阅读(643)  评论(1编辑  收藏  举报