【多项式求逆】[BZOJ3456]城市规划

题目描述

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

Input

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

Output

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

Sample Input

3

Sample Output

4

HINT

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

分析

令f(i)表示有i个点的无向连通图的个数，g(i)表示有i个点的无向图的个数。
显然

g(i)=2(i2)

我们通过枚举一号节点所在连通块的大小来转移。

g(i)=∑j=1i(i−1j−1)f(j)g(i−j)=∑j=1i(i−1)!(j−1)!(n−j)!f(j)g(i−j)g(i)(i−1)!=f(j)(j−1)!g(i−j)(i−j)!

令Ai=g(i)(i−1)!  Bi=f(j)(j−1)!  Ci=g(i−j)(i−j)!
特别地，根据式子可以看出，A0=0

A(x)=B(x)C(x)B(x)=A(x)÷C(x)

将其放在modxn+1意义下

B(x)=A(x)C(x)−1

然后就可以求出Bn，然后算出f(n)

代码

#include<cstdio>
#include<algorithm>
#define MOD 1004535809
#define g 3
#define MAXN 130000
using namespace std;
int N,F[MAXN*2+10000],G[MAXN*2+10000],C[MAXN*2+10000],n,fac[MAXN+10],inv[MAXN+10],GR[MAXN*2+10000],tmp[MAXN*2+10000];
void Read(int &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
        }
}
int quick_pow(int a,int b){
    int ret(1);
    while(b){
        if(b&1)
            ret=1ll*ret*a%MOD;
        a=1ll*a*a%MOD;
        b>>=1;
    }
    return ret;
}
void ntt(int *a,int N,int f){
    int i,j,k,t;
    for(i=1,j=0;i<N-1;i++){
        for(t=N;j^=t>>=1,~j&t;);
        if(i<j)
            swap(a[i],a[j]);
    }
    for(i=1;i<N;i<<=1){
        t=i<<1;
        int wn(quick_pow(g,(MOD-1)/t));
        for(j=0;j<N;j+=t){
            int w(1);
            for(k=0;k<i;k++,w=1ll*w*wn%MOD){
                int x(a[j+k]),y(1ll*w*a[j+k+i]%MOD);
                a[j+k]=(x+y)%MOD,a[j+k+i]=(x-y+MOD)%MOD;
            }
        }
    }
    if(f==-1){
        reverse(a+1,a+N);
        int inv=quick_pow(N,MOD-2);
        for(i=0;i<N;i++)
            a[i]=1ll*a[i]*inv%MOD;
    }
}
void Divide_Conqure(int n,int *A,int *B){
    if(n==1){
        B[0]=quick_pow(A[0],MOD-2);
        return;
    }
    int mid((n+1)>>1),i,N;
    Divide_Conqure(mid,A,B);
    for(N=1;N<(n<<1)-1;N<<=1);
    for(i=0;i<n;i++)
        tmp[i]=A[i];
    for(i=n;i<N;i++)
        tmp[i]=0;
    ntt(tmp,N,1);
    ntt(B,N,1);
    for(i=0;i<N;i++)
        B[i]=((B[i]*2-1ll*tmp[i]*B[i]%MOD*B[i]%MOD)+MOD)%MOD;
    ntt(B,N,-1);
    for(i=n;i<N;i++)
        B[i]=0;
}
void prepare(){
    fac[0]=inv[0]=1;
    for(int i=1;i<=n;i++){
        fac[i]=1ll*fac[i-1]*i%MOD;
        inv[i]=quick_pow(fac[i],MOD-2);
    }
}
void solve(){
    int i;
    for(i=1;i<=n;i++){
        G[i]=1ll*quick_pow(2,1ll*i*(i-1)/2%(MOD-1))*inv[i]%MOD;
        C[i]=1ll*quick_pow(2,1ll*i*(i-1)/2%(MOD-1))*inv[i-1]%MOD;
    }
    for(N=1;N<n*2+1;N<<=1);
    G[0]=1;
    ntt(C,N,1);
    Divide_Conqure(n+1,G,GR);
    ntt(GR,N,1);
    for(i=0;i<N;i++)
        F[i]=1ll*GR[i]*C[i]%MOD;
    ntt(F,N,-1);
}
int main()
{
    Read(n);
    prepare();
    solve();
    printf("%lld\n",1ll*F[n]*fac[n-1]%MOD);
}