题解 P5900 【无标号无根树计数】

· · 题解

题意：题如其名。

原题

这里提供一份被分治 \text{NTT} 两只 \log 吊着锤的单 \log 牛顿迭代题解。

首先是 \text{Euler} 变换：

让我们把一些无标号物品分成若干组（非空），大小为 k 的组内部有 f_k 种方案，求方案数。

对于无标号问题，如果出现大小与方案皆相同的两组，则它们没有区别。也就是说，对于大小、方案皆相同的几组，形成它们的方案数为 1 。

那么，枚举大小 k ，固定一种方案，则得到一个生成函数 \sum_{i=0}\limits x^{ki} = \dfrac{1}{1-x^k}

对于每个大小 k，方案数为 f_k ，那么最后总方案数的生成函数就是

\prod_{k=1} (\frac{1}{1-x^k})^{f_k}

即

\text{Euler}(f(x)) = \prod_{k=1} (\frac{1}{1-x^k})^{[x^k]f(x)}

我们可以把它化成另一个更实用的形式：

\text{Euler}(f(x)) = \exp (\sum_{k=1} \ln( (\frac{1}{1-x^k})^{[x^k]f(x)}))=\exp (\sum_{k=1}[x^k]f(x) \cdot \ln(\frac{1}{1-x^k}))

做过付公主的背包的同志应该知道：

\ln(\frac{1}{1-x^k})=\sum_{t=1} \frac{x^{kt}}{t}

如果没做过，可以先将其求导，一顿操作，然后积回来得到上式。

那么，

\text{Euler}(f(x))=\exp(\sum_{k=1}[x^k]f(x)\cdot (\sum_{t=1}\frac{x^{kt}}{t}))

改变求和顺序，先枚举 t ，可以得到：

\text{Euler}(f(x))=\exp(\sum_{t=1} \frac{f(x^t)}{t})

先来看看无标号有根树的计数：

g(x) = x \cdot \text{Euler}(g(x)) \frac{g(x)}{x} = \exp(\sum_{t=1}\frac{g(x^{t})}{t})

两边一起取 \ln，得到

\ln(\frac{g(x)}{x})=\sum_{t=1}\frac{g(x^t)}{t}

设函数 f ，且令

f(g(x))=\ln(\frac{g(x)}{x})-\sum_{t=1}\frac{g(x^t)}{t}

对这个东西上牛顿迭代：求 g(x) 使 f(g(x))\equiv 0 \pmod{x^n}

假设我们递归的求出了 g_0(x) ，使 f(g_0(x))\equiv 0 \pmod{x^{\lceil \frac{n}{2} \rceil}}

注意到 g(x) 和 g_0(x) 的前 \lceil \frac{n}{2} \rceil 项是一致的，而且 \dfrac{g(x^t)}{t} 只会跟 g(x) 的前 \lfloor \dfrac{n}{t} \rfloor 项有关。

也就是说，对于 t>=2 ，\dfrac{g(x^t)}{t} 已经被决定了。因为此处的 g 完全可以用 g_0 换掉。那么，令

\alpha = \sum_{t=2}{\frac{g_0(x^t)}{t}}

这个东西可以调和级数得到。那么：

f(g(x))=\ln(\frac{{g(x)}}{{x}})-g(x)-\alpha

其中 \alpha 与 g(x) 无关。\dfrac{1}{x} 与 g(x) 也无关。

f'(g(x))=\frac{1}{x} \cdot \frac{1}{\frac{g(x)}{x}}-1

上牛迭公式：

