一道图论计数好题：《两弹一星》

题意：定义一个简单无向图的权值为无环的联通块数x的m次方，求由n个点构成的有标号图的权值之和，n<=3e4，m<=20。

直接来说正解，先考虑每个图对答案的贡献：

$$x^m = \sum _{i=0} ^{m} S(m,i) \, i! \; C_{x}^{i}$$

由于m比较小，我们可以直接枚举i，然后对于不同的x，求和的前两项是确定的，问题转化为求最后一项的和。

这里组合数的意义就是，假定我有一个图，它的树联通块数为x，那么对于一个i，它的贡献就是从中任选出i个块的方案数。那么我们转化一下，如果我们确定了树的联通块大小，那么它的贡献就是剩下的点随便连边的方案数，因此我们有了一个做法，设f[i][j]表示i个点的带标号无向图，所有i个点构成了j个树联通块的方案数；g[i]表示i个点的生成树个数，就有如下转移：

$$f[i][j]=\sum _{k=1} ^{i} f[i-k][j-1] * g[k] *C_{i-1} ^{k-1}$$

显然这个式子可以O(n^3)计算，如果我们用FFT优化，只计算j<=m的部分，就可以做到O(mnlogn)的复杂度。

#pragma GCC optimize("-O3")
#include <bits/stdc++.h>
#define ll long long
#define P 998244353
#define G 3
#define mxn 100005
using namespace std;
int bit[32];
inline ll qpow(ll a,ll b=P-2){
    ll ret=1;
    while(b){
        if(b&1) ret=ret*a%P;
        a=a*a%P,b>>=1;
    }
    return ret;
}
int wn[2][18][mxn];
inline void init(int bt){
    ll t1,t2;
    for(int i=1;i<=bt;++i){
        wn[0][i][0]=wn[1][i][0]=1;
        t1=qpow(G,(P-1)/bit[i]),t2=qpow(G,P-1-((P-1)/bit[i]));
        for(int j=1;j<=bit[i-1];++j) wn[0][i][j]=wn[0][i][j-1]*t1%P,wn[1][i][j]=wn[1][i][j-1]*t2%P;
    }
}
inline void ntt(int *s,int n,int o){
    register int i,j,k,l,m,rng;
    for(i=0,j=0;i<n;++i){
        if(i>j) swap(s[i],s[j]);
        for(k=n>>1;(j^=k)<k;k>>=1);
    }
    register int t,u;
    register bool d=(o!=1);
    for(l=2,m=1,rng=1;l<=n;m=l,l<<=1,++rng){
        for(i=0;i<n;i+=l) for(j=0;j<m;++j) {
            t=s[i+j],u=1ll*wn[d][rng][j]*s[i+j+m]%P;
            s[i+j]=(t+u)%P,s[i+j+m]=(t-u+P)%P;
        }
    }
    if(o==-1){
        ll iv=qpow(n);
        for(i=0;i<n;++i) s[i]=1ll*s[i]*iv%P;
    }
}
 
int fac[mxn],inv[mxn];
int S[25][25],f[25][mxn],g[mxn];
int A[mxn],B[mxn];
int n,m;
inline ll C(ll N,ll M){
    if(M>N) return 0;
    return 1ll*fac[N]*inv[M]%P*1ll*inv[N-M]%P;
}
signed main(){
    for(int i=0;i<=30;++i) bit[i]=1<<i;
    scanf("%d%d",&n,&m);
    S[0][0]=1;
    fac[0]=inv[0]=1;
    for(int i=1;i<=n;++i) fac[i]=1ll*fac[i-1]*i%P;
    inv[n]=qpow(fac[n]);
    for(int i=n-1;i>=1;--i) inv[i]=1ll*inv[i+1]*(i+1)%P;
    for(int i=1;i<=20;++i) for(int j=1;j<=i;++j) S[i][j]=(S[i-1][j-1]+1ll*S[i-1][j]*j)%P;
    g[1]=1;
    for(int i=2;i<=n;++i) g[i]=1ll*qpow(i,i-2)*inv[i-1]%P;
    f[0][0]=1;
    int len=1,lim=min(m,n),bt=0;
    while(len<(n<<1)) len<<=1,++bt;
    init(bt);
    for(int i=0;i<=len;++i) B[i]=g[i];
    ntt(B,len,1);
    for(int j=1;j<=lim;++j){
        for(int i=0;i<=len;++i) A[i]=1ll*f[j-1][i]*inv[i]%P;
        ntt(A,len,1);
        for(int i=0;i<=len;++i) A[i]=1ll*A[i]*B[i]%P;
        ntt(A,len,-1);
        for(int i=j;i<=n;++i) f[j][i]=1ll*A[i]*fac[i-1]%P;
    }
    ll ans=0;
    for(int i=0;i<=m;++i){
        ll cr=1ll*S[m][i]*fac[i]%P,cur=0;
        for(int j=i;j<=n;++j) (cur+=1ll*C(n,j)*f[i][j]%P*qpow(2,(1ll*(n-j)*(n-j-1)/2))%P)%=P;
        (ans+=cr*cur%P)%=P;
    }
    printf("%lld",ans);
}

