「BZOJ 3684」大朋友和多叉树

BZOJ 3684: 大朋友和多叉树

题意

求包含 \(s\) 个叶子，非叶子节点的孩子数目在集合 \(D\) 中的有根树数量。

孩子之间有顺序，保证 \(1\notin D\)。

模质数 \(950009857=453\times 2^{21}+1\)

\(s,|D|\le 10^5\)

做法

好像就是板题

令 \(D(x)=\sum\limits_{i=2}^\infty [i\in D]x^i\)

可以发现答案的生成函数 \(f(x)=D(f(x))+x\)，\(f(x)\) 为 \(x-D(x)\) 的复合逆

所求即为 \([x^s]f(x)\)

根据拉格朗日反演公式

\[ [x^s]f(x) = \frac{1}{s} [x^{s-1}]\left(\frac{x-D(x)}{x}\right)^{-s} \]

可以在 \(\mathcal O(n\log n)\) 的复杂度内完成

代码

#include<cstdio>
#include<algorithm>
#include<cctype>
#include<string.h>
#include<cmath>
#include<vector>

using namespace std;
#define ull unsigned long long

inline char read() {
	static const int IN_LEN = 1000000;
	static char buf[IN_LEN], *s, *t;
	return (s==t?t=(s=buf)+fread(buf,1,IN_LEN,stdin),(s==t?-1:*s++):*s++);
}
template<class T>
inline void read(T &x) {
	static bool iosig;
	static char c;
	for (iosig=false, c=read(); !isdigit(c); c=read()) {
		if (c == '-') iosig=true;
		if (c == -1) return;
	}
	for (x=0; isdigit(c); c=read()) x=((x+(x<<2))<<1)+(c^'0');
	if (iosig) x=-x;
}

const unsigned N = 1<<18, P = 950009857;
struct Z{
	unsigned x;
	Z(const unsigned _x=0):x(_x){}
	inline Z operator +(const Z &rhs)const{ return x+rhs.x<P?x+rhs.x:x+rhs.x-P;}
	inline Z operator -(const Z &rhs)const{ return x<rhs.x?x-rhs.x+P:x-rhs.x;}
	inline Z operator -()const{ return x?P-x:0;}
	inline Z operator *(const Z &rhs)const{ return static_cast<ull>(x)*rhs.x%P;}
	inline Z operator +=(const Z &rhs){ return x=x+rhs.x<P?x+rhs.x:x+rhs.x-P, *this;}
	inline Z operator -=(const Z &rhs){ return x=x<rhs.x?x-rhs.x+P:x-rhs.x, *this;}
	inline Z operator *=(const Z &rhs){ return x=static_cast<ull>(x)*rhs.x%P, *this;}
};
int s, m;
vector<Z> d;

namespace Poly{
	Z w[N];// for DFT
	Z Inv[N];// for Integral

