【BZOJ3451】Normal【期望线性性】【点分治】【NTT卷积】

题意:随机分治中心点分治的期望操作次数

n ≤ 3 × 1 0 4 n\leq 3\times 10^4 n3×104

即求点分树的 siz 之和的期望

即祖孙关系对数期望

考虑一有序点对 ( u , v ) (u,v) (u,v) u u u 在点分树上是 v v v 祖先当且仅当 u u u u ∼ v u\sim v uv 路径上第一个被选为分治中心的,并且选择路径外的点是不影响的,所以概率为 1 dist ⁡ ( u , v ) + 1 \dfrac 1{\operatorname{dist}(u,v)+1} dist(u,v)+11

由期望线性性,答案为

∑ i = 1 n ∑ j = 1 n 1 dist ⁡ ( i , j ) + 1 \sum_{i=1}^n\sum_{j=1}^{n}\frac 1{\operatorname{dist}(i,j)+1} i=1nj=1ndist(i,j)+11

相当于是求距离为 0 ∼ n − 1 0\sim n-1 0n1 的点对数量,考虑点分治

计算当前答案时发现是个卷积形式,使用 NTT 优化

算出整个连通块的答案,然后容斥掉在同一个子树中的

因为深度是不会超过大小的,所以复杂度为 O ( n log ⁡ 2 n ) O(n\log ^2n) O(nlog2n)

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#include <vector>
#define MAXN 65536
using namespace std;
const int MOD=998244353;
inline int add(const int& x,const int& y){return x+y>=MOD? x+y-MOD:x+y;}
inline int dec(const int& x,const int& y){return x<y? x-y+MOD:x-y;}
typedef long long ll;
inline int qpow(int a,int p)
{
	int ans=1;
	while (p)
	{
		if (p&1) ans=(ll)ans*a%MOD;
		a=(ll)a*a%MOD,p>>=1;
	}
	return ans;
}
#define inv(x) qpow(x,MOD-2)
int rt[2][24];
int l,lim,r[MAXN];
inline void init(){lim=1<<l;for (int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));}
inline void NTT(int* a,int type)
{
	for (int i=0;i<lim;i++) if (i<r[i]) swap(a[i],a[r[i]]);
	for (int L=0;L<l;L++)
	{
		int mid=1<<L,len=mid<<1;
		ll Wn=rt[type][L+1];
		for (int s=0;s<lim;s+=len)
			for (int k=0,w=1;k<mid;k++,w=w*Wn%MOD)
			{
				int x=a[s+k],y=(ll)w*a[s+mid+k]%MOD;
				a[s+k]=add(x,y),a[s+mid+k]=dec(x,y);
			}
	}
	if (type)
	{
		int t=inv(lim);
		for (int i=0;i<lim;i++) a[i]=(ll)a[i]*t%MOD;
	}
}
inline double calc(int* a,int n)
{
	double ans=0;
	for (l=0;(1<<l)<=(n<<1);++l);
	init(),NTT(a,0);
	for (int i=0;i<lim;i++) a[i]=(ll)a[i]*a[i]%MOD;
	NTT(a,1);
	for (int i=0;i<lim;i++)
	{
		ans+=(double)a[i]/(i+1);
		a[i]=0;
	}
	return ans;
}
vector<int> e[MAXN];
double ans;
bool vis[MAXN];
int siz[MAXN],maxp[MAXN]={0x7fffffff},Rt;
int a[MAXN];
void findrt(int u,int f,int sum)
{
	siz[u]=1,maxp[u]=0;
	for (int i=0;i<(int)e[u].size();i++)
		if (e[u][i]!=f&&!vis[e[u][i]])
		{
			findrt(e[u][i],u,sum);
			siz[u]+=siz[e[u][i]];
			maxp[u]=max(maxp[u],siz[e[u][i]]);
		}
	maxp[u]=max(maxp[u],sum-siz[u]);
	if (maxp[u]<maxp[Rt]) Rt=u;
}
int dfs(int u,int f,int dep)
{
	siz[u]=1,++a[dep];
	int mx=dep;
	for (int i=0;i<(int)e[u].size();i++)
		if (e[u][i]!=f&&!vis[e[u][i]])
		{
			mx=max(mx,dfs(e[u][i],u,dep+1));
			siz[u]+=siz[e[u][i]];
		}
	return mx;
}
inline void calc()
{
	ans+=calc(a,dfs(Rt,0,0));
	for (int i=0;i<(int)e[Rt].size();i++)
		if (!vis[e[Rt][i]])
			ans-=calc(a,dfs(e[Rt][i],0,1));
}
void solve()
{
	int u=Rt;
	vis[u]=1;
	calc();
	for (int i=0;i<(int)e[u].size();i++)
		if (!vis[e[u][i]])
		{
			int v=e[u][i];
			Rt=0,findrt(v,0,siz[v]);
			solve();
		}
}
int main()
{
	rt[0][23]=qpow(3,119),rt[1][23]=inv(rt[0][23]);
	for (int i=22;i>=0;i--)
	{
		rt[0][i]=(ll)rt[0][i+1]*rt[0][i+1]%MOD;
		rt[1][i]=(ll)rt[1][i+1]*rt[1][i+1]%MOD;
	}
	int n;
	scanf("%d",&n);
	for (int i=1;i<n;i++)
	{
		int u,v;
		scanf("%d%d",&u,&v);
		++u,++v;
		e[u].push_back(v),e[v].push_back(u);
	}
	findrt(1,0,n);
	solve();
	printf("%.4f\n",ans);
	return 0;
} 
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

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

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值