【LOJ #2320】「清华集训 2017」生成树计数(生成函数+分治NTT+多项式Exp)

传送门

参考的就是lojloj讨论里面jokerjoker写的做法

首先答案和度数有关,考虑树的pruferprufer序列
先不考虑那个idim\sum_{i}d_i^m
对每个连通块构造生成函数
ft(x)=iati+1(i+1)mxii!f_t(x)=\sum_{i}a_t^{i+1}(i+1)^m\frac{x^i}{i!}
f(x)=i(i+1)mxii!f(x)=\sum_{i}(i+1)^m\frac{x^i}{i!}
那么ans=(n2)!iaiif(aix)[xn2]ans=(n-2)!\prod_ia_i\prod_{i}f(a_ix)[x^{n-2}]
考虑求if(aix)\prod_if(a_ix)
ff取对数,最后再ExpExp回去
Ln(f)=ifixiLn(f)=\sum_if_ix^i
那么就是
if(aix)\prod_if(a_ix)
=Exp(ifij=0(aix)j)=Exp(\sum_{i}f_i\sum_{j=0}^{\infty}(a_ix)^j)

考虑求出t=[1,n],iait\forall t=[1,n],求出\sum_{i}a_i^t
=i11aix=iji(1ajx)i(1ajx)=\sum_{i}\frac{1}{1-a_ix}=\frac{\sum_{i}\prod_{j\not=i}(1-a_jx)}{\prod_{i}(1-a_jx)}
这个先分治NTTNTT求出s(x)=i(1aix)s(x)=\prod_i(1-a_ix)
考虑上面一块的意义,就是每次去掉一个单项式后乘加起来
那么s[xi]s[x^i]系数就是选ii个单项式乘起来之和
考虑对于分子那一块,xix^i的会删去被nin-i个没有去掉选中的单项式的情况算到
所以分子=i(ni)sixi=\sum_i(n-i)s_ix^i
然后和jj对应的位置乘起来即可

考虑还有一个idim\sum_{i}d_i^m
分配律拆开
那么就是对于每一个tt,将(j+1)m(j+1)^m变成(j+1)2m(j+1)^{2m}

g(x)=i(i+1)2mxii!,h(x)=g(x)f(x)g(x)=\sum_{i}(i+1)^{2m}\frac{x^i}{i!},h(x)=\frac{g(x)}{f(x)}

那么可以变形成(if(aix))(ih(aix))(\prod_{i}f(a_ix))(\sum_{i}h(a_ix))
后面又可以用刚才分治NTTNTT的方法求
然后乘起来即可

