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

Posted 2021-03-12 cjoiershiina-mashiro

Link
下文中的点指的是题目给的(n)个连通块。
因为题目给的式子只与度数和点权相关，因此考虑Prufer序列。
枚举每个点在Prufer序列中的出现次数，那么此时的贡献就是Prufer序列的个数乘上该生成树的价值。
(ans=(n-2)!sumlimits_{sum d_i=n-2}prodlimits_{i=1}^nfrac{a_i^{d_i+1}}{d_i!}(d_i+1)^msumlimits_{i=1}^n(d_i+1)^m)
(ans=(n-2)!sumlimits_{i=1}^na_iprodlimits_{sum d_i=n-2}prodlimits_{i=1}^nfrac{a_i^{d_i}}{d_i!}(d_i+1)^msumlimits_{i=1}^n(d_i+1)^m)
(ans=(n-2)!prodlimits_{i=1}^na_isumlimits_{sum d_i=n-2}sumlimits_{i=1}^nfrac{a_i^{d_i}}{d_i!}(d_i+1)^{2m}prodlimits_{j e i}frac{a_j^{d_j}}{d_j!}(d_j+1)^m)
前面的((n-2)!prodlimits_{i=1}^na_i)可以直接算，那么我们要求的就是(sumlimits_{sum d_i=n-2}sumlimits_{i=1}^nfrac{a_i^{d_i}}{d_i!}(d_i+1)^{2m}prodlimits_{j e i}frac{a_j^{d_j}}{d_j!}(d_j+1)^m)，考虑以(d_i)为下标构造生成函数。
设(P(x)=sumfrac{(i+1)^m}{i!}x^i,Q(x)=sumfrac{(i+1)^{2m}}{i!}x^i,F(x)=sumlimits_{i=1}^nQ(a_ix)prodlimits_{j e i}P(a_ix))，那么我们有(ans=(n-2)!prodlimits_{i=1}^na_i[x^{n-2}]F(x))。
我们知道(F(x)=prodlimits_{i=1}^nP(a_ix)sumlimits_{i=1}^nfrac{Q(a_ix)}{P(a_ix)}=exp(sumlimits_{i=1}^n{ln P(a_ix)})sumlimits_{i=1}^nfrac{Q(a_ix)}{P(a_ix)})。
那么我们可以先求出(frac{Q(x)}{P(x)},ln P(x))，将(a_ix)代入并求和等价于(x^k)的系数要乘上(sumlimits_{i=1}^na_i^k)。
现在我们要求的变成了(forall kin[0,n-2],sumlimits_{i=1}^na_i^k)。
考虑答案的生成函数(F(x)=sumlimits_{j=0}^{n-2}sumlimits_{i=1}^na_i^jx^j=sumlimits_{i=1}^nfrac1{1-a_ix})。
我们知道有(ln'(frac1{1-a_ix})=frac{-a_i}{1-a_ix}=-a_isumlimits(a_ix)^j)。
因此我们可以考虑先求出(G(x)=sumlimits_{i=1}^n-a_isumlimits_{j=0}^{n-2}(a_ix)^j)，然后由(F(x)=-xG(x)+n)得出答案。
我们知道(G(x)=sumlimits_{i=1}ln'(frac1{1-a_ix})=(sumlimits_{i=1}ln(frac1{1-a_ix}))'=ln'(prodlimits_{i=1} frac1{1-a_ix}))，那么分治NTT就行了。
总的时间复杂度是(O(nlog n+nlog^2n))，实际上exp的一个(log)和分治的两个(log)差不多。

#include<cstdio>
#include<cstring>
#include<algorithm>
const int N=65537,P=998244353;
int read(){int x;scanf("%d",&x);return x;}
int inc(int a,int b){return a+=b-P,a+(a>>31&P);}
int dec(int a,int b){return a-=b,a+(a>>31&P);}
int mul(int a,int b){return 1ll*a*b%P;}
int pow(int a,int k){int r=1;for(;k;k>>=1,a=mul(a,a))if(k&1)r=mul(a,r);return r;}
int len,fac[N],inv[N],ifac[N],rev[N],w[N];
int getlen(int n){return 1<<(32-__builtin_clz(n));}
void init(int n)
{
    int lim=1<<(len=32-__builtin_clz(n)),g=pow(3,(P-1)/lim);
    w[lim>>1]=1,fac[0]=ifac[0]=inv[0]=fac[1]=ifac[1]=inv[1]=1;
    for(int i=1;i<lim;++i) rev[i]=(rev[i>>1]>>1)|(i&1? lim>>1:0);
    for(int i=(lim>>1)+1;i<lim;++i) w[i]=mul(w[i-1],g);
    for(int i=(lim>>1)-1;i;--i) w[i]=w[i<<1];
    for(int i=2;i<=lim;++i) fac[i]=mul(fac[i-1],i),ifac[i]=mul(ifac[i-1],inv[i]=mul(inv[P%i],P-P/i));
}
void NTT(int*a,int lim,int f)
{
    if(!~f) std::reverse(a+1,a+lim);
    for(int i=0,x=len-__builtin_ctz(lim);i<lim;++i) if(i<rev[i]>>x) std::swap(a[i],a[rev[i]>>x]);
    for(int i=1;i<lim;i<<=1) for(int j=0,d=i<<1;j<lim;j+=d) for(int k=0,x;k<i;++k) x=mul(a[i+j+k],w[i+k]),a[i+j+k]=dec(a[j+k],x),a[j+k]=inc(a[j+k],x);
    if(!~f) for(int i=0,x=P-(P-1)/lim;i<lim;++i) a[i]=mul(a[i],x);
}
void Inv(int*a,int*b,int deg)
{
    if(deg==1) return b[0]=pow(a[0],P-2),void();
    static int t[N];int lim=getlen(deg*2-2);
    Inv(a,b,(deg+1)>>1),memcpy(t,a,deg<<2),memset(t+deg,0,(lim-deg)<<2);
    NTT(t,lim,1),NTT(b,lim,1);
    for(int i=0;i<lim;++i) b[i]=mul(dec(2,mul(b[i],t[i])),b[i]);
    NTT(b,lim,-1),memset(b+deg,0,(lim-deg)<<2);
}
void Der(int*a,int*b,int deg){for(int i=1;i<deg;++i)b[i-1]=mul(a[i],i);b[deg-1]=0;}
void Int(int*a,int*b,int deg){for(int i=1;i<deg;++i)b[i]=mul(a[i-1],inv[i]);b[0]=0;}
void Ln(int*a,int*b,int deg)
{
    static int t[N];int lim=getlen(deg*2-2);
    Inv(a,t,deg),Der(a,b,deg),NTT(t,lim,1),NTT( b,lim,1);
    for(int i=0;i<lim;++i) t[i]=mul(t[i],b[i]);
    NTT(t,lim,-1),Int(t,b,deg),memset(t,0,lim<<2),memset(b+deg,0,(lim-deg)<<2);
}
void Exp(int*a,int*b,int deg)
{
    if(deg==1) return b[0]=1,void();
    static int t[N];int lim=getlen(deg*2-2);
    Exp(a,b,(deg+1)>>1),Ln(b,t,deg);
    for(int i=0;i<deg;++i) t[i]=dec(a[i],t[i]);
    memset(t+deg,0,(lim-deg)<<2),++t[0],NTT(t,lim,1),NTT(b,lim,1);
    for(int i=0;i<lim;++i) b[i]=mul(b[i],t[i]);
    NTT(b,lim,-1),memset(b+deg,0,(lim-deg)<<2),memset(t+deg,0,(lim-deg)<<2);
}
int n,m,ans,a[N],f[16][N],t[5][N];
void solve(int l,int r,int d)
{
    if(l==r) return f[d][0]=1,f[d][1]=P-a[l],void();
    int mid=(l+r)>>1,lim=getlen(r-l+1);
    solve(l,mid,d),solve(mid+1,r,d+1);
    memset(f[d]+mid-l+2,0,(lim-mid+l-2)<<2),memset(f[d+1]+r-mid+1,0,(lim-r+mid-1)<<2),NTT(f[d],lim,1),NTT(f[d+1],lim,1);
    for(int i=0;i<lim;++i) f[d][i]=mul(f[d][i],f[d+1][i]);
    NTT(f[d],lim,-1);
}
int main()
{
    n=read(),m=read(),init(n*2);int lim=getlen(n*2);
    if(n==1) return !printf("%d",m==0);
    for(int i=1;i<=n;++i) a[i]=read();
    solve(1,n,0),Ln(f[0],t[0],n+1),t[0][0]=n;
    for(int i=1;i<=n;++i) t[0][i]=dec(0,mul(i,t[0][i]));
    for(int i=0;i<n;++i) t[1][i]=mul(pow(i+1,m),ifac[i]),t[2][i]=mul(pow(i+1,2*m),ifac[i]);
    Ln(t[1],t[3],n),Inv(t[1],t[4],n),NTT(t[2],lim,1),NTT(t[4],lim,1);
    for(int i=0;i<lim;++i) t[2][i]=mul(t[2][i],t[4][i]);
    NTT(t[2],lim,-1),memset(t[1],0,lim<<2),memset(t[4],0,lim<<2);
    for(int i=0;i<n;++i) t[2][i]=mul(t[2][i],t[0][i]),t[1][i]=mul(t[3][i],t[0][i]);
    Exp(t[1],t[4],n),memset(t[2]+n,0,(lim-n)<<2),NTT(t[2],lim,1),NTT(t[4],lim,1);
    for(int i=0;i<lim;++i) t[2][i]=mul(t[2][i],t[4][i]);
    NTT(t[2],lim,-1),ans=mul(t[2][n-2],fac[n-2]);
    for(int i=1;i<=n;++i) ans=mul(ans,a[i]);
    printf("%d",ans);
}