[数学记录]P5900 无标号无根树计数

command_block · 2020-02-29 20:37:46 · 个人记录

神仙题.jpg = 大佬口中的经典题.bmp

广告 : 多项式计数杂谈

• 无标号有根树

设 \mathcal T 为无标号有根树的组合类，有 :

写成生成函数形式 ：

• 法1 : 牛顿迭代

有 T_*(x)={\bf T}(x)\pmod{x^n} 欲求 T(x)={\bf T}(x)\pmod{x^{2n}}

注意到，对于 T(x^k)\pmod{x^{2n}} ,当 k>2 时已经可以从 T^*(x) 中得到。

计算 P(x)=\exp\sum\limits_{i=2}\frac{T(x^i)}{i}

同样有 P^*(x)={\bf P}(x)\pmod{x^n} 欲求 P(x)={\bf P}(x)\pmod{x^{2n}}

注意到， \exp(a+b)=\exp(a)\exp(b)

则有 T(x)=P(x)\exp T(x)

接下来就可以牛顿迭代了,然而要写 \exp 所以在代码量和时效上都不优。

考虑去掉 \exp ,得 \ln T(x)=\ln P(x)+T(x)

即解 T(x)-\ln T(x)-\ln P(x)=0

这样子就不用EXP了,省下一堆常数&代码量。

• 法2 : 推式子+分治FFT

去掉\exp。

对 x 求导去掉 \ln ，注意链式法则。

• \dfrac{T(x^i)}{i}\dfrac{d}{dx}=\dfrac{T(x^i)}{i}\dfrac{d}{dx^i}\dfrac{dx^i}{dx}=x^{i-1}T'(x^i)

通分。

• 注意[x^n]xT'(x)=n[x^n]T(x)

设 G(x)=\sum\limits_{i=1}x^iT'(x^i) ,提取系数则有 G[n]=\sum\limits_{i|n}iT[i] ( 由于枚举约数,暴力计算是O(n\log n)的 )。

对 xT'(x)=T(x)+T(x)G(x) 提取系数:

剩下的就是经典的半在线卷积了。如何计算请见 : 半在线卷积小记

• 无标号无根树

我们考虑利用树的重心的性质，只考虑根是树的重心的方案数。

• ① n 是奇数

此时重心是惟一的，回忆重心的性质 : 每个子树的大小都不超过 \lfloor n/2\rfloor。

那么 “根不是重心” \ \Longleftrightarrow\ 有一个子树大小超过 \lfloor n/2\rfloor。

使用所有方案减去不合法的部分，可得

解释 : 把根和超大的子树之间的边隔开，分成两部分计算。

• ② n 是偶数

此时还能出现双重心的情况。即恰好有一条边将树分成两个大小为 n/2 的部分。

由于无标号，如果这条边两侧的树是等价的，那么根在那一半也是等价的，这部分没有重复计数。

反之，如果两个部分不同，就会计数两次。我们额外减去 \dbinom{T[n/2]}{2} 即可。

#include<algorithm>
#include<cstring>
#include<cstdio>
#define ll long long
#define ull unsigned long long
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n))
using namespace std;
const int _G=3,mod=998244353,Maxn=270000;
ll powM(ll a,int t=mod-2){
  ll ans=1;
  while(t){
    if(t&1)ans=ans*a%mod;
    a=a*a%mod;t>>=1;
  }return ans;
}
const int invG=powM(_G);
int tr[Maxn],tf;
void tpre(int n){
  if (tf==n)return ;tf=n;
  for(int i=0;i<n;i++)
    tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
}
void NTT(int *g,bool op,int n)
{
  tpre(n);
  static ull f[Maxn],w[Maxn]={1};
  for (int i=0;i<n;i++)f[i]=(((ll)mod<<5)+g[tr[i]])%mod;
  for(int l=1;l<n;l<<=1){
    ull tG=powM(op?_G:invG,(mod-1)/(l+l));
    for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
    for(int k=0;k<n;k+=l+l)
      for(int p=0;p<l;p++){
        int tt=w[p]*f[k|l|p]%mod;
        f[k|l|p]=f[k|p]+mod-tt;
        f[k|p]+=tt;
      }      
    if (l==(1<<10))
      for (int i=0;i<n;i++)f[i]%=mod;
  }if (!op){
    ull invn=powM(n);
    for(int i=0;i<n;++i)
      g[i]=f[i]%mod*invn%mod;
  }else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
void px(int *f,int *g,int n)
{for(int i=0;i<n;++i)f[i]=1ll*f[i]*g[i]%mod;}
int n,m,G[Maxn],F[Maxn],s1[Maxn],s2[Maxn],s3[Maxn];
void cdq(int l,int r)
{
  if (r==1)return ;
  if (l+1==r){
    if (l^1)F[l]=F[l]*powM(l-1)%mod;
    for (int i=l;i<n;i+=l)
      G[i]=(G[i]+1ll*l*F[l])%mod;
    return ;
  }int mid=(l+r)>>1,n=r-l;
  cdq(l,mid);
  if (l==0){
    cpy(s1,F+l,n/2);clr(s1+(n/2),n/2);NTT(s1,1,n);
    cpy(s2,G,n/2);clr(s2+(n/2),n/2);NTT(s2,1,n);
    px(s1,s2,n);NTT(s1,0,n);
    for (int i=n/2;i<n;i++)
      F[l+i]=(F[l+i]+s1[i])%mod;
  }else{
    cpy(s1,F+l,n/2);clr(s1+(n/2),n/2);NTT(s1,1,n);
    cpy(s2,G,n);NTT(s2,1,n);
    px(s1,s2,n);cpy(s3,s1,n);
    cpy(s1,G+l,n/2);clr(s1+(n/2),n/2);NTT(s1,1,n);
    cpy(s2,F,n);NTT(s2,1,n);
    px(s1,s2,n);for (int i=0;i<n;i++)s3[i]+=s1[i];
    NTT(s3,0,n);
    for (int i=n/2;i<n;i++)
      F[l+i]=(F[l+i]+s3[i])%mod;
  }cdq(mid,r);
}
int main()
{
  scanf("%d",&m);F[1]=1;
  for (n=1;n<=m;n<<=1);cdq(0,n);
  ll ans=F[m];
  for (int i=m/2+1;i<m;i++)
    ans=(ans-1ll*F[i]*F[m-i])%mod;
  if (!(m&1))ans=(ans-1ll*F[m/2]*(F[m/2]-1)/2)%mod;
  printf("%lld",(ans+mod)%mod);
  return 0;
}