OEIS首页的数列

参考博客

无标号无根树计数

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

设点数为nnn,已经求出iii个点的无标号有根树方案数为fif_ifi

nnn为奇数时

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

枚举这棵子树大小,答案等于fn−∑k=⌈n2⌉n−1fk×fn−kf_n-\sum_{k=\lceil \frac{n}{2} \rceil}^{n-1}f_k\times f_{n-k}fn−∑k=⌈2n⌉n−1fk×fn−k

nnn为偶数时

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

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

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

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

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

答案等于fn−∑k=n2+1n−1fk×fn−k−Cfn22f_n-\sum_{k=\frac{n}{2}+1}^{n-1}f_k\times f_{n-k}-C_{f_{\frac{n}{2}}}^{2}fn−∑k=2n+1n−1fk×fn−k−Cf2n2

无标号有根树计数

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

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

考虑倍增

直接倍增的做法

现在已经求出fff满足f(x)−x×euler(f(x))≡0modxnf(x)-x\times euler(f(x))\equiv 0\mod x^nf(x)−x×euler(f(x))≡0modxn

考虑如何求出FFF满足F(x)−x×euler(F(x))≡0modx2nF(x)-x\times euler(F(x))\equiv 0\mod x^{2n}F(x)−x×euler(F(x))≡0modx2n

容易发现,∀k≥2\forall k\geq 2∀k≥2,f(xk)modx2nf(x^k)\mod x^{2n}f(xk)modx2n都可以很方便的求出,且f(xk)≡F(xk)modx2nf(x^k)\equiv F(x^k)\mod x^{2n}f(xk)≡F(xk)modx2n

所以当fff为常数时,P(x)=xexp⁡∑k=2f(xk)kmodx2nP(x)=x\exp\sum_{k=2}\frac{f(x^k)}{k}\mod x^{2n}P(x)=xexp∑k=2kf(xk)modx2n也为常数

再带入EulerEulerEuler变换的柿子,显然F(x)−P(x)exp⁡F(x)≡0modx2nF(x)-P(x)\exp F(x)\equiv 0\mod x^{2n}F(x)−P(x)expF(x)≡0modx2n

F(x)exp⁡F(x)≡P(x)modx2n\frac{F(x)}{\exp F(x)}\equiv P(x)\mod x^{2n}expF(x)F(x)≡P(x)modx2n

再做牛顿迭代

现在已经求出ggg满足g(x)exp⁡g(x)≡P(x)modxn\frac{g(x)}{\exp g(x)}\equiv P(x)\mod x^nexpg(x)g(x)≡P(x)modxn

考虑如何求出GGG满足G(x)exp⁡G(x)≡P(x)modx2n\frac{G(x)}{\exp G(x)}\equiv P(x)\mod x^{2n}expG(x)G(x)≡P(x)modx2n

求的是H(x)=xexp⁡x−PH(x)=\frac{x}{\exp x}-PH(x)=expxx−P,H(F(x))≡0modx2nH(F(x))\equiv 0\mod x^{2n}H(F(x))≡0modx2n的FFF

在g(x)g(x)g(x)处将HHH泰勒展开

