[BZOJ 3456] 城市规划

3456. 城市规划

Description

 刚刚解决完电力网络的问题, 阿狸又被领导的任务给难住了.
 刚才说过, 阿狸的国家有n个城市, 现在国家需要在某些城市对之间建立一些贸易路线, 使得整个国家的任意两个城市都直接或间接的连通. 为了省钱, 每两个城市之间最多只能有一条直接的贸易路径. 对于两个建立路线的方案, 如果存在一个城市对, 在两个方案中是否建立路线不一样, 那么这两个方案就是不同的, 否则就是相同的. 现在你需要求出一共有多少不同的方案.
 好了, 这就是困扰阿狸的问题. 换句话说, 你需要求出n个点的简单(无重边无自环)无向连通图数目.
 由于这个数字可能非常大, 你只需要输出方案数mod 1004535809(479 * 2 ^ 21 + 1)即可.

Input

 仅一行一个整数n(<=130000)
 

Output

 仅一行一个整数, 为方案数 mod 1004535809.
 

Sample Input

3

Sample Output

4

Hint

 对于 100%的数据, n <= 130000

题解

首先我们可以发现, 含有 $n$ 个点的有标号简单无向图的个数是 $2^{{n choose 2}}$.

然后我们有标号简单无向图和有标号简单连通无向图的关系: 每个有标号简单无向图都可以划分为若干联通分量, 而这些联通分量就是有标号简单联通无向图.

如何体现这个划分的过程? 首先我们使用EGF来保持标号运算, 若 $n$ 个点的有标号简单联通无向图的个数的生成函数为 $F(x)$, 则 $frac {F(x)^k} {k!}$ 即为含有 $k$ 个联通分量的简单无向图个数的生成函数.

为啥最后有一个分母 $k!$ 呢? 我们把 $k$ 个 $F(x)$ 卷起来, 每个划分方案的每个排列都会被统计一次. 因为直接卷积时每个分量是有序合并的. 而我们在用联通无向图构造无向图的时候显然不能考虑联通分量的顺序, 所以我们除一个 $k!$.

设 $n$ 个点的简单无向图的个数的EGF为 $G(x)$, 我们就有了:

$$ G(x) = sum_{k=0}^infty frac {F(x)^k}{k!} $$

又根据函数 $exp(x)=e^x$ 的泰勒展开式:

$$ e^x=sum_{k=0}^infty frac{x^k}{k!} $$

我们将 $F(x)$ 代入 $x$ 中, 得到:

$$ G(x) = sum_{k=0}^infty frac {F(x)^k}{k!} = exp(F(x)) $$

于是我们就有了:

$$ F(x)=ln(G(x)) $$

所以我们构造出 $G(x)=sumlimits_{k=0}^infty 2^{{kchoose 2}} frac {x^k}{k!} $ 之后求它的 $ln$ 就可以得到 $F(x)$ 了.

代码实现

直接上全家桶了...然而并不全都用得到...所以码长比较怪异

#include <bits/stdc++.h>

const int G=3;
const int DFT=1;
const int IDFT=-1;
const int MAXN=550000;
const int MOD=1004535809;
const int INV2=(MOD+1)>>1;
const int PHI=MOD-1;

typedef long long int64;
typedef std::vector<int> Poly;

Poly Sqrt(Poly);
void Read(Poly&);
Poly Inverse(Poly);
Poly Ln(const Poly&);
Poly Exp(const Poly&);
void Print(const Poly&);
void NTT(Poly&,int,int);
Poly Pow(const Poly&,int);
Poly Integral(const Poly&);
Poly Derivative(const Poly&);
Poly operator*(Poly,Poly);
Poly operator/(Poly,Poly);
Poly operator%(Poly,Poly);
Poly operator+(const Poly&,const Poly&);
Poly operator-(const Poly&,const Poly&);

int n;
int rev[MAXN];
int fact[MAXN];

int NTTPre(int);
int Pow(int,int,int);
int Pow(int,int64,int);

