[BZOJ3456] 城市规划

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

Solution

设(f(i))为(i)个点的题目所求的答案，(g(i))为(i)个点随便连的方案数，那么我们枚举(1)号点所在联通块的大小，可得：

[g(n)=sum_{i=1}^nf(i)inom{n-1}{i-1}g(n-i)=2^{inom{n}{2}}=2^{n(n-1)/2} ]

容易知道这样可以不重不漏的算到所有状态。

我们把(g)带进去：

[2^{n(n-1)/2}=sum_{i=1}^nf(i)inom{n-1}{i-1}2^{inom{n-i}{2}} ]

两边除以((n-1)!)：

[frac{2^{n(n-1)/2}}{(n-1)!}=sum_{i=1}^nfrac{f(i)}{(i-1)!} imes frac{2^{inom{n-i}{2}}}{(n-i)!} ]

注意到这是个卷积的形式，用生成函数表示一下：

[A(x)=sum_{n=0}^inftyfrac{2^{n(n-1)/2}}{(n-1)!}x^n\ B(x)=sum_{n=0}^inftyfrac{2^{n(n-1)/2}}{n!}x^n\ F(x)=sum_{n=0}^inftyfrac{f(i)}{(n-1)!}x^n ]

那么：

[A(x)=F(x)B(x),F(x)=A(x)B^{-1}(x) ]

那么直接套多项式求逆的板子就好了。

#include<bits/stdc++.h>
using namespace std;
 
void read(int &x) {
    x=0;int f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}
 
void print(int x) {
    if(x<0) putchar('-'),x=-x;
    if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('
');}

#define lf double
#define ll long long 

const int maxn = 6e5+10;
const int inf = 1e9;
const lf eps = 1e-8;
const int mod = 1004535809;

int n,w[maxn],rw[maxn],pos[maxn],f[maxn],a[maxn],b[maxn],c[maxn],mxn,N,bit,fac[maxn],ifac[maxn],g[maxn];

int qpow(int aa,int x) {
	int res=1;
	for(;x;x>>=1,aa=1ll*aa*aa%mod) if(x&1) res=1ll*res*aa%mod;
	return res;
}

void prepare() {
	w[0]=1,w[1]=qpow(3,(mod-1)/mxn);
	for(int i=2;i<=mxn;i++) w[i]=1ll*w[i-1]*w[1]%mod;
	rw[0]=1,rw[1]=qpow(qpow(3,mod-2),(mod-1)/mxn);
	for(int i=2;i<=mxn;i++) rw[i]=1ll*rw[i-1]*rw[1]%mod;
}

void ntt(int *r,int op) {
	for(int i=1;i<N;i++) if(pos[i]>i) swap(r[i],r[pos[i]]);
	for(int i=1,d=mxn>>1;i<N;i<<=1,d>>=1) 
		for(int j=0;j<N;j+=i<<1)
			for(int k=0;k<i;k++) {
				int x=r[j+k],y=1ll*r[i+j+k]*(op==1?w:rw)[k*d]%mod;
				r[j+k]=(x+y)%mod,r[i+j+k]=(x-y+mod)%mod;
			}
	if(op==-1) {
		int inv=qpow(N,mod-2);
		for(int i=0;i<N;i++) r[i]=1ll*r[i]*inv%mod;
	}
}

int tmp1[maxn],tmp2[maxn];

void get_inv(int *r,int *s,int len) {
	if(len==1) return s[0]=qpow(r[0],mod-2),void();
	get_inv(r,s,len>>1);
	for(int i=0;i<len>>1;i++) tmp1[i]=s[i];
	for(int i=len>>1;i<len;i++) tmp1[i]=0;
	for(int i=0;i<len;i++) tmp2[i]=r[i];

	for(bit=0,N=1;N<=len;N<<=1,bit++);
	for(int i=1;i<N;i++) pos[i]=pos[i>>1]>>1|((i&1)<<(bit-1));
	for(int i=len;i<N;i++) tmp1[i]=tmp2[i]=0;
	ntt(tmp1,1),ntt(tmp2,1);
	for(int i=0;i<N;i++) s[i]=((2ll*tmp1[i]%mod-1ll*tmp2[i]*tmp1[i]%mod*tmp1[i]%mod)%mod+mod)%mod;
	ntt(s,-1);
	for(int i=len;i<N;i++) s[i]=0;
}

int main() {
	read(n);
	for(mxn=1;mxn<=n<<1;mxn<<=1);
	prepare();
	fac[0]=ifac[0]=c[0]=1;
	for(int i=1;i<=mxn;i++) fac[i]=1ll*fac[i-1]*i%mod;
	ifac[mxn]=qpow(fac[mxn],mod-2);
	for(int i=mxn-1;i;i--) ifac[i]=1ll*ifac[i+1]*(i+1)%mod;
	for(int i=1;i<=n;i++) g[i]=qpow(2,(1ll*i*(i-1)/2)%(mod-1));
	for(int i=1;i<=n;i++) c[i]=1ll*g[i]*ifac[i]%mod,b[i]=1ll*g[i]*ifac[i-1]%mod;
	
	get_inv(c,a,mxn>>1);
	for(bit=0,N=1;N<mxn;N<<=1,bit++);
	for(int i=1;i<N;i++) pos[i]=pos[i>>1]>>1|((i&1)<<(bit-1));
	ntt(a,1),ntt(b,1);
	for(int i=0;i<N;i++) a[i]=1ll*a[i]*b[i]%mod;
	ntt(a,-1);
	printf("%lld
",(1ll*a[n]*fac[n-1]%mod+mod)%mod);
	return 0;
}