无标号有根树计数与无标号无根树计数

OEIS首页的数列

参考博客

---

无标号无根树计数

考虑用总的无标号有根树方案数减去根不是重心的方案数,这样每个方案只会在根为重心时被统计一遍

设点数为$n$,已经求出$i$个点的无标号有根树方案数为$f_i$

$n$为奇数时

根据重心的性质,根节点不是重心,那么必然有一个子树大小$\ge \lceil \frac{n}{2}\rceil$

枚举这棵子树大小,答案等于$f_n-{\sum }_{k=\lceil \frac{n}{2}\rceil }^{n-1}{f}_k\times f_{n-k}$

$n$为偶数时

可能同时存在两个重心,其中一个重心为根

断开两个重心之间的边,产生两棵以重心为根的树

若这两棵树完全相同,那么在$f$中只会被统计一遍,不用减去

否则,这个方案在$f$中会被统计两遍,需要减去一遍,总的方案数是$C_{f_{\frac{n}{2}}}^2$种

两棵子树大小不相等的情况,按照原来的方式计算即可

答案等于$f_n-{\sum }_{k=\frac{n}{2}+1}^{n-1}{f}_k\times f_{n-k}-C_{f_{\frac{n}{2}}}^2$

无标号有根树计数

对于一个给定的根,这棵树的方案数相当于将除根外的所有点划分进若干棵不能为空的子树的方案数,那么构造方案数关于点数的$OGF$,设为$F$,那么有方程$F(x)=x\times euler(F(x))$

那么问题就是如何求解方程$F(x)-x\times euler(F(x))=0$

考虑倍增

直接倍增的做法

现在已经求出$f$满足$f(x)-x\times euler(f(x))\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^n$

考虑如何求出$F$满足$F(x)-x\times euler(F(x))\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^{2n}$

容易发现,$\forall k\ge 2$,$f(x^k)\mathrm{m}\mathrm{o}\mathrm{d}{x}^{2n}$都可以很方便的求出,且$f(x^k)\equiv F(x^k)\mathrm{m}\mathrm{o}\mathrm{d}{x}^{2n}$

所以当$f$为常数时,$P(x)=x\exp {\sum }_{k=2}\frac{f(x^k)}{k}\mathrm{m}\mathrm{o}\mathrm{d}{x}^{2n}$也为常数

再带入$Euler$变换的柿子,显然$F(x)-P(x)\exp F(x)\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^{2n}$