int main(){
    scanf("%d",&n);
    fact[0]=1;
    for(int i=1;i<=n;i++)
        fact[i]=1ll*i*fact[i-1]%MOD;
    Poly a(n+1);
    a[0]=1;
    for(int i=1;i<=n;i++)
        a[i]=1ll*Pow(2,1ll*i*(i-1)/2,MOD)*Pow(fact[i],MOD-2,MOD)%MOD;
    a=Ln(a);
    printf("%d
",int(1ll*a[n]*fact[n]%MOD));
    return 0;
}

void Read(Poly& a){
    for(auto& i:a)
        scanf("%d",&i);
}

void Print(const Poly& a){
    for(auto i:a)
        printf("%d ",i);
    puts("");
}

Poly Pow(const Poly& a,int k){
    Poly log=Ln(a);
    for(auto& i:log)
        i=1ll*i*k%MOD;
    return Exp(log);
}

Poly Sqrt(Poly a){
    int len=a.size();
    if(len==1)
        return Poly(1,int(sqrt(a[0])));
    else{
        Poly b=a;
        b.resize((len+1)>>1);
        b=Sqrt(b);
        b.resize(len);
        Poly inv=Inverse(b);
        int bln=NTTPre(inv.size()+a.size());
        NTT(a,bln,DFT);
        NTT(inv,bln,DFT);
        for(int i=0;i<bln;i++)
            a[i]=1ll*a[i]*INV2%MOD*inv[i]%MOD;
        NTT(a,bln,IDFT);
        for(int i=0;i<len;i++)
            b[i]=(1ll*b[i]*INV2%MOD+a[i])%MOD;
        return b;
    }
}

Poly Exp(const Poly& a){
    size_t len=1;
    Poly ans(1,1),one(1,1);
    while(len<(a.size()<<1)){
        len<<=1;
        Poly b=a;
        b.resize(len);
        ans=ans*(one-Ln(ans)+b);
        ans.resize(len);
    }
    ans.resize(a.size());
    return ans;
}

Poly Ln(const Poly& a){
    Poly ans=Integral(Derivative(a)*Inverse(a));
    ans.resize(a.size());
    return ans;
}

Poly Integral(const Poly& a){
    int len=a.size();
    Poly ans(len+1);
    for(int i=1;i<len;i++)
        ans[i]=1ll*a[i-1]*Pow(i,MOD-2,MOD)%MOD;
    return ans;
}

Poly Derivative(const Poly& a){
    int len=a.size();
    Poly ans(len-1);
    for(int i=1;i<len;i++)
        ans[i-1]=1ll*a[i]*i%MOD;
    return ans;
}

Poly operator/(Poly a,Poly b){
    int n=a.size()-1,m=b.size()-1;
    Poly ans(1);
    if(n>=m){
        std::reverse(a.begin(),a.end());
        std::reverse(b.begin(),b.end());
        b.resize(n-m+1);
        ans=Inverse(b)*a;
        ans.resize(n-m+1);
        std::reverse(ans.begin(),ans.end());
    }
    return ans;
}

Poly operator%(Poly a,Poly b){
    int n=a.size()-1,m=b.size()-1;
    Poly ans;
    if(n<m)
        ans=a;
    else
        ans=a-(a/b)*b;
    ans.resize(m);
    return ans;
}

Poly operator*(Poly a,Poly b){
    int len=a.size()+b.size()-1;
    int bln=NTTPre(len);
    NTT(a,bln,DFT);
    NTT(b,bln,DFT);
    for(int i=0;i<bln;i++)
        a[i]=1ll*a[i]*b[i]%MOD;
    NTT(a,bln,IDFT);
    a.resize(len);
    return a;
}

Poly operator+(const Poly& a,const Poly& b){
    Poly ans(std::max(a.size(),b.size()));
    std::copy(a.begin(),a.end(),ans.begin());
    for(size_t i=0;i<b.size();i++)
        ans[i]=(ans[i]+b[i])%MOD;
    return ans;
}

Poly operator-(const Poly& a,const Poly& b){
    Poly ans(std::max(a.size(),b.size()));
    std::copy(a.begin(),a.end(),ans.begin());
    for(size_t i=0;i<b.size();i++)
        ans[i]=(ans[i]+MOD-b[i])%MOD;
    return ans;
}

Poly Inverse(Poly a){
    int len=a.size();
    if(len==1)
        return Poly(1,Pow(a[0],MOD-2,MOD));
    else{
        Poly b(a);
        b.resize((len+1)>>1);
        b=Inverse(b);
        int bln=NTTPre(b.size()*2+a.size());
        NTT(a,bln,DFT);
        NTT(b,bln,DFT);
        for(int i=0;i<bln;i++)
            b[i]=(2ll*b[i]%MOD-1ll*b[i]*b[i]%MOD*a[i]%MOD+MOD)%MOD;
        NTT(b,bln,IDFT);
        b.resize(len);
        return b;
    }
}

void NTT(Poly& a,int len,int opt){
    a.resize(len);
    for(int i=0;i<len;i++)
        if(rev[i]>i)
            std::swap(a[i],a[rev[i]]);
    for(int i=1;i<len;i<<=1){
        int step=i<<1;
        int wn=Pow(G,(PHI+opt*PHI/step)%PHI,MOD);
        for(int j=0;j<len;j+=step){
            int w=1;
            for(int k=0;k<i;k++,w=1ll*w*wn%MOD){
                int x=a[j+k];
                int y=1ll*a[j+k+i]*w%MOD;
                a[j+k]=(x+y)%MOD;
                a[j+k+i]=(x-y+MOD)%MOD;
            }
        }
    }
    if(opt==IDFT){
        int inv=Pow(len,MOD-2,MOD);
        for(int i=0;i<len;i++)
            a[i]=1ll*a[i]*inv%MOD;
    }
}

int NTTPre(int n){
    int bln=1,bct=0;
    while(bln<n){
        bln<<=1;
        ++bct;
    }
    for(int i=0;i<bln;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bct-1));
    return bln;
}

inline int Pow(int a,int n,int p){
    int ans=1;
    while(n>0){
        if(n&1)
            ans=1ll*a*ans%p;
        a=1ll*a*a%p;
        n>>=1;
    }
    return ans;
}

inline int Pow(int a,int64 n,int p){
    int ans=1;
    while(n>0){
        if(n&1)
            ans=1ll*a*ans%p;
        a=1ll*a*a%p;
        n>>=1;
    }
    return ans;
}