g(x)\equiv g_0(x)-\frac{f(g_0(x))}{f'(g_0(x))} \pmod {x^n}

令 g_1(x)=\dfrac{g_0(x)}{x}，g_2(x)=\dfrac{g(x)}{x}

g(x)\equiv g_0(x)-\frac{\ln(g_1(x))-g_0(x)-\alpha}{\frac{1}{x}\cdot \frac{1}{g_1(x)}-1} \pmod{x^n}

两边一起除掉 x，得到：

g_2(x) \equiv g_1(x)-\frac{\ln(g_1(x))-g_0(x)-\alpha}{\frac{1}{g_1(x)}-x} \pmod{x^n}

于是，可以对 \dfrac{g(x)}{x} 进行牛顿迭代求解。从而得到 g(x)。

这里复杂度是 n \log n 的，但是常数巨大。

最后，考虑怎么由无标号有根树得到无标号无根树：让每种情况只在根为重心时被统计一次。

也就是说，用总方案数减去根不是重心的方案数。

既然根不是重心，那么根一定会有且只有一个子树的大小大于等于 \lfloor \dfrac{n}{2} \rfloor+1。

枚举这个大小 k ，减掉 [x^k]g(x)\cdot [x^{n-k}]g(x) 就好。

同时，如果 n 是个偶数，那它有可能会有两个子树大小为 \dfrac{n}{2} 的重心。

如果这两个重心的子树内形状一致，则会被正确统计一次。

否则，会统计两次。于是要减掉多出来的那一半，也就是

\frac{[x^{\frac{n}{2}}]g(x)\cdot ([x^{\frac{n}{2}}]g(x)-1)}{2}=\binom{[x^{\frac{n}{2}}]g(x)}{2}

至此结束。这个做法避免了 \exp 这个常数要人命的东西，然而常数依然很大，大得被两只 \log 吊起来锤。

（可能是我自带大常数……板子的实现太不优秀了）

剩下的，见代码：

#include <cstdio>
#include <cstring>

const int maxn = 6e5+5 , mod = 998244353;

template<typename T>
inline T max(const T &a,const T &b){
    return a>b?a:b;
}

template<typename T>
inline void swap(T &a,T &b){
    T temp=a;a=b;b=temp;
}

inline int mul(int a,int b){return 1ll*a*b%mod;}

int fastpow(int x,int y){
    if(y==0) return 1;
    int tmp=fastpow(x,y>>1);
    return y&1?mul(mul(tmp,tmp),x):mul(tmp,tmp);
}

namespace Poly{
    const int gate = 3,invg = fastpow(gate,mod-2);
    int lim,lg,sing[maxn];
    void init(int len){
        lim = 1,lg = 0;while(lim<len) lim<<=1,lg++;
        for(int i=1;i<lim;i++) sing[i]=(sing[i>>1]>>1)|((i&1)<<(lg-1));
    }
    void NTT(int *Steins,int type){
        for(int i=0;i<lim;i++) if(i<sing[i]) swap(Steins[i],Steins[sing[i]]);
        for(int len=1;len<lim;len<<=1){
            int unit=fastpow(type==1?gate:invg,(mod-1)/(len<<1));
            for(int i=0;i<lim;i+=(len<<1)){
                int w = 1;
                for(int k=0;k<len;k++,w=mul(w,unit)){
                    int x=Steins[i+k],y=mul(w,Steins[i+k+len]);
                    Steins[i+k]=(x+y)%mod,Steins[i+k+len]=(x-y)%mod;
                }
            }
        }
        if(type!=1){int Okabe=fastpow(lim,mod-2);for(int i=0;i<lim;i++) Steins[i]=mul(Steins[i],Okabe);}
    }
    int g[maxn];
    void get_inv(int *Alpha,int len,int *Beta){
        //使用前清空Beta数组
        if(len==1){Beta[0]=fastpow(Alpha[0],mod-2);return ;}
        get_inv(Alpha,(len+1)>>1,Beta);init(len+len-1);
        for(int i=0;i<len;i++) g[i]=Alpha[i];for(int i=len;i<lim;i++) g[i]=0;
        NTT(Beta,1),NTT(g,1);
        for(int i=0;i<lim;i++) Beta[i]=mul(Beta[i],2-mul(Beta[i],g[i]));
        NTT(Beta,-1);
        for(int i=len;i<lim;i++) Beta[i]=0;
    }
    void convolution(int *Alpha,int la,int *Beta,int lb,int *Zeta){
        init(la+lb-1);for(int i=lb;i<lim;i++) Beta[i]=0;
        for(int i=la;i<lim;i++) Alpha[i]=0;
        NTT(Alpha,1),NTT(Beta,1);
        for(int i=0;i<lim;i++) Zeta[i]=mul(Alpha[i],Beta[i]);
        NTT(Zeta,-1);
        for(int i=la+lb-1;i<lim;i++) Zeta[i]=0;
    }
    const int inv2=fastpow(2,mod-2);
    int t[maxn];
    void get_sqrt(int *Alpha,int len,int *Beta){
        //使用前清空Beta数组
        if(len==1){Beta[0]=1;return ;}
        get_sqrt(Alpha,(len+1)>>1,Beta);get_inv(Beta,len,t);
        init(len+len-1);
        for(int i=0;i<len;i++) g[i]=Alpha[i];for(int i=len;i<lim;i++) g[i]=0;
        NTT(g,1),NTT(t,1),NTT(Beta,1);
        for(int i=0;i<lim;i++) Beta[i]=mul(inv2,(Beta[i]+mul(g[i],t[i]))%mod);
        NTT(Beta,-1);
        for(int i=len;i<lim;i++) Beta[i]=0;for(int i=0;i<lim;i++) t[i]=0;
    }
    void derivative(int *Alpha,int len,int *Beta){
        for(int i=1;i<len;i++) Beta[i-1]=mul(i,Alpha[i]);
        Beta[len-1] = 0;
    }
    void integral(int *Alpha,int len,int *Beta){
        for(int i=1;i<len;i++) Beta[i]=mul(Alpha[i-1],fastpow(i,mod-2));
        Beta[0] = 0;
    }
    int f1[maxn],f2[maxn],f3[maxn];
    void get_ln(int *Alpha,int len,int *Beta){
        memset(f2,0,sizeof(f2));
        derivative(Alpha,len,f1),get_inv(Alpha,len,f2);
        convolution(f1,len,f2,len,f3);
        integral(f3,len,Beta);
    }
    int e1[maxn],e2[maxn];
    void get_exp(int *Alpha,int len,int *Beta){
        //使用前清空Beta数组
        if(len==1) return void(Beta[0]=1);
        get_exp(Alpha,(len+1)>>1,Beta);get_ln(Beta,len,e1);
        init(len+len-1);
        for(int i=0;i<len;i++) e2[i]=Alpha[i];for(int i=len;i<lim;i++) e2[i]=0;
        NTT(e2,1),NTT(Beta,1),NTT(e1,1);
        for(int i=0;i<lim;i++) Beta[i]=mul(Beta[i],(1-e1[i]+e2[i]));
        NTT(Beta,-1);
        for(int i=len;i<lim;i++) Beta[i]=0;for(int i=0;i<lim;i++) e1[i]=0;
    }
    int g1[maxn],g2[maxn],g3[maxn];
    void division(int *Alpha,int la,int *Beta,int lb,int *Zeta,int *Gamma){
        if(la<lb){
            for(int i=0;i<la;i++) Gamma[i] = Alpha[i];
            return ;
        }
        for(int i=0;i<la;i++) g1[i]=Alpha[la-i-1];for(int i=0;i<lb;i++) g2[i]=Beta[lb-i-1];
        for(int i=lb;i<=la-lb;i++) g2[i]=0;
        Poly::init((la-lb+1)<<1|1);for(int i=0;i<lim;i++) g3[i]=0;
        get_inv(g2,la-lb+1,g3);
        convolution(g1,la,g3,la-lb+1,Zeta);
        for(int i=0;i<=la-lb;i++) if(i<la-lb-i) swap(Zeta[i],Zeta[la-lb-i]);else break;
        for(int i=la-lb+1;i<lim;i++) Zeta[i]=0;
        for(int i=0;i<lb;i++) g2[i]=Beta[i];for(int i=0;i<=la-lb;i++) g3[i]=Zeta[i];
        convolution(g2,lb,g3,la-lb+1,g1);
        for(int i=0;i<lb;i++) Gamma[i]=(Alpha[i]-g1[i])%mod;
    }
}

int n,Alpha[maxn],Beta[maxn],Gamma[maxn];
int g1[maxn],g2[maxn];

void Newton(int len,int *Alpha){
    if(len==1) return void(Alpha[0]=1);
    Newton((len+1)>>1,Alpha);
    for(int i=0;i<len;i++) Gamma[i]=0;
    for(int i=1;i<len;i++){
        int Okabe = fastpow(i,mod-2);
        for(int j=1;j*i<len;j++) Gamma[j*i]=(Gamma[j*i]+mul(Okabe,Alpha[j-1]))%mod;
    }
    //注意i从1开始，这里将g_0也顺手合并了进去
    memset(g1,0,sizeof(g1));
    Poly::get_ln(Alpha,len,g1);for(int i=0;i<len;i++) g1[i]=(g1[i]-Gamma[i])%mod;
    memset(g2,0,sizeof(g2));Poly::get_inv(Alpha,len,g2);
    g2[1]=(g2[1]-1)%mod;
    memset(Beta,0,sizeof(Beta));Poly::get_inv(g2,len,Beta);
    Poly::convolution(g1,len,Beta,len,g1);
    for(int i=0;i<len;i++) Alpha[i]=(Alpha[i]-g1[i])%mod;
}

int main(){
    scanf("%d",&n);Newton(n,Alpha);//这里得到的是除以x之后的结果
    int answer = Alpha[n-1];
    for(int i=(n>>1)+1;i<=n-1;i++) answer = (answer-mul(Alpha[i-1],Alpha[n-i-1]))%mod;
    if(!(n&1)) answer = answer-mul(fastpow(2,mod-2),mul(Alpha[n/2-1],Alpha[n/2-1]-1));
    answer%=mod;
    printf("%d\n",(answer+mod)%mod);
    return 0;
}

最后的最后，关于另一篇牛顿迭代题解的未解之谜2，我相信那是它写法上出了什么意外，因为我这里数组开到 6e5 就过了。