bzoj 3456 城市规划-CFANZ编程社区

3456: 城市规划

Time Limit: 40 Sec Memory Limit: 256 MB
Submit: 732 Solved: 413
[Submit][Status][Discuss]
Description

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

Input

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

Output

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

Sample Input

Sample Output

HINT

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

Source

【分析】

一道经典(恶心)的图的计数问题…

题意：求n个点无向连通图个数

令f(i)为i个点无向连通图的个数,。

则f(i)=2C(i,2)−∑i−1j=1f(j)∗C(i−1,j−1)∗2C(i−j,2)

证明：考虑用i个点所有无向图的个数减去不连通的无向图的个数。枚举1号点所在的连通块的大小j，方案为f(j)，那么我们要从2~i的点中取j−1个点和1号点形成一个连通块，方案为C(i−1,j−1)，剩下的i−j个点之间任意连边。方案数位2C(i−j,2)

两边同时除以fac(i−1)，并展开C(i−1,j−1)，然后移项。

得到 2C(i,2)fac(i−1)=∑ij=1f(j)fac(j−1)∗2C(i−j,2)fac(i−j)

发现这他喵的是个卷积…我已经写不动公式了…太TM麻烦了…

然后构造三个多项式

A[i]=f(i)fac(i−1)

B[i]=2C(i,2)fac(i)

C[i]=2C(i,2)fac(i−1)

A∗B=C

A=C∗B−1

得到A中xn项的系数，再乘以fac(n−1)即为答案。
用到多项式求逆，不会的同学自行百度

【代码】

//bzoj 3456 城市规划 
#include<bits/stdc++.h>
#define M 130000
#define ll long long
#define fo(i,j,k) for(int i=j;i<=k;i++)
using namespace std;
const int mxn=600005;
const int mod=1004535809;
int n,m,L,N,ans;
int fac[mxn],inv[mxn];
int a[mxn],b[mxn],c[mxn],R[mxn];
inline void init()
{
    fac[0]=inv[0]=inv[1]=1;
    fo(i,1,M) fac[i]=(ll)fac[i-1]*i%mod;
    fo(i,2,M) inv[i]=(ll)(mod-mod/i)*inv[mod%i]%mod;
    fo(i,1,M) inv[i]=(ll)inv[i]*inv[i-1]%mod;
}
inline int power(int x,int k)
{
    int res=1;
    while(k)
    {
        if(k&1) res=(ll)res*x%mod;
        x=(ll)x*x%mod,k>>=1;
    }
    return res;
}
inline void NTT(int *a,int f)
{
    fo(i,0,n-1) if(i<R[i]) swap(a[i],a[R[i]]);
    for(int i=1;i<n;i<<=1)
    {
        int wn=power(3,(mod-1)/(i<<1));
        for(int j=0;j<n;j+=(i<<1))
        {
            int w=1;
            for(int k=0;k<i;k++,w=(ll)w*wn%mod)
            {
                int x=a[j+k],y=(ll)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 rev=power(n,mod-2);
        fo(i,0,n-1) a[i]=(ll)a[i]*rev%mod;
    }
}
inline void getinv(int Nt)
{
    if(Nt==1) {b[0]=1;return;}
    getinv(Nt+1>>1);
    n=Nt;int i,m=n+n;
    for(n=1,L=0;n<=m;n<<=1) L++;
//  printf("%d %d\n",Nt,n);
    fo(i,0,n-1) R[i]=(R[i>>1]>>1)|((i&1)<<L-1);
    fo(i,0,Nt-1) c[i]=a[i];fo(i,Nt,n) c[i]=0;
    NTT(b,1),NTT(c,1);
    fo(i,0,n)
    {
        b[i]=(ll)(2-(ll)b[i]*c[i]%mod)*b[i]%mod;
        if(b[i]<0) b[i]+=mod;
    }
    NTT(b,-1);
    fo(i,Nt,n) b[i]=0;
}
int main()
{
    init();
    scanf("%d",&N);
    fo(i,0,N)
      a[i]=(ll)power(2,(ll)i*(i-1)/2%(mod-1))*inv[i]%mod;
    m=N;for(n=1;n<=m;n<<=1) L++;
    getinv(n);
    memset(c,0,sizeof c);
    memset(a,0,sizeof a);
    m=N+N;for(n=1,L=0;n<=m;n<<=1) L++;
    fo(i,0,n-1) R[i]=(R[i>>1]>>1)|((i&1)<<L-1);
    fo(i,1,N)
      c[i]=(ll)power(2,(ll)i*(i-1)/2%(mod-1))*inv[i-1]%mod;
    NTT(b,1),NTT(c,1);
    fo(i,0,n) a[i]=(ll)b[i]*c[i]%mod;
    NTT(a,-1);

    printf("%lld\n",(ll)a[N]*fac[N-1]%mod);
    return 0;
}