BZOJ3456 城市规划 (CDQ+NTT)

BZOJ-3456 (CDQ+NTT)

题意：求\(n\)个有标号点联通图的方案数

考虑减去\(n\)个点不连通的方案数

对于当前的\(i\)个点枚举1号点所在连通块大小为\(j(1<j<i)\)，则方案数为\(C(i-1,j-1)\cdot dp_j\cdot 2^{(i-j)(i-j-1)/2}\)

即选出剩下的\(j-1\)个点，其他点随意

得到\(dp_i=2^{i(i-1)/2}-C(i-1,j-1)\cdot dp_j\cdot 2^{(i-j)(i-j-1)/2}\)

可以看到是一个与作差有关的转移，可以用\(CDQ+NTT\)优化

#include<bits/stdc++.h>
using namespace std;
 
#define reg register
#define pb push_back
typedef long long ll;
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)
 
template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); }
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); }
 
char IO;
int rd(){
    int s=0,f=0;
    while(!isdigit(IO=getchar())) if(IO=='-') f=1;
    do s=(s<<1)+(s<<3)+(IO^'0');
    while(isdigit(IO=getchar()));
    return f?-s:s;
}
 
const int N=2e5+10,P=1004535809;
 
int n;
ll dp[N],Inv[N],Fac[N];
ll qpow(ll x,ll k){
    ll res=1;
    for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
    return res;
}
 
int rev[N];
ll A[N],B[N],w[N],g[N];
void NTT(int n,ll *a,int f) {
    rep(i,0,n-1) if(i>rev[i]) swap(a[i],a[rev[i]]);
    for(reg int i=1;i<n;i<<=1) {
        int len=n/i/2;
        for(reg int l=0;l<n;l+=2*i) {
            for(reg int j=l;j<l+i;++j) {
                ll t=a[j+i]*w[(j-l)*len]%P;
                a[j+i]=a[j]-t; (a[j+i]<0&&(a[j+i]+=P));
                a[j]=a[j]+t; (a[j]>=P&&(a[j]-=P));
            }
        }
    }
    if(f==-1) {
        ll base=qpow(n,P-2);
        rep(i,0,n-1) a[i]=a[i]*base%P;
    }
}
 
 
void Solve(int l,int r) {
    if(l==r) return;
    if(r-l+1<=20) {
        rep(i,l,r) rep(j,i+1,r) dp[j]=(dp[j]-dp[i]*Fac[j-1]%P*Inv[i-1]%P*Inv[j-i]%P*g[j-i]%P+P)%P;
        return;
    } // 底层展开优化
    int mid=(l+r)>>1;
    Solve(l,mid);
    int R=1,cc=-1;
    while(R<=r-l+1) R<<=1,cc++;
    rep(i,1,R) rev[i]=(rev[i>>1]>>1)|((i&1)<<cc);
    rep(i,0,R) A[i]=B[i]=0;
    rep(i,l,mid) A[i-l]=dp[i]*Inv[i-1]%P;
    rep(i,1,r-l) B[i]=Inv[i]*g[i]%P;
    w[0]=1,w[1]=qpow(3,(P-1)/R); rep(i,2,R-1) w[i]=w[i-1]*w[1]%P;
    NTT(R,A,1),NTT(R,B,1);
    rep(i,0,R-1) A[i]=A[i]*B[i]%P;
    w[0]=1,w[1]=qpow((P+1)/3,(P-1)/R); rep(i,2,R-1) w[i]=w[i-1]*w[1]%P;
    NTT(R,A,-1);
    rep(i,mid+1,r) dp[i]=(dp[i]-Fac[i-1]*A[i-l]%P+P)%P;
    Solve(mid+1,r);
}
 
int main(){
    n=rd();
    rep(i,1,n) dp[i]=g[i]=qpow(2,1ll*i*(i-1)/2); //g[i]存下i个点随便取的方案数
    Inv[0]=Inv[1]=Fac[0]=Fac[1]=1;
    rep(i,2,n) {
        Inv[i]=(P-P/i)*Inv[P%i]%P;
        Fac[i]=Fac[i-1]*i%P;
    }
    rep(i,1,n) Inv[i]=Inv[i-1]*Inv[i]%P;
    Solve(1,n);
    printf("%lld\n",(dp[n]%P+P)%P);
}