bzoj

[bzoj 3456] 城市规划

求n个点的有标号简单 (无重边或自环) 无向连通图数目, 答案模 1004535809 (479*2^21 + 1). 以下讨论的均是有标号简单无向连通图.

首先得列出递推式. 一个想法是, 枚举点1所在连通块大小, 计算n个点无向非连通图的数目, 用所有的减去它.

设f(n)为n个点无向连通图数目, g(n)为n个点无向非连通图数目, h(n)为n个点所有无向图数目. 那么, 

f(n) + g(n) = h(n) = 2^(n(n-1)/2)
g(n) = sum { C(n-1, i-1) f(i) h(n-i) | 1 <= i < n }

隐约看到卷积的形式. 把组合数用阶乘展开, 并用f替换g: 

f(n) = h(n) - (n-1)! sum { f(i)/(i-1)! h(n-i)/(n-i)! | 1 <= i < n }

但这是递推式, 无法直接计算.

据说 分治FFT = FFT + CDQ分治, 好像可以......

欲知f[l,r], 先计算[l,m]对[l,m]的贡献 (递归), 再计算[l,m]对[m+1,r]的贡献 (FFT), 最后计算[m+1,r]对[m+1,r]的贡献 (递归).

FFT具体把哪两个数列卷起来呢?

记a[i] = f(i)/(i-1)!, b[i] = h(i)/i!, w = (r-l+1)/2. 

a[l] b[w] + a[l+1] b[w-1] + ... + a[l+w-1] b[1] -> f[l+w]
a[l] b[w+1] + a[l+1] b[w] + ... + a[l+w-1] b[2] -> f[l+w+1]
.
.
.
a[l] b[2w-1] + a[l+1] b[2w-2] + ... + a[l+w-1] b[w-1] -> f[l+2w-1]

于是, 应该把 

a[l], a[l+1], ..., a[l+w-1], 0, 0, ..., 0
b[1], b[2], ..., b[w], b[w+1], ..., b[2w-1]
他俩卷起来.

由于本题要求对 479*2^21 + 1 取模, FFT的所有运算在模意义下进行, 即NTT.

typedef long long ll;

const int MOD = 1004535809, L = 17, N = (1<<L) + 1;

inline ll fpm(ll x, ll n)
{
	ll y = 1;
	for (; n; n >>= 1, (x *= x) %= MOD)
		if (n & 1)
			(y *= x) %= MOD;
	return y;
}

inline ll inverse(ll x)
{
	return fpm(x, MOD-2);
}

namespace Trans {
	const ll g[2] = {3, 334845270};
	// w_n = g^((p-1)/n)

	int rev[N];
	ll w[2][L+1], half[L+1];

	void init(int n, int s)
	{
		rep (i, 0, s + 1)
			half[i] = inverse(1<<i);
		rep (i, 0, 2) {
			w[i][s] = fpm(g[i], (MOD-1)/n);
			per (j, s, 1)
				w[i][j-1] = w[i][j] * w[i][j] % MOD;
		}
	}

	void ntt(ll x[], int n, int s, int d=0)
	{
		// x c : 0 x -> x^T 0 -> c x^T
		rep (i, 1, n) {
			rev[i] = (rev[i>>1]>>1) | ((i&1)<<(s-1));
			if (rev[i] < i)
				swap(x[i], x[rev[i]]);
		}
		
		for (int m = 1, k = 1; m < n; m <<= 1, ++k) {
			ll w0 = w[d][k];
			for (int i = 0; i < n; i += m<<1) {
				ll _w = 1;
				rep (j, 0, m) {
					ll t = _w * x[i+m+j] % MOD;
					x[i+m+j] = (x[i+j] - t + MOD) % MOD;
					(x[i+j] += t) %= MOD;
					(_w *= w0) %= MOD;
				}
			}
		}

		if (d) rep (i, 0, n) (x[i] *= half[s]) %= MOD;
	}

	void convol(ll x[], ll y[], ll z[], int n, int s)
	{
		ntt(x, n, s);
		ntt(y, n, s);
		rep (i, 0, n) z[i] = x[i] * y[i] % MOD;
		ntt(z, n, s, 1);
	}
}

using Trans::convol;

ll fac[2][N], f[N], h[N];

void init(int n)
{
	fac[0][0] = 1;
	rep (i, 1, n+1) {
		fac[0][i] = fac[0][i-1] * i % MOD;
		h[i] = fpm(2, (ll)i*(i-1)/2);
	}
	fac[1][n] = inverse(fac[0][n]);
	per (i, n, 1) fac[1][i-1] = fac[1][i] * i % MOD;
}

void solve(int l, int r, int s)
{
	if (l == r) {
		f[l] = l == 1 ? 1 : (h[l] - fac[0][l-1] * f[l] % MOD + MOD) % MOD;
		return;
	}
	
	int m = (l+r)/2, w = r-m;
	solve(l, m, s-1);
	
	static ll x[N], y[N], z[N];

	rep (i, 0, w)
		x[i] = f[l+i] * fac[1][l+i-1] % MOD;
	rep (i, w, w<<1)
		x[i] = 0;
	rep (i, 0, w<<1)
		y[i] = h[i+1] * fac[1][i+1] % MOD;
	
	convol(x, y, z, w<<1, s);

	rep (i, 0, w)
		(f[m+1+i] += z[w-1+i]) %= MOD;

	solve(m+1, r, s-1);
}

int main()
{
	int n;
	scanf("%d", &n);
	int s = 0;
	while (1<<s < n) ++s;
	init(1<<s);
	Trans::init(1<<s, s);

	solve(1, 1<<s, s);

	printf("%lld\n", f[n]);

	return 0;
}