BZOJ 3451Normal (点分治+FFT)

题意

BZOJ3451

题解

Orz yyb题解

CODE

#pragma GCC optimize ("O2")
#include <bits/stdc++.h>
using namespace std;
const double Pi = acos(-1.0);
const int MAXN = 30005;

int fir[MAXN], to[MAXN<<1], nxt[MAXN<<1], cnt;
inline void link(int u, int v) {
	to[++cnt] = v; nxt[cnt] = fir[u]; fir[u] = cnt;
	to[++cnt] = u; nxt[cnt] = fir[v]; fir[v] = cnt;
}
bool ban[MAXN];
int Get_Size(int u, int ff) {
	int re = 1;
	for(int i = fir[u]; i; i = nxt[i])
		if(to[i] != ff && !ban[to[i]])
			re += Get_Size(to[i], u);
	return re;
}
inline int Get_G(int u, int ff, int Size, int &g) {
	int re = 1; bool flg = 1;
	for(int i = fir[u]; i; i = nxt[i])
		if(to[i] != ff && !ban[to[i]]) {
			int tmp = Get_G(to[i], u, Size, g); re += tmp;
			if(tmp<<1 > Size) flg = 0;
		}
	if((Size-re)<<1 > Size) flg = 0;
	if(flg) g = u; return re;
}
const int MAX = 70000;
struct cp {
	double x, y; cp(){} cp(double x, double y):x(x), y(y){}
	inline cp operator +(const cp &o)const { return cp(x + o.x, y + o.y); }
	inline cp operator -(const cp &o)const { return cp(x - o.x, y - o.y); }
	inline cp operator *(const cp &o)const { return cp(x*o.x-y*o.y, x*o.y+y*o.x); }
}A[MAX], B[MAX], W[MAX];
int rv[MAX];
inline void FFT(cp *a, int N, int flg) {
	for(int i = 0; i < N; ++i) if(rv[i] < i) swap(a[rv[i]], a[i]);
	cp u, v, w;
	for(int i = 2; i <= N; i<<=1)
		for(int j = 0; j < N; j += i)
			for(int k = j; k < j+i/2; ++k) {
				w = W[N/i*(k-j)]; w.y *= flg;
				u = a[k], v = a[k+i/2] * w;
				a[k] = u + v;
				a[k+i/2] = u - v;
			}
	if(flg == -1) for(int i = 0; i < N; ++i) a[i].x /= N;
}
inline void Mul(int *a, int *b, int n, int m, int *c) {
	int N, lg = 0; for(N = 1; N <= (n+m); ++lg, N<<=1);
	for(int i = 0; i < N; ++i) rv[i] = (rv[i>>1]>>1) | ((i&1)<<(lg-1));
	for(int i = 2; i <= N; i<<=1)
		for(int k = 0; k < i/2; ++k)
			W[N/i*k] = cp(cos(2*Pi/i*k), sin(2*Pi/i*k));
	for(int i = 0; i <= n; ++i) A[i].x = a[i];
	for(int i = 0; i <= m; ++i) B[i].x = b[i];
	FFT(A, N, 1), FFT(B, N, 1);
	for(int i = 0; i < N; ++i) A[i] = A[i] * B[i];
	FFT(A, N, -1);
	for(int i = 0; i <= n+m; ++i) c[i] = int(round(A[i].x));
	for(int i = 0; i < N; ++i) A[i] = B[i] = cp(0, 0);
}
int a[MAXN], b[MAXN], c[MAXN], la, lb, Ans[MAXN];
void dfs(int u, int ff, int d) {
	++b[d]; lb = max(lb, d);
	for(int i = fir[u]; i; i = nxt[i])
		if(to[i] != ff && !ban[to[i]])
			dfs(to[i], u, d+1);
}
void TDC(int u) {
	int Size = Get_Size(u, 0);
	Get_G(u, 0, Size, u);
	la = 0; a[0] = 1;
	for(int i = fir[u]; i; i = nxt[i])
		if(!ban[to[i]]) {
			lb = 0; dfs(to[i], u, 1);
			Mul(a, b, la, lb, c);
			for(int i = 0; i <= la+lb; ++i) Ans[i] += c[i], c[i] = 0;
			for(int i = 0; i <= lb; ++i) a[i] += b[i], b[i] = 0;
			la = max(la, lb);
		}
	for(int i = 0; i <= la; ++i) a[i] = 0;
	ban[u] = 1;
	for(int i = fir[u]; i; i = nxt[i])
		if(!ban[to[i]]) TDC(to[i]);
}
int n;
int main () {
    scanf("%d", &n);
	for(int i = 1, x, y; i < n; ++i)
		scanf("%d%d", &x, &y), link(x+1, y+1);
	TDC(1);
	double ans = n;
	for(int i = 1; i <= n; ++i) ans += 2.0 * Ans[i] / (i+1);
	printf("%.4f
", ans);
}

看似十分正常而且420ms过了美滋滋。但是实际上最坏情况是$O(n^2logn)$。可以发现我们点分治算答案的时候每一次的复杂度都与前面子树的最大深度也就是代码中的la。这样下来如果第一个子树深度为$O(n)$，每次都是$O(nlogn)$的，所以就$O(n^2logn)$了。

正确写法是把所有子树的多项式和的平方后减去平方的和，就是两两乘积的和了。

还有一种就是点分治时把儿子按深度排个序，然后每次每棵子树的时间复杂度只跟自己深度有关（因为前面的深度都小于等于自己的深度），然后深度小于等于子树大小。如果是子树大小的时间复杂度则容易证明是$O(nlogn)$的。总时间复杂度就是$O(nlog^2n)$