	inline Z Pow(Z x, int y=P-2){
		Z ans=1;
		for(; y; y>>=1, x=x*x) if(y&1) ans=ans*x;
		return ans;
	}
	inline void Init(){
		for(unsigned i=1; i<N; i<<=1){
			w[i]=1;
			Z t=Pow(7, (P-1)/i/2);
			for(unsigned j=1; j<i; ++j) w[i+j]=w[i+j-1]*t;
		}
		Inv[1]=1;
		for(unsigned i=2; i<N; ++i) Inv[i]=Inv[P%i]*(P-P/i);
	}
	inline int Get(int x){ int n=1; while(n<=x) n<<=1; return n;}
	inline void DFT(vector<Z> &f, int n){
		static ull F[N];
		if((int)f.size()!=n) f.resize(n);
		for(int i=0, j=0; i<n; ++i){
			F[i]=f[j].x;
			for(int k=n>>1; (j^=k)<k; k>>=1);
		}
		for(int i=1; i<n; i<<=1) for(int j=0; j<n; j+=i<<1){
			Z *W=w+i;
			ull *F0=F+j, *F1=F+j+i;
			for(int k=j; k<j+i; ++k, ++W, ++F0, ++F1){
				ull t=(*F1)*(W->x)%P;
				(*F1)=*F0+P-t, (*F0)+=t;
			}
		}
		for(int i=0; i<n; ++i) f[i]=F[i]%P;
	}
	inline void IDFT(vector<Z> &f, int n){
		f.resize(n), reverse(f.begin()+1, f.end());
		DFT(f, n);
		Z I=Pow(n);
		for(int i=0; i<n; ++i) f[i]=f[i]*I;
	}
	inline vector<Z> operator +(const vector<Z> &f, const vector<Z> &g){
		vector<Z> ans=f;
		for(unsigned i=0; i<f.size(); ++i) ans[i]+=g[i];
		return ans;
	}
	inline vector<Z> operator *(const vector<Z> &f, const vector<Z> &g){
		if((ull)f.size()*g.size()<=1000){
			vector<Z> ans;
			ans.resize(f.size()+g.size()-1);
			for(unsigned i=0; i<f.size(); ++i) for(unsigned j=0; j<g.size(); ++j)
				ans[i+j]+=f[i]*g[j];
			return ans;
		}
		static vector<Z> F, G;
		F=f, G=g;
		int p=Get(f.size()+g.size()-2);
		DFT(F, p), DFT(G, p);
		for(int i=0; i<p; ++i) F[i]*=G[i];
		IDFT(F, p);
		return F.resize(f.size()+g.size()-1), F;
	}
	vector<Z> &PolyInv(const vector<Z> &f, int n=-1){
		if(n==-1) n=f.size();
		if(n==1){
			static vector<Z> ans;
			return ans.clear(), ans.push_back(Pow(f[0])), ans;
		}
		vector<Z> &ans=PolyInv(f, (n+1)/2);
		vector<Z> tmp(&f[0], &f[0]+n);
		int p=Get(n*2-2);
		DFT(tmp, p), DFT(ans, p);
		for(int i=0; i<p; ++i) ans[i]=((Z)2-ans[i]*tmp[i])*ans[i];
		IDFT(ans, p);
		return ans.resize(n), ans;
	}
	// a=d*b+r
	inline void PolyDiv(const vector<Z> &a, const vector<Z> &b, vector<Z> &d, vector<Z> &r){
		if(b.size()>a.size()) return d.clear(), (void)(r=a);

		vector<Z> A=a, B=b, iB;
		int n=a.size(), m=b.size();
		reverse(A.begin(), A.end()), reverse(B.begin(), B.end());
		B.resize(n-m+1), iB=PolyInv(B, n-m+1);
		d=A*iB;
		d.resize(n-m+1), reverse(d.begin(), d.end());

		r=b*d, r.resize(m-1);
		for(int i=0; i<m-1; ++i) r[i]=a[i]-r[i];
	}
	inline vector<Z> Derivative(const vector<Z> &a){
		vector<Z> ans(a.size()-1);
		for(unsigned i=1; i<a.size(); ++i) ans[i-1]=a[i]*i;
		return ans;
	}
	inline vector<Z> Integral(const vector<Z> &a){
		vector<Z> ans(a.size()+1);
		for(unsigned i=0; i<a.size(); ++i) ans[i+1]=a[i]*Inv[i+1];
		return ans;
	}
	inline vector<Z> PolyLn(const vector<Z> &f){
		vector<Z> ans=Derivative(f)*PolyInv(f);
		ans.resize(f.size()-1);
		return Integral(ans);
	}
	vector<Z> PolyExp(const vector<Z> &f, int n=-1){
		if(n==-1) n=f.size();
		if(n==1) return {1};
		vector<Z> ans=PolyExp(f, (n+1)/2), tmp;
		ans.resize(n), tmp=PolyLn(ans);
		for(Z &i:tmp) i=-i;
		++tmp[0].x;
		ans=ans*(tmp+f);
		return ans.resize(n), ans;
	}
	inline vector<Z> PolyPower(const vector<Z> &f, int k){
		vector<Z> ans=PolyLn(f);
		for(Z &i:ans) i=i*k;
		return PolyExp(ans);
	}
}
int main() {
	Poly::Init();
	read(s), read(m), d.resize(s);
	d[0]=1;
	for(int i=1, x; i<=m; ++i) read(x), d[x-1]+=P-1;
	d=Poly::PolyPower(d, P-s);
	return printf("%d", (d[s-1]*Poly::Pow(s)).x), 0;
}