题解 P5900 【无标号无根树计数】
Rui_R
2020-12-21 20:01:52
题意：题如其名。
[原题](https://www.luogu.com.cn/problem/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$ 就过了。