H(G(x))≡H(g(x))+∑i=1Hi′(g(x))(G(x)−g(x))ii!modx2nH(G(x))\equiv H(g(x))+\sum_{i=1}\frac{H^{i'}(g(x))(G(x)-g(x))^i}{i!}\mod x^{2n}H(G(x))≡H(g(x))+i=1∑i!Hi′(g(x))(G(x)−g(x))imodx2n

≡H(g(x))+H′(g(x))×(G(x)−g(x))\equiv H(g(x))+H'(g(x))\times (G(x)-g(x))≡H(g(x))+H′(g(x))×(G(x)−g(x))

≡0\equiv 0≡0

G(x)=g(x)−H(g(x))H′(g(x))G(x)=g(x)-\frac{H(g(x))}{H'(g(x))}G(x)=g(x)−H′(g(x))H(g(x))

考虑如何求H′(g(x))H'(g(x))H′(g(x))

H′(x)=exp⁡x−xexp⁡x(exp⁡x)2=1−xexp⁡xH'(x)=\frac{\exp x-x\exp x}{(\exp x)^2}=\frac{1-x}{\exp x}H′(x)=(expx)2expx−xexpx=expx1−x

H′(g(x))=1−g(x)exp⁡g(x)H'(g(x))=\frac{1-g(x)}{\exp g(x)}H′(g(x))=expg(x)1−g(x)

(这里是对ggg求导而不是xxx)

G(x)=g(x)−g(x)exp⁡g(x)−P(x)1−g(x)exp⁡g(x)G(x)=g(x)-\frac{\frac{g(x)}{\exp g(x)}-P(x)}{\frac{1-g(x)}{\exp g(x)}}G(x)=g(x)−expg(x)1−g(x)expg(x)g(x)−P(x)

=g(x)+g(x)−P(x)×exp⁡g(x)g(x)−1=g(x)+\frac{g(x)-P(x)\times \exp g(x)}{g(x)-1}=g(x)+g(x)−1g(x)−P(x)×expg(x)

该做法是正确的,复杂度也是正确的,但是常数巨大,只能过50%50\%50%的点,O(nlog⁡n)O(n\log n)O(nlogn)的复杂度2×1052\times 10^52×105跑了25s25s25s

代码

/*
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⁡\expexp,而是倍增套倍增的层数太多了.
事实上,第一层倍增与第二层牛顿迭代完全可以合并.

回到一开始的问题.我们要求解的是方程F(x)−x×euler(F(x))≡0modxnF(x)-x\times euler(F(x))\equiv 0\mod x^nF(x)−x×euler(F(x))≡0modxn

现在已经求得方程f(x)−x×euler(f(x))≡0modxn/2f(x)-x\times euler(f(x))\equiv 0\mod x^{n/2}f(x)−x×euler(f(x))≡0modxn/2的解

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

将方程转化为F(x)−P(F(x))exp⁡F(x)≡0modxnF(x)-P(F(x))\exp F(x)\equiv 0\mod x^nF(x)−P(F(x))expF(x)≡0modxn

其中P(F(x))≡x×exp⁡∑k=2F(xk)kmodxnP(F(x))\equiv x\times \exp\sum_{k=2}\frac{F(x^k)}{k}\mod x^nP(F(x))≡x×exp∑k=2kF(xk)modxn

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

牛顿迭代求解方程F(x)−P(F(x))exp⁡F(x)≡0modxnF(x)-P(F(x))\exp F(x)\equiv 0\mod x^nF(x)−P(F(x))expF(x)≡0modxn

现在已经求得方程g(x)−P(F(x))exp⁡g(x)≡0modxn/2g(x)-P(F(x))\exp g(x)\equiv 0\mod x^{n/2}g(x)−P(F(x))expg(x)≡0modxn/2

显然,g(x)−P(f(x))exp⁡g(x)≡0modxn/2g(x)-P(f(x))\exp g(x)\equiv 0\mod x^{n/2}g(x)−P(f(x))expg(x)≡0modxn/2,同时f(x)−P(f(x))exp⁡f(x)≡0modxn/2f(x)-P(f(x))\exp f(x)\equiv 0\mod x^{n/2}f(x)−P(f(x))expf(x)≡0modxn/2

不难发现,g(x)≡f(x)modxn/2g(x)\equiv f(x)\mod x^{n/2}g(x)≡f(x)modxn/2

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

该做法是正确的,复杂度也是正确的,但是常数也巨大,也只能过50%50\%50%的点,不过O(nlog⁡n)O(n\log n)O(nlogn)的复杂度2×1052\times 10^52×105跑了14s14s14s,比之前快了不少

代码

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⁡\lnln的做法

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

回到一开始的问题.我们要求解的是方程F(x)−x×euler(F(x))≡0modxnF(x)-x\times euler(F(x))\equiv 0\mod x^nF(x)−x×euler(F(x))≡0modxn

现在已经求得方程f(x)−x×euler(f(x))≡0modxn/2f(x)-x\times euler(f(x))\equiv 0\mod x^{n/2}f(x)−x×euler(f(x))≡0modxn/2的解

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

将方程转化为F(x)≡P(F(x))exp⁡F(x)modxnF(x)\equiv P(F(x))\exp F(x)\mod x^nF(x)≡P(F(x))expF(x)modxn

其中P(F(x))≡x×exp⁡∑k=2F(xk)kmodxnP(F(x))\equiv x\times \exp\sum_{k=2}\frac{F(x^k)}{k}\mod x^nP(F(x))≡x×exp∑k=2kF(xk)modxn

对这个柿子两边取ln⁡\lnln
,ln⁡F(x)≡∑k=2F(xk)k+ln⁡x+F(x)modxn\ln F(x)\equiv \sum_{k=2}\frac{F(x^k)}{k}+\ln x+F(x)\mod x^{n}lnF(x)≡∑k=2kF(xk)+lnx+F(x)modxn

显然,ln⁡f(x)≡∑k=2F(xk)k+ln⁡x+f(x)modxn/2\ln f(x)\equiv \sum_{k=2}\frac{F(x^k)}{k}+\ln x+f(x)\mod x^{n/2}lnf(x)≡∑k=2kF(xk)+lnx+f(x)modxn/2

设P(x)=∑k=2F(xk)k+ln⁡xP(x)=\sum_{k=2}\frac{F(x^k)}{k}+\ln xP(x)=∑k=2kF(xk)+lnx

设G(x)=ln⁡x−x−PG(x)=\ln x-x-PG(x)=lnx−x−P

G′(x)=1x−1G'(x)=\frac{1}{x}-1G′(x)=x1−1
G′(f(x))=1f(x)−1G'(f(x))=\frac{1}{f(x)}-1G′(f(x))=f(x)1−1

F(x)=f(x)−ln⁡f(x)−f(x)−P(x)1f(x)−1F(x)=f(x)-\frac{\ln f(x)-f(x)-P(x)}{\frac{1}{f(x)}-1}F(x)=f(x)−f(x)1−1lnf(x)−f(x)−P(x)

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

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

F(x)=f(x)−ln⁡a(x)+ln⁡x−f(x)−∑k=2f(xk)k−ln⁡x1f(x)−1F(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)=f(x)−f(x)1−1lna(x)+lnx−f(x)−∑k=2kf(xk)−lnx
=f(x)−ln⁡a(x)−∑k=1f(xk)k1f(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)1−1lna(x)−∑k=1kf(xk)

F(x)=f(x)+ln⁡a(x)−∑k=1f(xk)kf(x)−1×f(x)F(x)=f(x)+\frac{\ln a(x)-\sum_{k=1}\frac{f(x^k)}{k}}{f(x)-1}\times f(x)F(x)=f(x)+f(x)−1lna(x)−∑k=1kf(xk)×f(x)

A(x)=a(x)−ln⁡a(x)−∑k=1f(xk)k1a(x)−xA(x)=a(x)-\frac{\ln a(x)-\sum_{k=1}\frac{f(x^k)}{k}}{\frac{1}{a(x)}-x}A(x)=a(x)−a(x)1−xlna(x)−∑k=1kf(xk)

这样就不用求exp⁡\expexp了

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;
}

未解之谜

按照柿子F(x)=f(x)+ln⁡f(x)x−∑k=1f(xk)kf(x)−1×f(x)F(x)=f(x)+\frac{\ln \frac{f(x)}{x}-\sum_{k=1}\frac{f(x^k)}{k}}{f(x)-1}\times f(x)F(x)=f(x)+f(x)−1lnxf(x)−∑k=1kf(xk)×f(x)牛顿迭代没法过,但是按照柿子A(x)=a(x)−ln⁡a(x)−∑k=1f(xk)k1a(x)−xA(x)=a(x)-\frac{\ln a(x)-\sum_{k=1}\frac{f(x^k)}{k}}{\frac{1}{a(x)}-x}A(x)=a(x)−a(x)1−xlna(x)−∑k=1kf(xk)牛顿迭代可以过,非常奇怪
nnn最大会等于2182^{18}218,但是把数组开到8×1058\times 10^58×105没法过,开到16×10516\times 10^516×105能过

无标号有根树计数与无标号无根树计数相关推荐

生成函数Euler变换学习笔记（无标号有根树计数）
众所周知,对于有标号计数的指数型生成函数 f(x)f(x)f(x),将其任意地进行无顺序的组合,得到的生成函数是exp(f(x))exp(f(x))exp(f(x)). 而对于无标号计数的这样的组合, ...
android 无埋点简书,无埋点README
无埋点编码规范无埋点方案基于窗口回调(Window.Callback)机制.BaseActivity中集成了自动打点相关逻辑.但由于dialog和activity实现机制不一样.为了dialog同样 ...
android studio logcat 无筛选显示全部日志无应用包名区分
android studio logcat 无筛选显示全部日志无应用包名区分不显示所有应用出现这个情况后很多同学无法解决,重启adb,重启studio,重启电脑,都是没用的... 其实是有个开 ...
ArcGIS\QGIS无插件加载（无偏移）MapBox高清影像图
喜欢就关注我们吧! 首先介绍一下MapBOX. Mapbox 是用于移动和 Web 应用程序的位置数据平台.用户可以使用Mapbox Studio创建一个自定义.交互式的地图,然后可以将这些自定义的地 ...
何水无鱼？何山无石？何人无父？何女无夫？何树无枝？何城无市？
何水无鱼? 共有八个答案: 子. 南无阿弥托佛南水无鱼,无山无石,阿人无父,弥女无夫,陀树无枝,佛城无市.南无阿弥陀佛. 何水无鱼?何山无石?何人无父?何女无夫?何树无枝?何城无市?原于释迦凡尘 ...
姓名学中萍字无根怎么解释_无根Buildah的工作原理：在非特权环境中构建容器
姓名学中萍字无根怎么解释在以前的文章中,包括无根Podman如何工作? ,我谈到了Podman ,该工具使用户可以管理Pod,容器和容器图像. Buildah是用于构建与Podman互补的Open ...
chm打开秒退_Mac_Mac电脑程序无响应怎么办?Mac程序无响应解决方法，虽然Mac电脑一向以运行稳定、 - phpStudy...
Mac电脑程序无响应怎么办?Mac程序无响应解决方法虽然Mac电脑一向以运行稳定.流畅而著称,但Mac电脑运行时间长了,难免也会遇到程序卡死无响应.一直"转菊花"的情况,可能是由 ...
java无侵入埋点 point_无侵入埋点
埋点是一种了解用户行为,分析用户行为,提高用户体验的一种方式. 常见的解决方案有三种,代码埋点.可视化埋点.和无埋点三种. 代码埋点主要就是通过手写代码的方式来埋点,能很精确的在需要埋点的地方,添加代 ...
分享20个无版权的高清无码图库站
今天这组网站比较有特色,有专门分享美食图片的,有专门分享复古图片的,各领风骚,质量都是一顶一的棒.下面就是20个无版权的高清无码图库站,记得收藏啊. 您可能感兴趣的相关文章 35款精致的 CSS3 和 ...
苹果xr怎么截屏_手机资讯:iPhone XR更新系统后无信号怎么办iPhone XR无信号解决办法...
如今使用IT数码设备的小伙伴们是越来越多了,那么IT数码设备当中是有很多知识的,这些知识很多小伙伴一般都是不知道的,就好比最近就有很多小伙伴们想要知道iPhone XR更新系统后无信号怎么办iPhon ...

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