生成函数Euler变换学习笔记(无标号有根树计数)

本文探讨了无标号有根树的计数问题,通过引入Eular变换,详细推导了指数型生成函数的计算公式,并提供了一种高效的算法实现,复杂度达到O(nlog^2n)。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

众所周知,对于有标号计数的指数型生成函数 f(x)f(x)f(x),将其任意地进行无顺序的组合,得到的生成函数是exp(f(x))exp(f(x))exp(f(x))

而对于无标号计数的这样的组合,我们就需要引入一个叫 Eular\text{Eular}Eular 变换的东西 E(f(x))\mathcal{E}(f(x))E(f(x))

(这个符号打法:\mathcal(E)

为了方便,本文中记fif_ifi表示[xi]f(x)[x^i]f(x)[xi]f(x)。两种表述可能混用。

我们先引入一个实际问题:nnn个点无标号有根树的计数。两棵树被认为相同,当且仅当存在一种重标号方案,使得它们的根标号相同且边集相等。

对于一棵树,我们去掉它的根,然后就分成若干棵独立的子树。所以对无标号有根树的OGF\text{OGF}OGF F(x)F(x)F(x),有

F(x)=x⋅E(F(x))F(x)=x\cdot\mathcal E(F(x))F(x)=xE(F(x))

先思考这个问题的主要矛盾在什么地方:对于分出来的若干子树,如果有一些子树完全一致,那么它们只会贡献一次方案。但如果只是大小相同而本质不同,仍然算作不同的方案。

换句话说,最终的方案只和每种本质不同的子树使用的个数有关。考虑背包。

每种本质不同的大小为kkk的树构造完全背包的生成函数,即

11−xk\frac{1}{1-x^k}1xk1

那么我们可以得到一个简单粗暴的定义式:

E(F(x))=∏i=1∞1(1−xi)Fi\mathcal E(F(x))=\prod_{i=1}^{\infin}\frac 1{(1-x^i)^{F_i}}E(F(x))=i=1(1xi)Fi1

FiF_iFi是因为有FiF_iFi种本质不同大小为iii的物品

然后开始愉快地推式子

G(x)=E(F(x))=∏i=1∞1(1−xi)FiG(x)=\mathcal E(F(x))=\prod_{i=1}^{\infin}\frac 1{(1-x^i)^{F_i}}G(x)=E(F(x))=i=1(1xi)Fi1

既然有连乘,先求个ln⁡\lnln

ln⁡(G(x))=−∑i=1∞Filn⁡(1−xi)\ln (G(x))=-\sum_{i=1}^{\infin}F_i\ln(1-x^i)ln(G(x))=i=1Filn(1xi)

ln⁡\lnln还是看不顺眼,求个导

G′(x)G(x)=∑i=1∞Fiixi−11−xi\frac{G'(x)}{G(x)}=\sum_{i=1}^{\infin}F_i\frac{ix^{i-1}}{1-x^i}G(x)G(x)=i=1Fi1xiixi1

又看见了喜闻乐见的完全背包,直接展开

G′(x)G(x)=∑i=1∞Fi⋅ixi−1∑j=0∞xij\frac{G'(x)}{G(x)}=\sum_{i=1}^{\infin}F_i\cdot ix^{i-1}\sum_{j=0}^{\infin}x^{ij}G(x)G(x)=i=1Fiixi1j=0xij

G′(x)G(x)=∑i=1∞Fi⋅i∑j=0∞xi(j+1)−1\frac{G'(x)}{G(x)}=\sum_{i=1}^{\infin}F_i\cdot i\sum_{j=0}^{\infin}x^{i(j+1)-1}G(x)G(x)=i=1Fiij=0xi(j+1)1

G′(x)G(x)=∑i=1∞Fi⋅i∑j=1∞xij−1\frac{G'(x)}{G(x)}=\sum_{i=1}^{\infin}F_i\cdot i\sum_{j=1}^{\infin}x^{ij-1}G(x)G(x)=i=1Fiij=1xij1

再积分回来

ln⁡(G(x))=∑i=1∞Fi⋅i∑j=1∞xijij\ln(G(x))=\sum_{i=1}^{\infin}F_i\cdot i\sum_{j=1}^{\infin}\frac{x^{ij}}{ij}ln(G(x))=i=1Fiij=1ijxij

整理一下

ln⁡(G(x))=∑i=1∞Fi∑j=1∞xijj\ln(G(x))=\sum_{i=1}^{\infin}F_i\sum_{j=1}^{\infin}\frac{x^{ij}}jln(G(x))=i=1Fij=1jxij

ln⁡(G(x))=∑j=1∞1j∑i=1∞Fixij\ln(G(x))=\sum_{j=1}^{\infin}\frac 1j\sum_{i=1}^{\infin}F_ix^{ij}ln(G(x))=j=1j1i=1Fixij

ln⁡(G(x))=∑j=1∞F(xj)j\ln(G(x))=\sum_{j=1}^{\infin}\frac{F(x^j)}jln(G(x))=j=1jF(xj)

ln⁡(G(x))=∑i=1∞F(xi)i\ln(G(x))=\sum_{i=1}^{\infin}\frac{F(x^i)}iln(G(x))=i=1iF(xi)

得到 Eular\text{Eular}Eular 变换的计算式(雾)

E(F(x))=exp⁡(∑i=1∞F(xi)i)\mathcal E(F(x))=\exp(\sum_{i=1}^{\infin}\frac{F(x^i)}{i})E(F(x))=exp(i=1iF(xi))

然后回到原问题

F(x)=x⋅E(F(x))F(x)=x\cdot\mathcal E(F(x))F(x)=xE(F(x))

F(x)=x⋅exp⁡(∑i=1∞F(xi)i)F(x)=x\cdot \exp(\sum_{i=1}^{\infin}\frac{F(x^i)}{i})F(x)=xexp(i=1iF(xi))

两边求导(如果高数不好建议自己多写几步)

F′(x)=exp⁡(∑i=1∞F(xi)i)+x⋅exp⁡(∑i=1∞F(xi)i)⋅∑i=1∞ixi−1⋅F′(xi)iF'(x)=\exp(\sum_{i=1}^{\infin}\frac{F(x^i)}{i})+x\cdot \exp(\sum_{i=1}^{\infin}\frac{F(x^i)}{i})\cdot\sum_{i=1}^{\infin}\frac{ix^{i-1}\cdot F'(x^i)}iF(x)=exp(i=1iF(xi))+xexp(i=1iF(xi))i=1iixi1F(xi)

F′(x)=exp⁡(∑i=1∞F(xi)i)+F(x)∑i=1∞xi−1⋅F′(xi)F'(x)=\exp(\sum_{i=1}^{\infin}\frac{F(x^i)}{i})+F(x)\sum_{i=1}^{\infin}x^{i-1}\cdot F'(x^i)F(x)=exp(i=1iF(xi))+F(x)i=1xi1F(xi)

两边同时乘上xxx

x⋅F′(x)=F(x)+F(x)∑i=1∞xi⋅F′(xi)x\cdot F'(x)=F(x)+F(x)\sum_{i=1}^{\infin}x^i\cdot F'(x^i)xF(x)=F(x)+F(x)i=1xiF(xi)

#undef G

G(x)=∑i=1∞xi⋅F′(xi)G(x)=\sum_{i=1}^{\infin}x^i\cdot F'(x^i)G(x)=i=1xiF(xi)

x⋅F′(x)=F(x)+F(x)G(x)x\cdot F'(x)=F(x)+F(x)G(x)xF(x)=F(x)+F(x)G(x)

考虑求这东西的第xnx^nxn

nFn=Fn+∑k=1n−1FkGn−knF_n=F_n+\sum_{k=1}^{n-1}F_kG_{n-k}nFn=Fn+k=1n1FkGnk

Fn=1n−1∑k=1n−1FkGn−kF_{n}=\frac 1{n-1}\sum_{k=1}^{n-1}F_kG_{n-k}Fn=n11k=1n1FkGnk

如果我们知道GGG,就可以直接分治 NTT 算出FFF

现在考虑这个GGG

G(x)=∑i=1∞xi⋅F′(xi)G(x)=\sum_{i=1}^{\infin}x^i\cdot F'(x^i)G(x)=i=1xiF(xi)

Gn=∑i=1∞[xn−i]F′(xi)G_n=\sum_{i=1}^{\infin}[x^{n-i}]F'(x^i)Gn=i=1[xni]F(xi)

经过大胆的观察和仔细的猜想

Gn=∑d∣n(n−dd+1)Fn−dd+1G_n=\sum_{d\mid n}(\frac {n-d}d+1)F_{\frac {n-d}d+1}Gn=dn(dnd+1)Fdnd+1

Gn=∑d∣n(nd)FndG_n=\sum_{d\mid n}(\frac nd)F_{\frac nd}Gn=dn(dn)Fdn

Gn=∑d∣ndFdG_n=\sum_{d\mid n}d F_dGn=dndFd

然后每求出一个FFF,就可以调和级数复杂度求出之前的GGG,就可以分治 NTT 了。

但是分治时如果L=1L=1L=1GGG可能不够用。解决办法是先不计算不够的地方的贡献,在后面再把它加上。

具体的讲,如果我们要计算F1∼midF_{1\sim mid}F1midFmid+1∼RF_{mid+1\sim R}Fmid+1R的贡献,我们只计算iii不超过midmidmidGiG_iGi的贡献,这样少算了Gmid+1,R−1G_{mid+1,R-1}Gmid+1,R1的贡献。而这部分贡献可以通过RRR之前交换FFFGGG做一次类似的卷积计算出来。详见代码。

复杂度O(nlog⁡2n)O(n\log ^2n)O(nlog2n)

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#define MAXN 524288
using namespace std;
const int MOD=998244353;
typedef long long ll;
inline int add(const int& x,const int& y){return x+y>=MOD? x+y-MOD:x+y;}
inline int dec(const int& x,const int& y){return x<y? x-y+MOD:x-y;}
inline int qpow(int a,int p)
{
	int ans=1;
	while (p)
	{
		if (p&1) ans=(ll)ans*a%MOD;
		a=(ll)a*a%MOD;p>>=1;
	}
	return ans;
}
#define inv(x) qpow(x,MOD-2)
int rt[2][24],r[MAXN],l,lim;
inline void init(){lim=1<<l;for (int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));}
inline void NTT(int* a,int type)
{
	for (int i=0;i<lim;i++) if (i<r[i]) swap(a[i],a[r[i]]);
	for (int L=0;L<l;L++)
	{
		int mid=1<<L,len=mid<<1;
		ll Wn=rt[type][L+1];
		for (int s=0;s<lim;s+=len)
			for (int k=0,w=1;k<mid;k++,w=w*Wn%MOD)
			{
				int x=a[s+k],y=(ll)w*a[s+mid+k]%MOD;
				a[s+k]=add(x,y),a[s+mid+k]=dec(x,y);
			}
	}
	if (type)
	{
		int t=inv(lim);
		for (int i=0;i<lim;i++) a[i]=(ll)a[i]*t%MOD;
	}
}
int f[MAXN],g[MAXN],n;
int a[MAXN],b[MAXN];
void solve(int L,int R)
{
	if (L==R)
	{
		if (L>1) f[L]=(ll)f[L]*inv(L-1)%MOD;
		int t=(ll)L*f[L]%MOD;
		for (int i=L;i<=n;i+=L) g[i]=add(g[i],t);
		return;
	}
	int mid=(L+R)>>1;
	solve(L,mid);
	int len=R-L;
	for (l=0;(1<<l)<=(len<<1);l++);
	init();
	if (L==1)
	{
		for (int i=0;i<lim;i++) a[i]=b[i]=0;
		for (int i=L;i<=mid;i++) a[i]=f[i],b[i]=g[i];
		NTT(a,0);NTT(b,0);
		for (int i=0;i<lim;i++) a[i]=(ll)a[i]*b[i]%MOD;
		NTT(a,1);
		for (int i=mid+1;i<=R;i++) f[i]=add(f[i],a[i]);
	}
	else
	{
		for (int i=0;i<lim;i++) a[i]=b[i]=0;
		for (int i=L;i<=mid;i++) a[i-L]=f[i];
		for (int i=0;i<=len;i++) b[i]=g[i];
		NTT(a,0);NTT(b,0);
		for (int i=0;i<lim;i++) a[i]=(ll)a[i]*b[i]%MOD;
		NTT(a,1);
		for (int i=mid+1;i<=R;i++) f[i]=add(f[i],a[i-L]);
		for (int i=0;i<lim;i++) a[i]=b[i]=0;
		for (int i=L;i<=mid;i++) a[i-L]=g[i];
		for (int i=0;i<=len;i++) b[i]=f[i];
		NTT(a,0);NTT(b,0);
		for (int i=0;i<lim;i++) a[i]=(ll)a[i]*b[i]%MOD;
		NTT(a,1);
		for (int i=mid+1;i<=R;i++) f[i]=add(f[i],a[i-L]);
	}
	
	solve(mid+1,R);
}
int main()
{
	rt[0][23]=qpow(3,119),rt[1][23]=inv(rt[0][23]);
	for (int i=22;i>=0;i--)
	{
		rt[0][i]=(ll)rt[0][i+1]*rt[0][i+1]%MOD;
		rt[1][i]=(ll)rt[1][i+1]*rt[1][i+1]%MOD;
	}
	scanf("%d",&n);
	f[1]=1;
	solve(1,n);
	int ans=f[n];
	printf("%d\n",ans);
	return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值