[清华集训2017]生成树计数
一个用斯特林数的\(O(nm\log^2n)\)做法
考虑我们要求的答案即
至于为什么有这个\(\prod_{i=1}^na_i^{d_i}\)，是因为在第\(i\)个联通块里连的这\(d_i\)条边的端点可以是\(a_i\)个点中随便选，于是就是\(a_i^{d_i}\)
考虑暴力化柿子
其实就是从后面那个\(\prod\)里拆了个\(d_i^ma_i^{d_i}\)出来
这老写\(\sum_T\)多难看，把生成树转化成\(\rm prufer\)序列，设\(c_i\)表示第\(i\)个点在\(\rm prufer\)序列中出现的次数，显然有\(d_i=c_i+1\)，有
\(\{c_i\}\)对应的\(\rm prufer\)序列共有\(\frac{(n-2)!}{\prod c_i!}\)，于是我们能进一步写成
不妨将\(\frac{1}{\prod c_i!}\)分配进去，则有
不难设两个\(\rm EGF\)，\({\rm A}_i(x)=\sum_{j=0}^{\infty}\frac{(j+1)^{2m}a_i^{j+1}}{j!}x^j,{\rm B}_i(x)=\sum_{j=0}^{\infty}\frac{(j+1)^{m}a_i^{j+1}}{j!}x^j\)
这样我们的答案就能写成
即那些个\(\rm EGF\)卷起来的第\(n-2\)次项系数
现在就是神奇的魔法了，以\({\rm B}_i(x)\)为例，我们先积分一下，得到的是
用第二类斯特林数将\(j^m\)展开
注意到\(\frac{j^{\underline k}}{j!}\)其实就是\(\frac{!}{(j-k)!}\)，于是我们从后面提一个\(a_i^{k}x_k\)出来，并将后面的柿子泰勒展开一下
我们对于\(\int {\rm B}_i(x)\rm dx\)得到了一个优美的形式，那么我们再求导回去，我们可以看成是\(\sum_{k=0}^m\begin{Bmatrix}m\\k\end{Bmatrix}a_i^{k}x^k\)乘\(e^{a_ix}\)这两个函数相乘求导，那么根据求导的法则\([f(x)g(x)]'=f(x)g'(x)+f'(x)g(x)\)，就有
整理一下就有
根据第二类斯特林数的递推式，我们惊奇的发现\((k+1)\begin{Bmatrix}m\\k+1\end{Bmatrix}+\begin{Bmatrix}m\\k\end{Bmatrix}=\begin{Bmatrix}m+1\\k+1\end{Bmatrix}\)
于是我们得到一个极其优美的形式
类似的
再设\({\alpha_i(x)}=\frac{{\rm A}_i(x)}{e^{a_ix}},{\beta _i(x)}=\frac{{\rm B}_i(x)}{e^{a_ix}}\)
回到最开始的式子，不妨将\(\sum_{i=1}^n{\rm A}_i(x)\times \prod_{i\neq j}{\rm B}_j(x)\)写成\(e^{x\sum a}\sum_{i=1}^n{\alpha}_i(x)\times \prod_{i\neq j}{\beta}_j(x)\)，之后我们就可以进一步写成
\(\sum_{i=1}^n\frac{{\alpha}_i(x)}{{\beta}_i(x)}\)并没有想到什么好做法，那就一边分治一边通分吧，之后分治通分到最后发现分母就是\(\prod_{i=1}^n{\beta}_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