$$\frac{F(x)}{\exp F(x)}\equiv P(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^{2n}$$

再做牛顿迭代

现在已经求出$g$满足$\frac{g(x)}{\exp g(x)}\equiv P(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^n$

考虑如何求出$G$满足$\frac{G(x)}{\exp G(x)}\equiv P(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^{2n}$

求的是$H(x)=\frac{x}{\exp x}-P$,$H(F(x))\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^{2n}$的$F$

在$g(x)$处将$H$泰勒展开

$$H(G(x))\equiv H(g(x))+\sum _{i=1}\frac{H^{i^{\prime }}(g(x))(G(x)-g(x){)}^i}{i!}\mathrm{m}\mathrm{o}\mathrm{d}{x}^{2n}$$

$$\equiv H(g(x))+H^{\prime }(g(x))\times (G(x)-g(x))$$

$$\equiv 0$$

$$G(x)=g(x)-\frac{H(g(x))}{H^{\prime }(g(x))}$$

考虑如何求$H^{\prime }(g(x))$

$$H^{\prime }(x)=\frac{\exp x-x\exp x}{(\exp x{)}^2}=\frac{1-x}{\exp x}$$

$$H^{\prime }(g(x))=\frac{1-g(x)}{\exp g(x)}$$

(这里是对$g$求导而不是$x$)

$$G(x)=g(x)-\frac{\frac{g(x)}{\exp g(x)}-P(x)}{\frac{1-g(x)}{\exp g(x)}}$$

$$=g(x)+\frac{g(x)-P(x)\times \exp g(x)}{g(x)-1}$$

该做法是正确的,复杂度也是正确的,但是常数巨大,只能过$50\%$的点,$O(n\log n)$的复杂度$2\times 10^5$跑了$25s$

代码

/*
void checker(long long * f,long long *P,int n)
{
	static long long pd[800005],pdd[800005];
	memset(pd,0,sizeof(pd));
	memset(pdd,0,sizeof(pdd));
	EXP(f,pd,n);
	INV(pd,pdd,n);
	MUL(f,pdd,pd,n,n);
	printf("checker::f:");for(int i=0;i<n;i++)printf("%lld ",f[i]);printf("\n");
	printf("checker::P:");for(int i=0;i<n;i++)printf("%lld ",P[i]);printf("\n");
	printf("checker::frac{f}{exp f}-P:");for(int i=0;i<n;i++)printf("%lld ",pd[i]-P[i]);printf("\n");
}
*/
void Newton(long long * P,long long * G,int n)//答案存在G里
{
	if(n==1)return G[0]=P[0],void();
	static long long g[800005];
	Newton(P,g,n>>1);	
	static long long Exp[800005];
	EXP(g,Exp,n);
	MUL(Exp,P,Exp,n,n);
	for(int i=0;i<n;i++)Exp[i]=(g[i]-Exp[i]+p)%p;
	static long long Inv[800005];
	for(int i=0;i<n;i++)Inv[i]=0;
	g[0]=(g[0]+p-1)%p;
	INV(g,Inv,n);
	g[0]=(g[0]+1)%p;
	MUL(Exp,Inv,Exp,n,n);
	for(int i=0;i<n;i++)
		G[i]=(g[i]+Exp[i])%p;
}
long long inv[800005];
void Solve(long long * f,int n)//先要把n补成2的次幂 
{
	if(n==1)return f[0]=0,void();//\mod x意义下只要f\equiv 0即可
	static long long g[800005];
	Solve(g,n>>1);
	static long long P[800005];
	for(int i=0;i<n;i++)P[i]=0;
	for(int k=2;k<n;k++)
		for(int i=0;i*k<n;i++)
			P[i*k]=(P[i*k]+g[i]*inv[k])%p;
	static long long PP[800005];
	EXP(P,PP,n);
	for(int i=n-1;i>=1;i--)
		PP[i]=PP[i-1];
	PP[0]=0;
	for(int i=0;i<n;i++)f[i]=0;
	Newton(PP,f,n);	
}
/*
200000
174218497
*/

基于直接倍增做法的大优化

参考题解,可以发现以上做法常数巨大的原因并不是因为需要做$\exp$,而是倍增套倍增的层数太多了.
事实上,第一层倍增与第二层牛顿迭代完全可以合并.

回到一开始的问题.我们要求解的是方程$F(x)-x\times euler(F(x))\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^n$

现在已经求得方程$f(x)-x\times euler(f(x))\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n/2}$的解

显然,$F(x)$与$f(x)$的前$n/2$项是相同的

将方程转化为$F(x)-P(F(x))\exp F(x)\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^n$

其中$P(F(x))\equiv x\times \exp {\sum }_{k=2}\frac{F(x^k)}{k}\mathrm{m}\mathrm{o}\mathrm{d}{x}^n$

显然,$P(F(x))=P(f(x))$

牛顿迭代求解方程$F(x)-P(F(x))\exp F(x)\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^n$

现在已经求得方程$g(x)-P(F(x))\exp g(x)\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n/2}$

显然,$g(x)-P(f(x))\exp g(x)\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n/2}$,同时$f(x)-P(f(x))\exp f(x)\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n/2}$

不难发现,$g(x)\equiv f(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n/2}$

直接套用之前牛顿迭代的柿子
$$F(x)=f(x)+\frac{f(x)-P(f(x))\times \exp f(x)}{f(x)-1}$$

该做法是正确的,复杂度也是正确的,但是常数也巨大,也只能过$50\%$的点,不过$O(n\log n)$的复杂度$2\times 10^5$跑了$14s$,比之前快了不少

代码

void Solve(long long * f,int n)//先要把n补成2的次幂 
{
	if(n==1)return f[0]=0,void();//\mod x意义下只要f\equiv 0即可
	static long long g[800005];
	Solve(g,n>>1);
	static long long P[800005];
	for(int i=0;i<n;i++)P[i]=0;
	for(int k=2;k<n;k++)
		for(int i=0;i*k<n;i++)
			P[i*k]=(P[i*k]+g[i]*inv[k])%p;
	static long long PP[800005];
	EXP(P,PP,n);
	for(int i=n-1;i>=1;i--)
		P[i]=PP[i-1];
	P[0]=0;
	static long long Exp[800005];
	EXP(g,Exp,n);
	static long long Inv[800005];
	for(int i=0;i<n;i++)Inv[i]=0;
	g[0]=(g[0]+p-1)%p;
	INV(g,Inv,n);
	g[0]=(g[0]+1)%p;
	MUL(P,Exp,P,n,n);
	for(int i=0;i<n;i++)
		P[i]=(g[i]-P[i]+p)%p;
	MUL(P,Inv,P,n,n);
	for(int i=0;i<n;i++)
		f[i]=(g[i]+P[i])%p;
}

基于两边取$\ln$的做法

以上做法因为做了很多次$\exp$,常数十分吓人.如果我们可以把$\exp$换成$\ln$,常数会优很多

