BZOJ3456：城市规划 （多项式求逆/分治+NTT）

权限题传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=3456

---

题目大意：给定n，求n个点的简单无向连通图的数量（每个点标号不同，图无重边，无自环）。$n<=1.3*10^5$，模1004535809。

---

题目分析：这题一看就很套路。首先n个点的简单无向图有$2^{C_{n}^{2}}$种，记为g[n]。然后根据容斥原理，我们枚举1号点所在的连通块大小k$(0<k<n)$，以减去不合法的情况。记答案为f[n]，则1号点的连通块有f[k]种情况，其它点的连边有g[n-k]种情况，1号点所在连通块的点编号有$C_{n-1}^{k-1}$种选择，于是：

$$f[n]=g[n]-\sum_{k=1}^{n-1}f[k]*g[n-k]*C_{n-1}^{k-1}$$

经过简单的推导，我们可以得到：

$$f[n]=g[n]-(n-1)!*\sum_{k=1}^{n-1}\frac{f[k]}{(k-1)!}*\frac{g[n-k]}{(n-k)!}$$

由于计算f[n]需要知道f[1~n-1]，而计算贡献又是个卷积的形式，这已经可以用分治+NTT解决了，时间复杂度$O(n\log^2(n))$。然而我们还可以继续推导倒。

同样地，通过考虑1号节点所在连通块大小，我们可以得到g[n]和f数组的关系：

$$g[n]=(n-1)!*\sum_{k=1}^{n}\frac{f[k]}{(k-1)!}*\frac{g[n-k]}{(n-k)!}$$

两边同时除以$(n-1)!$，可得：

$$\frac{g[n]}{(n-1)!}=\sum_{k=1}^{n}\frac{f[k]}{(k-1)!}*\frac{g[n-k]}{(n-k)!}$$

记$A[n]=\frac{g[n]}{(n-1)!},B[k]=\frac{f[k]}{(k-1)!},C[n-k]=\frac{g[n-k]}{(n-k)!}$。只要我们强行令B[0]=0（但C[0]不为0），就可以让式子中的k从0开始枚举，这样就完全是一个卷积了。由于$A=B*C$，所以$B=A*C^{-1}$，多项式求逆即可（注意A[0]根据原式算出来应该也为0）。时间复杂度$O(n\log(n))$。

---

CODE：

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
using namespace std;

const int maxn=1000000;
const long long M=1004535809;
const long long g=3;
typedef long long LL;

LL A[maxn];
LL B[maxn];

int Rev[maxn];
int N,Lg;

LL F[maxn];
LL G[maxn];

int fac[maxn];
int nfac[maxn];
int n;

LL Pow(LL x,LL y)
{
    if (!y) return 1LL;
    LL temp=Pow(x,y>>1);
    temp=temp*temp%M;
    if (y&1) temp=temp*x%M;
    return temp;
}

void DFT(LL *a,int f)
{
    for (int i=0; i<N; i++)
        if (i<Rev[i]) swap(a[i],a[ Rev[i] ]);

    for (int len=2; len<=N; len<<=1)
    {
        int mid=(len>>1);
        LL e=Pow(g,(M-1)/len);
        if (f==-1) e=Pow(e,M-2);

        for (LL *p=a; p!=a+N; p+=len)
        {
            LL wn=1;
            for (int i=0; i<mid; i++)
            {
                LL temp=wn*p[mid+i]%M;
                p[mid+i]=(p[i]-temp+M)%M;
                p[i]=(p[i]+temp)%M;
                wn=wn*e%M;
            }
        }
    }
}

void NTT(bool rev)
{
    for (int i=0; i<N; i++)
    {
        Rev[i]=0;
        for (int j=0; j<Lg; j++)
            if (i&(1<<j)) Rev[i]|=(1<<(Lg-j-1));
    }

    DFT(A,1);
    DFT(B,1);
    if (rev) for (int i=0; i<N; i++) A[i]=A[i]*A[i]%M;
    for (int i=0; i<N; i++) A[i]=A[i]*B[i]%M;
    DFT(A,-1);

    LL inv=Pow(N,M-2);
    for (int i=0; i<N; i++) A[i]=A[i]*inv%M;
}

void Poly_Rev(int m)
{
    if (m==1) F[0]=Pow(G[0],M-2);
    else
    {
        Poly_Rev(m>>1);

        N=2*m;
        Lg=0;
        while ((1<<Lg)!=N) Lg++;

        for (int i=0; i<m; i++) A[i]=F[i],B[i]=G[i];
        for (int i=m; i<N; i++) A[i]=B[i]=0;
        NTT(true);

        for (int i=0; i<m; i++) F[i]=2LL*F[i]%M,F[i]=(F[i]-A[i]+M)%M;
    }
}

int main()
{
    freopen("3456.in","r",stdin);
    freopen("3456.out","w",stdout);

    scanf("%d",&n);

    N=1;
    while (N<=n+2) N<<=1;
    fac[0]=1;
    for (LL i=1; i<N; i++) fac[i]=fac[i-1]*i%M;
    for (int i=0; i<N; i++) nfac[i]=Pow(fac[i],M-2);

    for (LL i=0; i<N; i++) G[i]=Pow(2, i*(i-1LL)/2LL )*nfac[i]%M;
    Poly_Rev(N);

    N=1,Lg=0;
    while (N<2*n+3) N<<=1,Lg++;
    for (int i=0; i<=n; i++)
    {
        A[i]=F[i];
        if (i) B[i]=Pow(2, i*(i-1LL)/2LL )*nfac[i-1]%M;
        else B[i]=0;
    }
    for (int i=n+1; i<N; i++) A[i]=B[i]=0;
    NTT(false);
    for (int i=1; i<=n; i++)
    {
        A[i]=A[i]*fac[i-1]%M;
        if (i==n) printf("%I64d\n",A[i]);
    }

    return 0;
}