BZOJ3451 Normal 期望、点分治、NTT

BZOJCH传送门

题目大意：给出一棵树，求对其进行随机点分治的复杂度期望

---

可以知道一个点的贡献就是其点分树上的深度，也就是这个点在点分树上的祖先数量+1。

根据期望的线性性，考虑一个点对((x,y))在何时(x)能够是(y)的祖先，那么在(x)到(y)的路径上的所有点中(x)必须要是第一个被选中的点，否则要么点(y)变成点(x)的祖先，要么点(x)和点(y)会被分在不同的子树中。那么我们要求的就是(sum frac{1}{dis_{x,y}})，其中(dis_{x,y})表示(x)到(y)路径上的点的数量。直接使用(NTT)统计路径长度即可。

注意在某一个分治中心进行合并时必须按照深度从小到大合并，否则复杂度是不正确的。

#include<bits/stdc++.h>
#define INF 0x7fffffff
using namespace std;

inline int read(){
	int a = 0;
	char c = getchar();
	while(!isdigit(c))
		c = getchar();
	while(isdigit(c)){
		a = a * 10 + c - 48;
		c = getchar();
	}
	return a;
}

const int MAXN = 1e5 + 10 , MOD = 998244353 , G = 3 , INV = 332748118;
struct Edge{
	int end , upEd;
}Ed[MAXN << 1];
int head[MAXN] , ind[MAXN] , size[MAXN] , ans[MAXN];
int N , need , nSize , mSize , mInd , cntEd;
bool vis[MAXN];
long long inv_need;
vector < vector < int > > v;

inline void addEd(int a , int b){
	Ed[++cntEd].end = b;
	Ed[cntEd].upEd = head[a];
	head[a] = cntEd;
}

inline int poww(long long a , int b){
	int times = 1;
	while(b){
		if(b & 1)
			times = times * a % MOD;
		a = a * a % MOD;
		b >>= 1;
	}
	return times;
}

void NTT(vector < int > &v , int type){
	for(int i = 0 ; i < need ; ++i)
		if(i < ind[i])
			swap(v[i] , v[ind[i]]);
	for(int i = 1 ; i < need ; i <<= 1){
		int wn = poww(type == 1 ? G : INV , (MOD - 1) / (i << 1));
		for(int j = 0 ; j < need ; j += i << 1){
			long long w = 1;
			for(int k = 0 ; k < i ; ++k , w = w * wn % MOD){
				int x = v[j + k] , y = v[i + j + k] * w % MOD;
				v[j + k] = x + y >= MOD ? x + y - MOD : x + y;
				v[i + j + k] = x < y ? x - y + MOD : x - y;
			}
		}
	}
}

void getSize(int x){
	vis[x] = 1;
	++nSize;
	for(int i = head[x] ; i ; i = Ed[i].upEd)
		if(!vis[Ed[i].end])
			getSize(Ed[i].end);
	vis[x] = 0;
}

void getRoot(int x){
	vis[x] = 1;
	int maxN = 0;
	size[x] = 1;
	for(int i = head[x] ; i ; i = Ed[i].upEd)
		if(!vis[Ed[i].end]){
			getRoot(Ed[i].end);
			maxN = max(maxN , size[Ed[i].end]);
			size[x] += size[Ed[i].end];
		}
	maxN = max(maxN , nSize - size[x]);
	if(maxN < mSize){
		mSize = maxN;
		mInd = x;
	}
	vis[x] = 0;
}

void calc(vector < int > &v , int x , int l){
	while(v.size() <= l)
		v.push_back(0);
	++v[l];
	vis[x] = 1;
	for(int i = head[x] ; i ; i = Ed[i].upEd)
		if(!vis[Ed[i].end])
			calc(v , Ed[i].end , l + 1);
	vis[x] = 0;
}

bool cmp(vector < int > a , vector < int > b){
	return a.size() < b.size();
}

void solve(int x){
	nSize = 0;
	mSize = INF;
	getSize(x);
	getRoot(x);
	x = mInd;
	vis[x] = 1;
	v.clear();
	vector < int > cur , ano;
	for(int i = head[x] ; i ; i = Ed[i].upEd)
		if(!vis[Ed[i].end]){
			cur.clear();
			calc(cur , Ed[i].end , 1);
			v.push_back(cur);
		}
	cur.clear();
	sort(v.begin() , v.end());
	need = 1;
	for(int i = 0 ; i < v.size() ; ++i){
		while(need < cur.size() + v[i].size())
			need <<= 1;
		inv_need = poww(need , MOD - 2);
		for(int j = 1 ; j < need ; ++j)
			ind[j] = (ind[j >> 1] >> 1) | (j & 1 ? need >> 1 : 0);
		while(cur.size() < need)
			cur.push_back(0);
		while(v[i].size() < need)
			v[i].push_back(0);
		ano.clear();
		NTT(cur , 1);
		NTT(v[i] , 1);
		for(int j = 0 ; j < need ; ++j){
			ano.push_back(1ll * cur[j] * v[i][j] % MOD);
			cur[j] = cur[j] + v[i][j] >= MOD ? cur[j] + v[i][j] - MOD : cur[j] + v[i][j];
		}
		NTT(ano , -1);
		NTT(cur , -1);
		for(int j = 0 ; j < need ; ++j){
			ans[j] += ano[j] * inv_need % MOD;
			cur[j] = cur[j] * inv_need % MOD;
		}
		while(cur.back() == 0)
			cur.pop_back();
	}
	for(int i = 1 ; i < need ; ++i)
		ans[i] += cur[i];
	for(int i = head[x] ; i ; i = Ed[i].upEd)
		if(!vis[Ed[i].end])
			solve(Ed[i].end);
}

int main(){
	N = read();
	for(int i = 1 ; i < N ; ++i){
		int a = read() , b = read();
		addEd(a , b);
		addEd(b , a);
	}
	solve(1);
	long double sum = 0;
	for(int i = 1 ; i <= N ; ++i)
		sum += ans[i] * 2.0 / (i + 1);
	cout << fixed << setprecision(4) << sum + N;
	return 0;
}