回到一开始的问题.我们要求解的是方程$F(x)-x\times euler(F(x))\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^n$

现在已经求得方程$f(x)-x\times euler(f(x))\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n/2}$的解

显然,$F(x)$与$f(x)$的前$n/2$项是相同的

将方程转化为$F(x)\equiv P(F(x))\exp F(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^n$

其中$P(F(x))\equiv x\times \exp {\sum }_{k=2}\frac{F(x^k)}{k}\mathrm{m}\mathrm{o}\mathrm{d}{x}^n$

对这个柿子两边取$\ln$
,$\ln F(x)\equiv {\sum }_{k=2}\frac{F(x^k)}{k}+\ln x+F(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^n$

显然,$\ln f(x)\equiv {\sum }_{k=2}\frac{F(x^k)}{k}+\ln x+f(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n/2}$

设$P(x)={\sum }_{k=2}\frac{F(x^k)}{k}+\ln x$

设$G(x)=\ln x-x-P$

$$G^{\prime }(x)=\frac{1}{x}-1$$
$$G^{\prime }(f(x))=\frac{1}{f(x)}-1$$

$$F(x)=f(x)-\frac{\ln f(x)-f(x)-P(x)}{\frac{1}{f(x)}-1}$$

容易发现,$f(x)$的常数项为$0$且一次方项不为$0$,并不能求$\ln$,同时$P(x)$中的$\ln x$也不能求出

那么,设$a(x)=\frac{f(x)}{x},A(x)=\frac{F(x)}{x}$

$$F(x)=f(x)-\frac{\ln a(x)+\ln x-f(x)-{\sum }_{k=2}\frac{f(x^k)}{k}-\ln x}{\frac{1}{f(x)}-1}$$
$$=f(x)-\frac{\ln a(x)-{\sum }_{k=1}\frac{f(x^k)}{k}}{\frac{1}{f(x)}-1}$$

$$F(x)=f(x)+\frac{\ln a(x)-{\sum }_{k=1}\frac{f(x^k)}{k}}{f(x)-1}\times f(x)$$

$$A(x)=a(x)-\frac{\ln a(x)-{\sum }_{k=1}\frac{f(x^k)}{k}}{\frac{1}{a(x)}-x}$$

这样就不用求$\exp$了

void Solve(long long * f,int n)//先要把n补成2的次幂 
{
	if(n==1)return f[0]=1,void();
	static long long g[1600005];
	Solve(g,n>>1);
	static long long P[1600005];
	for(int i=0;i<n;i++)P[i]=0;
	for(int k=1;k<n;k++)
		for(int i=1;i*k<n;i++)
			P[i*k]=(P[i*k]+g[i-1]*inv[k])%p;
	static long long Ln[1600005];
	LN(g,Ln,n);
	for(int i=0;i<n;i++)Ln[i]=(Ln[i]-P[i]+p)%p;
	static long long Inv[1600005],IInv[1600005];
	for(int i=0;i<n;i++)Inv[i]=0;
	for(int i=0;i<n;i++)IInv[i]=0;
	INV(g,Inv,n);
	Inv[1]=(Inv[1]-1+p)%p;
	INV(Inv,IInv,n);
	MUL(Ln,IInv,Ln,n,n);
	for(int i=0;i<n;i++)
		f[i]=(g[i]-Ln[i]+p)%p;
}
long long f[1600005];
int main()
{
	int n;
	read(n);
	for(int i=1;i<=1600000;i++)inv[i]=power(i);
	int N=1;
	while(N<n)N<<=1;
	Solve(f,N);
	for(int i=n;i>=1;i--)f[i]=f[i-1];
	f[0]=0;
	if(n&1)
	{
		long long ans=f[n];
		for(int i=n/2+1;i<=n-1;i++)
			ans=(ans+p-f[i]*f[n-i]%p)%p;
		printf("%lld\n",ans);
	}
	else
	{
		long long ans=(f[n]-(f[n/2]-1)*f[n/2]/2%p+p)%p;
		for(int i=n/2+1;i<=n-1;i++)
			ans=(ans+p-f[i]*f[n-i]%p)%p;
		printf("%lld\n",ans);
	}
	return 0;
}

未解之谜

1. 按照柿子$F(x)=f(x)+\frac{\ln \frac{f(x)}{x}-{\sum }_{k=1}\frac{f(x^k)}{k}}{f(x)-1}\times f(x)$牛顿迭代没法过,但是按照柿子$A(x)=a(x)-\frac{\ln a(x)-{\sum }_{k=1}\frac{f(x^k)}{k}}{\frac{1}{a(x)}-x}$牛顿迭代可以过,非常奇怪
2. $n$最大会等于$2^{18}$,但是把数组开到$8\times 10^5$没法过,开到$16\times 10^5$能过