【BZOJ3451】【Tyvj1953】—Normal(点分治+NTT)

传送门

Description

某天WJMZBMR学习了一个神奇的算法:树的点分治!
这个算法的核心是这样的:
消耗时间=0
Solve(树 a)
消耗时间 += a 的 大小
如果 a 中 只有 1 个点
退出
否则在a中选一个点x,在a中删除点x
那么a变成了几个小一点的树,对每个小树递归调用Solve
我们注意到的这个算法的时间复杂度跟选择的点x是密切相关的。
如果x是树的重心,那么时间复杂度就是O(nlogn)
但是由于WJMZBMR比较傻逼,他决定随机在a中选择一个点作为x!
Sevenkplus告诉他这样做的最坏复杂度是O(n^2)
但是WJMZBMR就是不信>_<。。。
于是Sevenkplus花了几分钟写了一个程序证明了这一点。。。你也试试看吧_
现在给你一颗树,你能告诉WJMZBMR他的傻逼算法需要的期望消耗时间吗?(消耗时间按在Solve里面的那个为标准)

Input

第一行一个整数n,表示树的大小
接下来n-1行每行两个数a,b,表示a和b之间有一条边
注意点是从0开始标号的

Output
一行一个浮点数表示答案
四舍五入到小数点后4位
如果害怕精度跪建议用long double或者extended

Sample Input

3

0 1

1 2

Sample Output

5.6667

HINT

n<=30000

考虑计算一个点 i i i j j j的贡献
那显然只有 j j j为分治中心的时候 i i i对它才有贡献
考虑如果分治中心不在路径 i i i~ j j j内那显然不会产生任何影响
而如果一次中心把 i i i j j j隔开了那就没有贡献了
那也就是说 i i i j j j的时候是分治中心是 j j j,而不是 i i i~ j j j路径上的其他点
也就是说概率是 1 d i s ( i , j ) \frac 1 {dis(i,j)} dis(i,j)1
那答案就是 ∑ i = 1 n ∑ j = 1 n 1 d i s ( i , j ) \sum_{i=1}^{n}\sum_{j=1}^{n}\frac{1}{dis(i,j)} i=1nj=1ndis(i,j)1

而也就是可以求出每个距离出现了多少次,除以距离就可以了
而这个是可以通过点分治+ N T T NTT NTT O ( n l o g 2 n ) O(nlog^2n) O(nlog2n)的时间求出来

复杂度 O ( n l o g 2 n ) O(nlog^2n) O(nlog2n)

#include<bits/stdc++.h>
using namespace std;
const int RLEN=1<<20|1;
#define ll long long
inline char gc(){
	static char ibuf[RLEN],*ib,*ob;
	(ib==ob)&&(ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
	return (ib==ob)?EOF:*ib++;
}
#define gc getchar
inline int read(){
	char ch=gc();
	int res=0,f=1;
	while(!isdigit(ch))f^=ch=='-',ch=gc();
	while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
	return f?res:-res;
}
const int mod=998244353,g=3;
const int M=120005,N=30005;
inline int add(int a,int b){
	return a+b>=mod?a+b-mod:a+b;
}
inline int dec(int a,int b){
	return a>=b?a-b:a-b+mod;
}
inline int mul(int a,int b){
	return 1ll*a*b>=mod?1ll*a*b%mod:a*b;
}
inline int ksm(int a,int b,int res=1){
	for(;b;b>>=1,a=mul(a,a))(b&1)?(res=mul(res,a)):0;return res;
}
int rev[M],A[M],lim,tim;
inline void ntt(int *f,int kd){
	for(int i=0;i<lim;i++)if(i<rev[i])swap(f[i],f[rev[i]]);
	int bas=(kd==1)?g:((mod+1)/3);
	for(int mid=1;mid<lim;mid<<=1){
		int now=ksm(bas,(mod-1)/(mid<<1));
		for(int i=0;i<lim;i+=(mid<<1)){
			int w=1;
			for(int j=0;j<mid;j++,w=mul(w,now)){
				int a0=f[i+j],a1=mul(w,f[i+j+mid]);
				f[i+j]=add(a0,a1),f[i+j+mid]=dec(a0,a1);
			}
		}
	}
	if(kd==-1)for(int i=0,inv=ksm(lim,mod-2);i<lim;i++)f[i]=mul(f[i],inv);
}
int dep[N],val[N],ans[N],maxn,rt,mx,siz[N],son[N],adj[N],nxt[N<<1],to[N<<1],vis[N],cnt,tot;
inline void addedge(int u,int v){
	nxt[++cnt]=adj[u],adj[u]=cnt,to[cnt]=v;
}
void getrt(int u,int fa){
	siz[u]=1,son[u]=0;
	for(int e=adj[u];e;e=nxt[e]){
		int v=to[e];
		if(vis[v]||v==fa)continue;
		getrt(v,u),siz[u]+=siz[v];
		if(siz[v]>son[u])son[u]=siz[v];
	}
	son[u]=max(son[u],maxn-siz[u]);
	if(son[u]<son[rt])rt=u;
}
void getdep(int u,int fa){
	val[++tot]=dep[u],mx=max(mx,dep[u]);
	for(int e=adj[u];e;e=nxt[e]){
		int v=to[e];
		if(vis[v]||v==fa)continue;
		dep[v]=dep[u]+1;
		getdep(v,u);
	}
}
void calc(int u,int l,int f){
	dep[u]=l,mx=tot=0;
	getdep(u,0);
	lim=1,tim=0;
	while(lim<=2*mx)lim<<=1,tim++;
	for(int i=0;i<lim;i++)A[i]=0,rev[i]=(rev[i>>1]>>1)|((i&1)<<(tim-1));
	for(int i=1;i<=tot;i++)A[val[i]]++;
	ntt(A,1);
	for(int i=0;i<lim;i++)A[i]=mul(A[i],A[i]);
	ntt(A,-1);
	for(int i=0;i<=2*mx;i++)ans[i]+=f*A[i];
}
void solve(int u){
	vis[u]=1;
	calc(u,0,1);
	for(int e=adj[u];e;e=nxt[e]){
		int v=to[e];
		if(vis[v])continue;
		calc(v,1,-1),maxn=siz[v];
		getrt(v,rt=0);
		solve(rt);
	}
}
int n;
int main(){
	maxn=son[0]=n=read();
	for(int i=1;i<n;i++){
		int u=read()+1,v=read()+1;
		addedge(u,v),addedge(v,u);
	}
	getrt(1,0);
	solve(rt);
	long double res=0;
	for(int i=0;i<n;i++)res+=(long double)ans[i]/(i+1);
	printf("%.4Lf",res);
} 
  • 0
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值