【BZOJ3456】城市规划（分治NTT）

题面

求出 $n$ 个点的简单(无重边无自环)无向连通图数目取模 $1004535809$ （一个质数，$479\ast 2^{21}+1$）的结果。

$n\le 130000$.

题解

直接求不好着手，用暴力差得远，那我们就绞尽脑汁想一想递推式。

令 $f(x)$ 表示 $x$ 个点的简单无向连通图数目，怎么从更小的状态转移过来？

我们不妨从编号最大的点考虑，想办法把整张图分裂成两个子图方便转移。

如果直接把最大的点的邻接边都删掉，有可能会分离出多个连通块，如果只删一条邻接边，有可能还是联通的。

因此有这样一种策略，如果删掉最大的点，一定会剩下一个或多个连通块。我们要选择其中的一个，把它完全分离出去。不妨就选择第二大的点所在的连通块。那么就断掉最大点与该连通块之间的所有边，使整张图变成两个连通块（第二大点所在连通块，以及其他点包括最大点所在的连通块），得到 $f(x)$ 的前驱状态。于是，可以枚举第二大的点所在的连通块大小，然后转移：
$$f(x)=\sum _{i=1}^{x-1}f(i)\cdot f(x-i)\cdot \left(\genfrac{}{}{0ex}{}{x-2}{x-i-1}\right)\cdot (2^{x-i}-1)$$

求和符号后的四项分别是第二大的点所在连通块方案数、剩下点所在连通块方案数、节点编号合并方案数、连边方案数。

接着就可以变成下面的式子：
$$\begin{gathered} f(x)=\sum _{i=1}^{x-1}f(i)\cdot f(x-i)\cdot \frac{(x-2)!}{(x-i-1)!(i-1)!}\cdot (2^{x-i}-1) \\ =(x-2)!\sum _{i=1}^{x-1}\left(\frac{f(i)}{(i-1)!}\right)\cdot \left(\frac{f(x-i)}{(x-i-1)!}\cdot (2^{x-i}-1)\right) \end{gathered}$$

这就是一个自卷积形式，可以用分治NTT解决。

时间复杂度 $O(n{\log }^{{}_2}n)$ 。

CODE

$1004535809$ 的原根为 3 。

#include<map>
#include<set>
#include<cmath>
#include<stack>
#include<vector>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
#define MAXN 130005
#define LL long long
#define ENDL putchar('\n')
#define DB double
#define lowbit(x) (-(x) & (x))
#define FI first
#define SE second
int xchar() {
	static const int maxn = 1000000;
	static char b[maxn];
	static int pos = 0,len = 0;
	if(pos == len) pos = 0,len = fread(b,1,maxn,stdin);
	if(pos == len) return -1;
	return b[pos ++];
}
//#define getchar() xchar()
LL read() {
	LL f = 1,x = 0;int s = getchar();
	while(s < '0' || s > '9') {if(s<0)return -1;if(s=='-')f=-f;s = getchar();}
	while(s >= '0' && s <= '9') {x = (x<<1) + (x<<3) + (s^48);s = getchar();}
	return f*x;
}
void putpos(LL x) {if(!x)return ;putpos(x/10);putchar((x%10)^48);}
void putnum(LL x) {
	if(!x) {putchar('0');return ;}
	if(x<0) putchar('-'),x = -x;
	return putpos(x);
}
void AIput(LL x,int c) {putnum(x);putchar(c);}

