【bzoj3456】城市规划 dp+多项式求逆

Description

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

Input

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

Output

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

Sample Input

4

Sample Output

38

HINT

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

Sol

设i个顶点的简单连通图的个数为(f[i])，则我们容斥一下，得到：

(f[n]=2^{frac{n(n-1)}{2}}-sum_{i=1}^{n-1}f[i]*C_{n-1}^{i-1}*2^{frac{(n-i)(n-i-1)}{2}})

意义其实就是所有情况减去联通个数不足n的情况，注意后面要枚举哪些点和1号节点连通，所以要乘个组合数。

我们把组合数拆开，然后两边同除以((n-1)!)得到：

(sum_{i=1}^{n}frac{f[i]*2^{frac{(n-i)(n-i-1)}{2}}}{(i-1)!*(n-i)!}=frac{2^{frac{n(n-1)}{2}}}{(i-1)!})

如果我们把式子左边拆成两半，会发现这是一个卷积的形式，所以设这些多项式分别为(A,B,C)，则有(A*B=C)，所以我们通过(C)和(B^{-1})相乘，就得到了(A)，之后取(A)的第n项然后乘个阶乘就是答案。

其实如果(f[i])表示i+1个点的答案的话，好像会简洁很多，实质没有变化（我就是这么写的）。

Code

#include <bits/stdc++.h>
using namespace std;
int n,a[262145],b[262145],c[262145],d[262145],P=1004535809,i,j,k,fac[262145],ifa[262145],inv[262145],wn,w,t,len;
int ksm(int a,long long b){int res=1;for(;b;b>>=1,a=1ll*a*a%P) if(b&1) res=1ll*res*a%P;return res;}
void ntt(int *a,int n,int op)
{
    for(i=k=0;i<n;i++){if(i>k) swap(a[i],a[k]);for(j=(n>>1);(k^=j)<j;j>>=1);}
    for(k=2,wn=ksm(3,op==1?(P-1)/k:P-1-(P-1)/k);k<=n;k<<=1,wn=ksm(3,op==1?(P-1)/k:P-1-(P-1)/k))
        for(i=0,w=1;i<n;i+=k,w=1) for(j=0;j<(k>>1);j++,w=1ll*w*wn%P)
            t=1ll*a[i+j+(k>>1)]*w%P,a[i+j+(k>>1)]=(a[i+j]-t+P)%P,a[i+j]=(a[i+j]+t)%P;
    if(op==-1) for(t=ksm(n,P-2),i=0;i<n;i++) a[i]=1ll*a[i]*t%P;
}
void getinv(int *a,int *b,int *c,int n)
{
    if(n==1){b[0]=ksm(a[0],P-2);return;}
    getinv(a,b,c,n>>1);
    memcpy(c,a,sizeof(int)*n);memset(c+n,0,sizeof(int)*n);
    ntt(c,n<<1,1);ntt(b,n<<1,1);
    for(int i=0;i<(n<<1);i++) b[i]=(b[i]*2ll%P-1ll*b[i]*b[i]%P*c[i]%P+P)%P;
    ntt(b,n<<1,-1);memset(b+n,0,sizeof(int)*n);
}
int main()
{
    scanf("%d",&n);for(len=1;len<=n;len<<=1);
    inv[0]=inv[1]=fac[0]=fac[1]=ifa[0]=ifa[1]=1;
    for(int i=2;i<=len;i++) inv[i]=P-1ll*(P/i)*inv[P%i]%P,fac[i]=1ll*fac[i-1]*i%P,ifa[i]=1ll*ifa[i-1]*inv[i]%P;
    for(int i=0;i<len;i++) d[i]=1ll*ksm(2,1ll*i*(i+1)/2)*ifa[i]%P;
    for(int i=1;i<len;i++) a[i]=1ll*ksm(2,1ll*i*(i-1)/2)*ifa[i]%P;
    a[0]=1;getinv(a,b,c,len);ntt(d,len<<1,1);ntt(b,len<<1,1);
    for(int i=0;i<(len<<1);i++) b[i]=1ll*b[i]*d[i]%P;
    ntt(b,len<<1,-1);printf("%d
",1ll*b[n-1]*fac[n-1]%P);
}