#include<bits/stdc++.h>
using namespace std;
#define cs const
#define re register
#define pb push_back
#define pii pair<int,int>
#define fi first
#define se second
#define ll long long
#define poly vector<int>
#define bg begin
cs int RLEN=1<<20|1;
inline char gc(){
	static char ibuf[RLEN],*ib,*ob;
	(ib==ob)&&(ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
	return (ib==ob)?EOF:*ib++;
}
inline int read(){
	char ch=gc();
	int res=0;bool f=1;
	while(!isdigit(ch))f^=ch=='-',ch=gc();
	while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
	return f?res:-res;
}
template<class tp>inline void chemx(tp &a,tp b){a<b?a=b:0;}
template<class tp>inline void chemn(tp &a,tp b){a>b?a=b:0;}
cs int mod=998244353,G=3;
inline int add(int a,int b){return (a+=b)>=mod?(a-mod):a;}
inline int dec(int a,int b){a-=b;return a+(a>>31&mod);}
inline int mul(int a,int b){static ll r;r=1ll*a*b;return r>=mod?(r%mod):r;}
inline void Add(int&a,int b){a+=b,a>=mod?(a-=mod):0;}
inline void Dec(int &a,int b){a-=b,a+=a>>31&mod;}
inline void Mul(int &a,int b){static ll r;r=1ll*a*b;a=(r>=mod)?(r%mod):r;}
inline int ksm(int a,int b,int res=1){for(;b;b>>=1,Mul(a,a))(b&1)&&(Mul(res,a),1);return res;}
inline int Inv(int x){return ksm(x,mod-2);}
cs int C=18;
poly w[C+1];
int rev[(1<<C)|5];
inline void init_w(){
	for(int i=1;i<=C;i++)w[i].resize(1<<(i-1));
	int wn=ksm(G,(mod-1)/(1<<C));
	w[C][0]=1;
	for(int i=1;i<(1<<(C-1));i++)w[C][i]=mul(w[C][i-1],wn);
	for(int i=C-1;i;i--)
	for(int j=0;j<(1<<(i-1));j++)w[i][j]=w[i+1][j<<1];
}
cs int N=30005;
int fac[N],ifac[N],inv[N];
inline void init_inv(){
	cs int len=N-5;
	fac[0]=ifac[0]=inv[0]=inv[1]=1;
	for(int i=1;i<=len;i++)fac[i]=mul(fac[i-1],i);
	ifac[len]=Inv(fac[len]);
	for(int i=len-1;i;i--)ifac[i]=mul(ifac[i+1],i+1);
	for(int i=2;i<=len;i++)inv[i]=mul(mod-mod/i,inv[mod%i]);
}
inline void init_rev(int lim){
	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)*(lim>>1));
}
inline void ntt(poly &f,int lim,int kd){
	for(int i=0;i<lim;i++)if(i>rev[i])swap(f[i],f[rev[i]]);
	for(int mid=1,l=1,a0,a1;mid<lim;mid<<=1,l++)
	for(int i=0;i<lim;i+=mid<<1)
	for(int j=0;j<mid;j++)
	a0=f[i+j],a1=mul(w[l][j],f[i+j+mid]),f[i+j]=add(a0,a1),f[i+j+mid]=dec(a0,a1);
	if(kd==-1){
		reverse(f.bg()+1,f.bg()+lim);
		for(int i=0,inv=Inv(lim);i<lim;i++)Mul(f[i],inv);
	}
}
inline poly operator *(poly a,poly b){
	int deg=a.size()+b.size()-1,lim=1;
	if(deg<=32){
		poly c(deg,0);
		for(int i=0;i<a.size();i++)
		for(int j=0;j<b.size();j++)
		Add(c[i+j],mul(a[i],b[j]));
		return c;
	}
	while(lim<deg)lim<<=1;
	init_rev(lim);
	a.resize(lim),ntt(a,lim,1);
	b.resize(lim),ntt(b,lim,1);
	for(int i=0;i<lim;i++)Mul(a[i],b[i]);
	ntt(a,lim,-1),a.resize(deg);
	return a;
}
inline poly Inv(poly a,int deg){
	poly b(1,Inv(a[0])),c;
	int n=a.size();
	for(int lim=4;lim<(deg<<2);lim<<=1){
		c.resize(lim>>1);
		for(int i=0;i<(lim>>1);i++)c[i]=(i<n)?a[i]:0;
		c.resize(lim),b.resize(lim);
		init_rev(lim);
		ntt(b,lim,1),ntt(c,lim,1);
		for(int i=0;i<lim;i++)Mul(b[i],dec(2,mul(b[i],c[i])));
		ntt(b,lim,-1),b.resize(lim>>1);
	}
	b.resize(deg);return b;
}
inline poly deriv(poly a){
	for(int i=0;i+1<a.size();i++)a[i]=mul(a[i+1],i+1);
	a.pop_back();return a;
}
inline poly integ(poly a){
	a.pb(0);
	for(int i=(int)a.size()-1;i;i--)a[i]=mul(a[i-1],inv[i]);
	a[0]=0;return a;
}
inline poly operator +(poly a,poly b){
	a.resize(max(a.size(),b.size()));
	for(int i=0;i<b.size();i++)Add(a[i],b[i]);return a;
}
inline void Mult(poly &f,poly &g,int lim){
	for(int i=0;i<lim;i++)Mul(f[i],g[i]);
}
inline poly Ln(poly a,int deg){
	a=integ(Inv(a,deg)*deriv(a)),a.resize(deg);return  a;
}
inline poly Exp(poly a,int deg){
	poly b(1,1),c;int n=a.size();
	for(int lim=2;lim<(deg<<1);lim<<=1){
		c=Ln(b,lim);
		for(int i=0;i<lim;i++)c[i]=dec(i<n?a[i]:0,c[i]);
		c[0]++,b=b*c,b.resize(lim);
	}
	b.resize(deg);return b;
}
poly s2[N<<2];
#define lc (u<<1)
#define rc ((u<<1)|1)
#define mid ((l+r)>>1)
void build(int u,int l,int r,int *v){
	if(l==r){s2[u].resize(2),s2[u][0]=1,s2[u][1]=dec(0,v[l]);return;}
	build(lc,l,mid,v),build(rc,mid+1,r,v);
	s2[u]=s2[lc]*s2[rc];
}
#undef lc
#undef rc
#undef mid
int ans,n,m,a[N];
int main(){
	#ifdef Stargazer
	freopen("lx.cpp","r",stdin);
	#endif
	init_w(),init_inv(),n=read(),m=read(),ans=fac[n-2];
	for(int i=1;i<=n;i++)a[i]=read(),Mul(ans,a[i]);
	poly f(n-1,0),f0(n-1,0);
	for(int i=0,t;i<=n-2;i++)t=ksm(i+1,m),f[i]=mul(t,ifac[i]),f0[i]=mul(mul(t,t),ifac[i]);
	poly g=f0*Inv(f,n-1);f=Ln(f,n-1);
//	for(int i=0;i<f0.size();i++)cout<<g[i]<<" ";puts("");
	build(1,1,n,a);
	poly t(n-1,0);
	for(int i=0;i<n-1;i++)t[i]=mul(s2[1][i],n-i);
	poly h=Inv(s2[1],n+1)*t;h.resize(n-1);
	Mult(g,h,n-1),Mult(f,h,n-1);
	f=Exp(f,n-1)*g;
	Mul(ans,f[n-2]);
	cout<<ans;
}
posted @ 2019-11-07 23:38  Stargazer_cykoi  阅读(157)  评论(0编辑  收藏  举报