【BZOJ 3456】城市规划

http://www.lydsy.com/JudgeOnline/problem.php?id=3456
(f(n))表示n个点有标号无向连通图的数目。
dp:(f(n)=2^{nchoose 2}-sumlimits_{i=1}^{n-1}f(i){n-1choose i-1}2^{n-ichoose 2})
这是一个可以用分治FFT(O(nlog^2n))做的式子。
移项,分配阶乘使之变为卷积的形式:$$sum_{i=0}nfrac{f(i)}{(i-1)!} imesfrac{2{n-ichoose 2}}{(n-i)!}=frac{2^{nchoose 2}}{(n-1)!}$$
(当(i=0)时默认(frac{f(0)}{(0-1)!}=0)
然后可以多项式求逆一波。
设多项式(A(x))在模(x^n)意义下的逆多项式为(B_n(x)),可以在任意一篇博客上找到推导过程,这里直接写结论:

[B_n(x)equiv B_{leftlceilfrac n2 ight ceil}(x)left(2-A(x)B_{leftlceilfrac n2 ight ceil}(x) ight)pmod {x^n} ]

时间复杂度(O(nlog n))

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;

const int p = 1004535809;
const int N = 131072 << 1;

int ipow(int a, int b) {
	int r = 1, w = a;
	while (b) {
		if (b & 1) r = 1ll * r * w % p;
		w = 1ll * w * w % p;
		b >>= 1;
	}
	return r;
}

int n, rev[N];
ll G[33], nG[33], f[N], ni[N], nifrac[N], t[N];

void DFT(ll *a, int n, int flag) {
	for (int i = 0; i < n; ++i) if (i < rev[i]) swap(a[i], a[rev[i]]);
	int now = -1;
	for (int m = 2; m <= n; m <<= 1) {
		int mid = m >> 1; ++now;
		ll wn = flag == 1 ? G[now] : nG[now];
		for (int i = 0; i < n; i += m) {
			ll w = 1;
			for (int j = 0; j < mid; ++j) {
				ll u = a[i + j], v = a[i + j + mid] * w % p;
				a[i + j] = (u + v) % p;
				a[i + j + mid] = (u - v + p) % p;
				w = w * wn % p;
			}
		}
	}
	
	if (flag == -1) {
		ll nii = ipow(n, p - 2);
		for (int i = 0; i < n; ++i)
			(a[i] *= nii) %= p;
	}
}

void INV(ll *A, ll *B, int n) {
	if (n == 1) {B[0] = ipow(A[0], p - 2); return;}
	
	INV(A, B, (n + 1) >> 1);
	int len = 1, bl = -1, nn = (n << 1) - 1;
	for (; len < nn; len <<= 1, ++bl);
	for (int i = 1; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << bl);
	for (int i = 0; i < n; ++i) t[i] = A[i];
	for (int i = n; i < len; ++i) t[i] = 0;
	
	DFT(t, len, 1); DFT(B, len, 1);
	for (int i = 0; i < len; ++i) B[i] = B[i] * ((2 - t[i] * B[i] % p + p) % p) % p;
	DFT(B, len, -1);
	for (int i = n; i < len; ++i) B[i] = 0;
}

ll A[N], B[N], C[N], nA[N];

int main() {
	scanf("%d", &n);
	if (n <= 2) {puts("1"); return 0;}
	
	int len = 1, bl = -1, nn = ((n + 1) << 1) - 1;
	for (; len < nn; len <<= 1, ++bl);
	G[bl] = ipow(3, (p - 1) / len); nG[bl] = ipow(G[bl], p - 2);
	for (int i = bl - 1; i >= 0; --i) G[i] = G[i + 1] * G[i + 1] % p, nG[i] = nG[i + 1] * nG[i + 1] % p;
	
	ni[1] = 1; nifrac[0] = nifrac[1] = 1;
	for (int i = 2; i <= n; ++i) {
		ni[i] = (p - p / i) * ni[p % i] % p;
		nifrac[i] = nifrac[i - 1] * ni[i] % p;
	}
	
	A[0] = 1;
	ll last = 1, C = 1;
	for (int i = 1; i <= n; ++i) {
		A[i] = last * nifrac[i] % p;
		B[i] = last * nifrac[i - 1] % p;
		last = last * ((C <<= 1) %= p) % p;
	}
	
	INV(A, nA, n + 1);
	
	for (int i = 1; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << bl);
	DFT(nA, len, 1); DFT(B, len, 1);
	for (int i = 0; i < len; ++i) (B[i] *= nA[i]) %= p;
	DFT(B, len, -1);
	
	ll noww = 1;
	for (int i = 2; i <= n; ++i) {
		(noww *= (i - 1)) %= p;
		(B[i] *= noww) %= p;
	}
	
	printf("%lld
", B[n]);
	return 0;
}
原文地址:https://www.cnblogs.com/abclzr/p/6809759.html