[清华集训2017]生成树计数

一个用斯特林数的(O(nmlog^2n))做法

考虑我们要求的答案即

[sum_{T}(prod_{i=1}^n d_i^{m})(sum_{i=1}^nd_{i}^m)(prod_{i=1}^na_i^{d_i}) ]

至于为什么有这个(prod_{i=1}^na_i^{d_i})，是因为在第(i)个联通块里连的这(d_i)条边的端点可以是(a_i)个点中随便选，于是就是(a_i^{d_i})

考虑暴力化柿子

[Ans=sum_Tsum_{i=1}^nd_i^{m}prod_{j=1}^nd_j^ma_j^{d_j}=sum_Tsum_{i=1}^nd_i^{2m} imes a_i^{d_i}prod_{i eq j}d_j^ma_j^{d_j} ]

其实就是从后面那个(prod)里拆了个(d_i^ma_i^{d_i})出来

这老写(sum_T)多难看，把生成树转化成( m prufer)序列，设(c_i)表示第(i)个点在( m prufer)序列中出现的次数，显然有(d_i=c_i+1)，有

[Ans=sum_{sum c=n-2}sum_{i=1}^n(c_i+1)^{2m}a_i^{c_i+1}prod_{i eq j}(c_j+1)^ma_j^{c_j+1} ]

({c_i})对应的( m prufer)序列共有(frac{(n-2)!}{prod c_i!})，于是我们能进一步写成

[frac{(n-2)!}{prod c_i!}sum_{i=1}^n(c_i+1)^{2m}a_i^{c_i+1}prod_{i eq j}(c_j+1)^ma_j^{c_j+1} ]

不妨将(frac{1}{prod c_i!})分配进去，则有

[(n-2)!sum_{i=1}^nfrac{(c_i+1)^{2m}a_i^{c_i+1}}{c_i!}prod_{i eq j}frac{(c_j+1)^ma_j^{c_j+1}}{c_j!} ]

不难设两个( m EGF)，({ m A}_i(x)=sum_{j=0}^{infty}frac{(j+1)^{2m}a_i^{j+1}}{j!}x^j,{ m B}_i(x)=sum_{j=0}^{infty}frac{(j+1)^{m}a_i^{j+1}}{j!}x^j)

这样我们的答案就能写成

[(n-2)![n-2](sum_{i=1}^n{ m A}_i(x) imes prod_{i eq j}{ m B}_j(x)) ]

即那些个( m EGF)卷起来的第(n-2)次项系数

现在就是神奇的魔法了，以({ m B}_i(x))为例，我们先积分一下，得到的是

[sum_{j=0}^{infty}frac{(j+1)^ma_i^{j+1}}{(j+1)!}x^{j+1}=sum_{j=1}^{infty}frac{j^ma_i^{j}}{j!}x^{j} ]

用第二类斯特林数将(j^m)展开

[sum_{j=1}^{infty}frac{a_i^{j}}{j!}x^{j}sum_{k=0}^megin{Bmatrix}m\kend{Bmatrix}j^{underline k}=sum_{k=0}^megin{Bmatrix}m\kend{Bmatrix}sum_{j=1}^{infty}frac{a_i^{j}j^{underline k}}{j!}x^{j} ]

注意到(frac{j^{underline k}}{j!})其实就是(frac{!}{(j-k)!})，于是我们从后面提一个(a_i^{k}x_k)出来，并将后面的柿子泰勒展开一下

[sum_{k=0}^megin{Bmatrix}m\kend{Bmatrix}a_i^{k}x^ksum_{j=1}^{infty}frac{a_i^{j-k}}{(j-k)!}x^{j-k}=sum_{k=0}^megin{Bmatrix}m\kend{Bmatrix}a_i^{k}x^ke^{a_ix} ]

