BZOJ 3451 Tyvj1953 Normal（点分治+FFT）

题目

题解：
这个概率统计神了。
考虑在以$i$为根的点分区域内（即我下一步将要选出来分治的点是$i$时的树），$j$被计算到贡献的概率。
那么这个概率等于$(i,j)$的路径上$i$是第一个被选出来当分治中心的点的概率。
在这之前路径上的点都没有被选，最后一次选出$i$的那次，路径上的点都在同一个点分区域内。
那么路径上的点被选到的概率都一样，所以我们要求的概率就是$\frac{1}{dis(i,j)+1}$
那么只需要求出$dis(i,j)=k$的方案数，这用一个树上点分治+FFT就可以简单计算了。
注意点分治时处理子树的复杂度只能和当前子树大小有关（就是在指指点点那些每次把当前子树和前面前缀的多项式相乘的写法），不然可以卡到$O(n^2\log n)$。
写边分治就可以不用管了。
很奇怪的是，学校OJ上FFT300ms，NTT150ms，BZOJ上FFT600ms，NTT1200ms。

AC Code：

#include<bits/stdc++.h>
#define maxn 70005
#define pb push_back
#define Pi 3.1415926535897932384626433832795
using namespace std;

int n;
vector<int>G[maxn];

struct cplx{
	double r,i;
	cplx(double r=0,double i=0):r(r),i(i){}
	cplx operator +(const cplx &B)const{ return cplx(r+B.r,i+B.i); }
	cplx operator -(const cplx &B)const{ return cplx(r-B.r,i-B.i); }
	cplx operator *(const cplx &B)const{ return cplx(r*B.r-i*B.i,r*B.i+i*B.r); }
	cplx conj()const{ return cplx(r,-i); }
}w[maxn]={1};
int lg[maxn],r[maxn],wlen;
void Init(int n){
	for(wlen=1;n>=2*wlen;wlen<<=1);
	for(int i=1;i<=wlen;i++) w[i]=cplx(cos(Pi*i/wlen),sin(Pi*i/wlen));
	for(int i=2;i<=2*wlen;i++) lg[i] = lg[i>>1] + 1;
}

void FFT(cplx *A,int n,int tp){
	int lgn = lg[n];
	for(int i=1;i<n;i++) r[i] = (r[i>>1]>>1)|((i&1)<<(lgn-1));
	for(int i=1;i<n;i++) if(i<r[i]) swap(A[i],A[r[i]]);
	for(int L=2;L<=n;L<<=1)
		for(int l=L>>1,st=0,inc=wlen/l;st<n;st+=L)
			for(int k=st,x=0;k<st+l;k++,x+=inc){
				cplx t = (tp == 1 ? w[x] : w[x].conj()) * A[k+l];
				A[k+l] = (A[k] - t) , A[k] = (A[k] + t);
			}
	if(tp==-1)
		for(int i=0;i<n;i++)
			A[i].r /= n;
}

int Min=0x3f3f3f3f,rt=-1;
int siz[maxn],vis[maxn];

void dfs(int u,int tsz,int ff){
	int Max = 0;siz[u] = 1;
	for(int i=0,v;i<G[u].size();i++)
		if(!vis[v=G[u][i]] && v!=ff)
			dfs(v,tsz,u),
			Max = max(Max , siz[v]),
			siz[u] += siz[v];
	Max = max(Max , tsz - siz[u]);
	if(Max < Min)
		rt = u , Min = Max;
}

int gert(int u,int tsz){
	Min=0x3f3f3f3f;
	dfs(u,tsz,0);
	return rt;
}

int mxd;
cplx ar[3][maxn],ret[maxn];

void ser(int u,int ff,int dis){
	mxd = max(mxd , dis);
	ar[0][dis].r++ , ar[2][dis].r++ , siz[u] = 1;
	for(int i=0,v;i<G[u].size();i++)
		if(!vis[v=G[u][i]] && v!=ff){
			ser(v,u,dis+1);
			siz[u] += siz[v];
		}
}

int ct = 0;
void Solve(int u){
	vis[u] = 1;
	ar[0][0] = 1;
	int mcd = 0;
	for(int i=0,v;i<G[u].size();i++)
		if(!vis[v=G[u][i]]){
			mxd = 0;
			ser(v,u,1);
			int len = 1<<(lg[2*mxd]+1);
			FFT(ar[2],len,1);
			for(int i=0;i<len;i++)
				ar[2][i] = ar[2][i] * ar[2][i];
			FFT(ar[2],len,-1);
			for(int i=0;i<len;i++)
				ar[1][i] = (ar[1][i] + ar[2][i]),
				ar[2][i] = 0;
			mcd = max(mcd , mxd);
		}
	int len = 1<<(lg[2*mcd]+1);
	FFT(ar[0],len,1);
	for(int i=0;i<len;i++)
		ar[0][i] = ar[0][i] * ar[0][i];
	FFT(ar[0],len,-1);
	for(int i=0;i<len;i++)
		ret[i] = ret[i] + ar[0][i] - ar[1][i],
		ar[0][i] = ar[1][i] = 0;
		
	for(int i=0,v;i<G[u].size();i++)
		if(!vis[v=G[u][i]])
			Solve(gert(v,siz[v]));
}

int main(){
	scanf("%d",&n);
	Init(2*n);
	for(int i=1,u,v;i<n;i++){
		scanf("%d%d",&u,&v);
		G[u].pb(v),G[v].pb(u);
	}
	Solve(gert(1,n));
	double ans = 0;
	for(int i=0;i<=n;i++)
		if(ret[i].r)
			ans += ret[i].r / (i+1.0);
	printf("%.4lf",ans);
}