#pragma GCC optimize("-O3")

#include <bits/stdc++.h>

#define ll long long

#define P 998244353

#define G 3

#define mxn 100005

using namespace std;

int bit[32];

inline ll qpow(ll a,ll b=P-2){

ll ret=1;

while(b){

if(b&1) ret=ret*a%P;

a=a*a%P,b>>=1;

}

return ret;

}

int wn[2][18][mxn];

inline void init(int bt){

ll t1,t2;

for(int i=1;i<=bt;++i){

wn[0][i][0]=wn[1][i][0]=1;

t1=qpow(G,(P-1)/bit[i]),t2=qpow(G,P-1-((P-1)/bit[i]));

for(int j=1;j<=bit[i-1];++j) wn[0][i][j]=wn[0][i][j-1]*t1%P,wn[1][i][j]=wn[1][i][j-1]*t2%P;

}

inline void ntt(int *s,int n,int o){

for(i=0,j=0;i<n;++i){

if(i>j) swap(s[i],s[j]);

for(k=n>>1;(j^=k)<k;k>>=1);

}

for(l=2,m=1,rng=1;l<=n;m=l,l<<=1,++rng){

for(i=0;i<n;i+=l) for(j=0;j<m;++j) {

t=s[i+j],u=1ll*wn[d][rng][j]*s[i+j+m]%P;

s[i+j]=(t+u)%P,s[i+j+m]=(t-u+P)%P;

}

if(o==-1){

ll iv=qpow(n);

for(i=0;i<n;++i) s[i]=1ll*s[i]*iv%P;

}

int fac[mxn],inv[mxn];

int S[25][25],f[25][mxn],g[mxn];

int A[mxn],B[mxn];

int n,m;

inline ll C(ll N,ll M){

if(M>N) return 0;

return 1ll*fac[N]*inv[M]%P*1ll*inv[N-M]%P;

}

signed main(){

for(int i=0;i<=30;++i) bit[i]=1<<i;

scanf("%d%d",&n,&m);

S[0][0]=1;

fac[0]=inv[0]=1;

for(int i=1;i<=n;++i) fac[i]=1ll*fac[i-1]*i%P;

inv[n]=qpow(fac[n]);

for(int i=n-1;i>=1;--i) inv[i]=1ll*inv[i+1]*(i+1)%P;

for(int i=1;i<=20;++i) for(int j=1;j<=i;++j) S[i][j]=(S[i-1][j-1]+1ll*S[i-1][j]*j)%P;

g[1]=1;

for(int i=2;i<=n;++i) g[i]=1ll*qpow(i,i-2)*inv[i-1]%P;

f[0][0]=1;

int len=1,lim=min(m,n),bt=0;

while(len<(n<<1)) len<<=1,++bt;

init(bt);

for(int i=0;i<=len;++i) B[i]=g[i];

ntt(B,len,1);

for(int j=1;j<=lim;++j){

for(int i=0;i<=len;++i) A[i]=1ll*f[j-1][i]*inv[i]%P;

ntt(A,len,1);

for(int i=0;i<=len;++i) A[i]=1ll*A[i]*B[i]%P;

ntt(A,len,-1);

for(int i=j;i<=n;++i) f[j][i]=1ll*A[i]*fac[i-1]%P;

}

ll ans=0;

for(int i=0;i<=m;++i){

ll cr=1ll*S[m][i]*fac[i]%P,cur=0;

for(int j=i;j<=n;++j) (cur+=1ll*C(n,j)*f[i][j]%P*qpow(2,(1ll*(n-j)*(n-j-1)/2))%P)%=P;

(ans+=cr*cur%P)%=P;

}

printf("%lld",ans);

}

题面（虽然没卵用）：

本文隐藏内容登陆后才可以浏览

Narcissus' Blog

一道图论计数好题：《两弹一星》

24.4.18 考试小记

JSOI2018 Day2 部分题解

Scott

叨叨几句... NOTHING

取消回复

Narcissus' Blog