【BZOJ 3684】大朋友和多叉树

传送门

---

problem

我们的大朋友很喜欢计算机科学，而且尤其喜欢多叉树。

对于一棵带有正整数点权的有根多叉树，如果它满足这样的性质，我们的大朋友就会将其称作神犇：

1. 点权为 $1$ 的结点是叶子结点；
2. 对于任一点权大于 $1$ 的结点 $u$，$u$ 的孩子数目 $de{g}_u$ 属于集合 $D$，且 $u$ 的点权等于这些孩子结点的点权之和。

给出一个整数 $n$，请求出根节点权值为 $n$ 的神犇多叉树的个数对 $950009857$（$453\times 2^{21}+1$，一个质数）取模后的值。

数据范围：$1\le m<n\le 10^5$（$m$ 为 $D$ 集合大小），$2\le D_i\le n$。

---

solution

设 $D$ 的生成函数是 $C(x)$，那么有 $C(x)={\sum }_{i\in D}{x}^i$。

又设树的生成函数为 $F(x)$，有：

$$F(x)=x+\sum _{i\in D}{F}^i(x)$$

就相当于是枚举儿子数量 $i$，有 $i$ 个儿子就有 $F^i(x)$ 的贡献，单出来的那个 $x$ 是单点成树的情况。

那么，我们现在要求的是 $[x^n]F(x)$（就是 $F(x)$ 中 $x^n$ 这一项的系数）。

拉格朗日反演

如果多项式 $f$ 和 $g$ 满足 $g(f(x))=x$，那么一定有 $f(g(x))=x$，这时我们称 $f$ 和 $g$ 互为复合逆。

若 $f$ 和 $g$ 互为复合逆，我们有：

$$[x^n]f(x)=\frac{1}{n}[x^{-1}]\frac{1}{g^n(x)}$$

---

然后回到这道题来。

看一下 $C(x)$ 和 $F(x)$ 的式子，容易得到 $F(x)=x+C(F(x))$，即 $F(x)-C(F(x))=x$。

令 $G(x)=x-C(x)$ 的话，发现上式就是 $G(F(x))=x$，那就可以拉格朗日反演了：

$$[x^n]F(x)=\frac{1}{n}[x^{-1}]\frac{1}{G^n(x)}$$

看到 $[x^{-1}]$，$x^{-1}$ 项的系数？不妨对右边乘一个 $x^n$，然后变成：

$$\begin{array}{rl}[x^n]F(x)& =\frac{1}{n}[x^{n-1}]\frac{x^n}{G^n(x)}\\ & =\frac{1}{n}[x^{n-1}](\frac{x}{G(x)}{)}^n\end{array}$$

根据 $G(x)=x-C(x)$，那么 $\frac{x}{G(x)}=\frac{1}{1-{\sum }_{i\in D}{x}^{i-1}}$。剩下的就比较简单了，直接多项式快速幂+求逆就可以了。

---

code

#include<cstdio>
#include<vector>
#include<cstring>
#include<algorithm>
#define N 1000005
#define P 950009857
using namespace std;
const int g=7;
typedef vector<int> poly;
int n,m,pos[N],inv[N];
int add(int x,int y)  {return x+y>=P?x+y-P:x+y;}
int dec(int x,int y)  {return x-y< 0?x-y+P:x-y;}
int mul(int x,int y)  {return 1ll*x*y%P;}
int power(int a,int b,int ans=1){
	for(;b;b>>=1,a=mul(a,a))
		if(b&1)  ans=mul(ans,a);
	return ans;
}
int C=21,*w[22];
void prework(){
	for(int i=1;i<=C;++i)
		w[i]=new int[1<<(i-1)];
	int now=power(g,(P-1)/(1<<C));
	w[C][0]=1;
	for(int i=1;i<(1<<(C-1));++i)  w[C][i]=mul(w[C][i-1],now);
	for(int i=C-1;i;--i)
		for(int j=0;j<(1<<(i-1));++j)
			w[i][j]=w[i+1][j<<1];
}
void init_inv(){
	inv[0]=inv[1]=1;
	for(int i=2;i<N;++i)  inv[i]=mul(P-P/i,inv[P%i]);
}
void init_pos(int limit){
	for(int i=0;i<limit;++i)  pos[i]=(pos[i>>1]>>1)|((i&1)*(limit>>1));
}
void NTT(poly &f,int limit,int type){
	for(int i=0;i<limit;++i)
		if(pos[i]>i)  swap(f[i],f[pos[i]]);
	for(int mid=1,l=1;mid<limit;mid<<=1,++l){
		for(int i=0;i<limit;i+=(mid<<1)){
			for(int j=0;j<mid;++j){
				int p0=f[i+j],p1=mul(f[i+j+mid],w[l][j]);
				f[i+j]=add(p0,p1),f[i+j+mid]=dec(p0,p1);
			}
		}
	}
	if(type==-1&&(reverse(f.begin()+1,f.begin()+limit),1)){
		int inv=power(limit,P-2);
		for(int i=0;i<limit;++i)  f[i]=mul(f[i],inv);
	}
}
poly operator*(poly A,int B){
	for(int i=0;i<A.size();++i)  A[i]=mul(A[i],B);
	return A;
}
poly operator*(poly A,poly B){
	int lim=1,len=A.size()+B.size()-2;
	while(lim<=len)  lim<<=1;init_pos(lim);
	A.resize(lim),NTT(A,lim,1);
	B.resize(lim),NTT(B,lim,1);
	for(int i=0;i<lim;++i)  A[i]=mul(A[i],B[i]);
	NTT(A,lim,-1),A.resize(len+1);
	return A;
}
poly Inv(poly A,int len){
	poly C,B(1,power(A[0],P-2));
	for(int lim=4;lim<(len<<2);lim<<=1){
		init_pos(lim);
		C=A,C.resize(lim>>1);
		C.resize(lim),NTT(C,lim,1);
		B.resize(lim),NTT(B,lim,1);
		for(int i=0;i<lim;++i)  B[i]=mul(B[i],dec(2,mul(B[i],C[i])));
		NTT(B,lim,-1),B.resize(lim>>1);
	}
	B.resize(len);return B;
}
poly Deriv(poly A){
	for(int i=0;i<A.size()-1;++i)  A[i]=mul(A[i+1],i+1);
	A.pop_back();return A;
}
poly Integ(poly A){
	A.push_back(0);
	for(int i=A.size()-1;i;--i)  A[i]=mul(A[i-1],inv[i]);
	A[0]=0;return A;
}
poly Ln(poly A,int len){
	A=Integ(Deriv(A)*Inv(A,len)),A.resize(len);
	return A;
}
poly Exp(poly A,int len){
	poly B(1,1),C;A.resize(len<<1);
	for(int lim=2;lim<(len<<1);lim<<=1){
		C=Ln(B,lim);
		for(int i=0;i<lim;++i)  C[i]=dec(A[i],C[i]);
		C[0]=add(C[0],1),B=B*C;
		B.resize(lim);
	}
	B.resize(len);return B;
}
poly Power(poly A,int len,int k){
	A=Exp(Ln(A,len)*k,len),A.resize(len);
	return A;
}
poly G;
int main(){
	scanf("%d%d",&n,&m);
	prework(),init_inv();
	G.resize(n);
	for(int i=1,x;i<=m;++i)  scanf("%d",&x),G[x-1]=P-1;
	G[0]++,G=Power(Inv(G,n),n,n);
	printf("%d\n",mul(G[n-1],power(n,P-2)));
	return 0;
}