我们对于(int { m B}_i(x) m dx)得到了一个优美的形式，那么我们再求导回去，我们可以看成是(sum_{k=0}^megin{Bmatrix}m\kend{Bmatrix}a_i^{k}x^k)乘(e^{a_ix})这两个函数相乘求导，那么根据求导的法则([f(x)g(x)]'=f(x)g'(x)+f'(x)g(x))，就有

[sum_{k=0}^megin{Bmatrix}m\kend{Bmatrix}a_i^{k}kx^{k-1}e^{a_ix}+sum_{k=0}^megin{Bmatrix}m\kend{Bmatrix}a_i^{k}x^ka_ie^{a_ix} ]

整理一下就有

[e^{a_ix}(sum_{k=1}^m(k+1)egin{Bmatrix}m\k+1end{Bmatrix}a_i^{k+1}x^k+sum_{k=0}^megin{Bmatrix}m\kend{Bmatrix}a_i^{k+1}x^k) ]

根据第二类斯特林数的递推式，我们惊奇的发现((k+1)egin{Bmatrix}m\k+1end{Bmatrix}+egin{Bmatrix}m\kend{Bmatrix}=egin{Bmatrix}m+1\k+1end{Bmatrix})

于是我们得到一个极其优美的形式

[{ m B}_i(x)=e^{a_ix}sum_{k=0}^megin{Bmatrix}m+1\k+1end{Bmatrix}a_i^{k+1}x^k ]

类似的

[{ m A}_i(x)=e^{a_ix}sum_{k=0}^{2m}egin{Bmatrix}2m+1\k+1end{Bmatrix}a_i^{k+1}x^k ]

再设({alpha_i(x)}=frac{{ m A}_i(x)}{e^{a_ix}},{eta _i(x)}=frac{{ m B}_i(x)}{e^{a_ix}})

回到最开始的式子，不妨将(sum_{i=1}^n{ m A}_i(x) imes prod_{i eq j}{ m B}_j(x))写成(e^{xsum a}sum_{i=1}^n{alpha}_i(x) imes prod_{i eq j}{eta}_j(x))，之后我们就可以进一步写成

[prod_{i=1}^n{eta}_i(x) imes sum_{i=1}^nfrac{{alpha}_i(x)}{{eta}_i(x)} ]

(sum_{i=1}^nfrac{{alpha}_i(x)}{{eta}_i(x)})并没有想到什么好做法，那就一边分治一边通分吧，之后分治通分到最后发现分母就是(prod_{i=1}^n{eta}_i(x))，于是我们只需要分子就好了

代码鸽了，有空来补

UPD on 2019.12.28 原来的式子最后一步钦定错了，好像导致一位老鸽调了一天，手动谢罪

来补上代码了

#include<bits/stdc++.h>
#define re register
const int mod=998244353;
const int G[2]={3,(mod+1)/3};
inline int read() {
	char c=getchar();int x=0;while(c<'0'||c>'9')c=getchar();
	while(c>='0'&&c<='9')x=(x<<3)+(x<<1)+c-48,c=getchar();return x;
}
const int Maxlen=65539;
const int maxn=3e4+5;
std::vector<int> p[maxn*3],q[maxn*3];
int n,m,a[maxn],S[105][105],sum,len,inv[maxn],__og[2][50];
int rev[Maxlen],A[Maxlen],C[Maxlen],B[Maxlen],D[Maxlen];
inline int dqm(int x) {return x<0?x+mod:x;}
inline int qm(int x) {return x>=mod?x-mod:x;}
inline int ksm(int a,int b) {
	int S=1;for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)S=1ll*S*a%mod;return S;
}
inline void NTT(int *f,int o) {
	for(re int i=0;i<len;i++)if(i<rev[i])std::swap(f[i],f[rev[i]]);
	for(re int w=0,i=2;i<=len;i<<=1,++w) {
		int ln=i>>1,og1;
		if(!__og[o][w])og1=__og[o][w]=ksm(G[o],(mod-1)/i);else og1=__og[o][w];
		for(re int t,og=1,l=0;l<len;l+=i,og=1)
			for(re int x=l;x<l+ln;++x) {
				t=1ll*f[x+ln]*og%mod,og=1ll*og1*og%mod;
				f[x+ln]=dqm(f[x]-t);f[x]=qm(f[x]+t);
			}
	}
	if(!o)return;int Inv=ksm(len,mod-2);for(re int i=0;i<len;i++)f[i]=1ll*f[i]*Inv%mod;
}
void cdq(int l,int r,int k) {
	if(l==r) {
		int nw=a[l];
		for(re int i=0;i<=m+m;i++)p[k].push_back(1ll*nw*S[m+m+1][i+1]%mod),nw=1ll*nw*a[l]%mod;
		nw=a[l];for(re int i=0;i<=m;i++)q[k].push_back(1ll*nw*S[m+1][i+1]%mod),nw=1ll*nw*a[l]%mod;
		return;
	}
	int mid=l+r>>1;
	cdq(l,mid,k<<1),cdq(mid+1,r,k<<1|1);
	int lenl=p[k<<1].size(),lenr=p[k<<1|1].size();
	len=1;while(len<=lenl+lenr+1) len<<=1;
	for(re int i=0;i<len;i++)rev[i]=rev[i>>1]>>1|((i&1)?len>>1:0);
	for(re int i=0;i<p[k<<1].size();++i) A[i]=p[k<<1][i];
	for(re int i=0;i<p[k<<1|1].size();++i) B[i]=p[k<<1|1][i];
	for(re int i=0;i<q[k<<1].size();++i) C[i]=q[k<<1][i];
	for(re int i=0;i<q[k<<1|1].size();++i) D[i]=q[k<<1|1][i];
	NTT(A,0),NTT(B,0),NTT(C,0),NTT(D,0);
	for(re int i=0;i<len;i++) A[i]=qm(1ll*A[i]*D[i]%mod+1ll*B[i]*C[i]%mod),C[i]=1ll*C[i]*D[i]%mod;
	NTT(A,1),NTT(C,1);
	for(re int i=0;i<len;i++) {
		if(i>n-2) break;
		p[k].push_back(A[i]);
		q[k].push_back(C[i]);
	}
	for(re int i=0;i<len;i++) A[i]=B[i]=C[i]=D[i]=0;
	int lp=p[k].size(),lq=q[k].size();lp--;lq--;
	for(re int i=lp;i>=0;--i) 
		if(!p[k][i]) p[k].pop_back();else break;
	for(re int i=lq;i>=0;--i)
		if(!q[k][i]) q[k].pop_back();else break;
}
signed main() {
	n=read(),m=read();S[0][0]=1;inv[0]=1;
	for(re int i=1;i<=m+m+1;i++)
		for(re int j=1;j<=i;++j) S[i][j]=qm(S[i-1][j-1]+1ll*S[i-1][j]*j%mod);
	for(re int i=1;i<=n;i++) a[i]=read(),sum=qm(sum+a[i]);
	if(n==1) return puts("1"),0;
	inv[1]=1;for(re int i=2;i<=n;i++) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
	cdq(1,n,1);for(re int i=0;i<p[1].size();++i) A[i]=p[1][i];int ans=0;
	for(re int i=1;i<=n-2;i++) inv[i]=1ll*inv[i]*inv[i-1]%mod;
	for(re int nw=1,i=0;i<=n-2;i++,nw=1ll*nw*sum%mod)B[i]=1ll*inv[i]*nw%mod;
	for(re int i=0;i<=n-2;i++)ans=qm(ans+1ll*A[i]*B[n-2-i]%mod);
	for(re int i=1;i<=n-2;i++)ans=1ll*ans*i%mod;
	std::cout<<ans<<std::endl;return 0;
}
//g++ lg4002.cpp -o lg4002