const int MOD = 1004535809;
const int pmr = 3;
int n,m,s,o,k;
int fac[MAXN],inv[MAXN],invf[MAXN],pw2[MAXN];
int qkpow(int a,int b) {
	int res = 1;
	while(b > 0) {
		if(b & 1) res = res *1ll* a % MOD;
		a = a *1ll* a % MOD; b >>= 1; 
	}return res;
}
int rev[MAXN<<2],xm[MAXN<<2],om;
void NTT(int *s,int n,int op) {
	for(int i = 1;i < n;i ++) {
		rev[i] = (rev[i>>1]>>1) | ((i&1) ? (n>>1):0);
		if(rev[i] < i) swap(s[rev[i]],s[i]);
	}
	om = qkpow(pmr,(MOD-1)/n); xm[0] = 1;
	if(op < 0) om = qkpow(om,MOD-2);
	for(int i = 1;i <= n;i ++) xm[i] = xm[i-1]*1ll*om%MOD;
	for(int k = 2,t = n>>1;k <= n;k <<= 1,t >>= 1) {
		for(int j = 0;j < n;j += k) {
			for(int i = j,l = 0;i < j+(k>>1);i ++,l += t) {
				int A = s[i],B = s[i+(k>>1)];
				s[i] = (0ll+A + xm[l]*1ll*B) % MOD;
				s[i+(k>>1)] = ((0ll+A-xm[l]*1ll*B) % MOD+MOD)%MOD;
			}
		}
	}
	if(op < 0) {
		int invn = qkpow(n,MOD-2);
		for(int i = 0;i < n;i ++) s[i] = s[i] *1ll* invn % MOD;
	}return ;
}
int A[MAXN<<2],B[MAXN<<2],f[MAXN<<2];
void solve(int l,int r) {
	if(l == r) {
		if(l == 1) f[l] = 1;
		return ;
	}
	int mid = (l + r) >> 1;
	solve(l,mid);
	int le1 = mid-l+1,le2 = min(r-l,mid);
	int le = 1; while(le < le1+le2) le <<= 1;
	for(int i = 0;i < le;i ++) A[i] = B[i] = 0;
	for(int i = l;i <= mid;i ++) A[i-l] = f[i]*1ll*invf[i-1]%MOD;
	for(int j = 1;j <= le2;j ++) B[j] = f[j]*1ll*invf[j-1]%MOD*(pw2[j]+MOD-1)%MOD;
	NTT(A,le,1); NTT(B,le,1);
	for(int i = 0;i < le;i ++) A[i] = A[i] *1ll* B[i] % MOD;
	NTT(A,le,-1);
	for(int i = 0;i < le;i ++) {
		if(i+l > mid && i+l <= r) (f[i+l] += fac[i+l-2]*1ll*A[i]%MOD) %= MOD;
	}
	le1 = mid-l+1,le2 = min(r-l,l-1);
	le = 1; while(le < le1+le2) le <<= 1;
	for(int i = 0;i < le;i ++) A[i] = B[i] = 0;
	for(int i = l;i <= mid;i ++) A[i-l] = f[i]*1ll*invf[i-1]%MOD*(pw2[i]+MOD-1)%MOD;
	for(int j = 1;j <= le2;j ++) B[j] = f[j]*1ll*invf[j-1]%MOD;
	NTT(A,le,1); NTT(B,le,1);
	for(int i = 0;i < le;i ++) A[i] = A[i] *1ll* B[i] % MOD;
	NTT(A,le,-1);
	for(int i = 0;i < le;i ++) {
		if(i+l > mid && i+l <= r) (f[i+l] += fac[i+l-2]*1ll*A[i]%MOD) %= MOD;
	}
	solve(mid+1,r);
	return ;
}
int main() {
	n = read();
	fac[0]=fac[1]=inv[0]=inv[1]=invf[0]=invf[1]=pw2[0]=1;
	pw2[1] = 2;
	for(int i = 2;i <= n;i ++) {
		fac[i] = fac[i-1] *1ll* i % MOD;
		inv[i] = (MOD-inv[MOD%i]) *1ll* (MOD/i) % MOD;
		invf[i] = invf[i-1] *1ll* inv[i] % MOD;
		pw2[i] = (pw2[i-1]<<1) % MOD;
	}
	solve(1,n);
	AIput(f[n],'\n');
	return 0;
}