【CodeChef】PRIMEDST Prime Distance On Tree-点分治+FFT

传送门:vjudge

---

题解

休咸板题 点分治+FFT 暴力算每个质数的贡献

---

代码

老年选手敲了半个小时，CE了一次(某函数没加type)，然后就过了

#include<bits/stdc++.h>
using namespace std;
const int N=5e4+10;
typedef double db;
typedef long long ll;

int n,t[N];ll ans;
bool pri[N];int p[N],num;
int head[N],to[N<<1],nxt[N<<1],tot;

inline int pre(int n)
{
	int i,j;pri[1]=true;
	for(i=2;i<=n;++i){
		if(!pri[i]) p[++num]=i;
		for(j=1;j<=num && i*p[j]<=n;++j){
			pri[i*p[j]]=true;
			if(i%p[j]==0) break;
		}
	}
}

inline void lk(int u,int v)
{to[++tot]=v;nxt[tot]=head[u];head[u]=tot;}

namespace FFT{
const db pi=acos(-1.0);

int L,len,rv[N<<2];

struct cc{
    db r,i;
    cc(db r_=0,db i_=0):r(r_),i(i_){};
    cc operator +(const cc&ky){return cc(r+ky.r,i+ky.i);}
    cc operator -(const cc&ky){return cc(r-ky.r,i-ky.i);}
    cc operator *(const cc&ky){return cc(r*ky.r-i*ky.i,r*ky.i+i*ky.r);}
    cc operator /(int ky){return cc(r/(db)ky,i/(db)ky);}
}f[N<<2],zr;

inline void init(int n)
{
	int i;L=0;len=1;
	for(;len<=n;len<<=1) L++;
	for(i=0;i<len;++i) rv[i]=(rv[i>>1]>>1)|((i&1)<<(L-1));
}

inline void fft(cc *e,int pr){
	int i,j,k;cc x,y,pd,ori;
    for(i=0;i<len;++i) if(i<rv[i]) swap(e[i],e[rv[i]]);
    for(i=1;i<len;i<<=1){
    	ori=cc(cos(pi/i),pr*sin(pi/i));
    	for(j=0;j<len;j+=(i<<1))
    	 for(pd=cc(1,0),k=0;k<i;++k,pd=pd*ori){
    	 	x=e[j+k];y=pd*e[i+j+k];
    	 	e[j+k]=x+y;e[i+j+k]=x-y;
		 }
	}
	if(pr==1) return;
	for(i=0;i<len;++i) e[i]=e[i]/len;
}

inline void cal(int n,int pr)
{
	init((n<<1));int i;ll j;
	for(i=0;i<=n;++i) f[i]=cc((db)t[i],0);
	for(i=n+1;i<len;++i) f[i]=zr;
	fft(f,1);
	for(i=0;i<len;++i) f[i]=f[i]*f[i];
	fft(f,-1);n<<=1;
	for(i=2;i<=n;++i) if(!pri[i]){
		j=(ll)(f[i].r+0.5);if(!(i&1)) j-=t[i>>1];
		ans+=pr*(j>>1);
	}
}
}

namespace DFZ{
int S,MN,sz[N],mx[N],rt;
int dep[N],mxdep;
bool vs[N];

void fdrt(int x,int fr)
{
	sz[x]=1;mx[x]=0;int i,j;
	for(i=head[x];i;i=nxt[i]){
		j=to[i];if(vs[j] || j==fr) continue;
		fdrt(j,x);sz[x]+=sz[j];
		mx[x]=max(mx[x],sz[j]);
	}
	mx[x]=max(mx[x],S-sz[x]);
	if(mx[x]<MN) MN=mx[x],rt=x;
}

void dfs(int x,int fr)
{
	int i,j;mxdep=max(mxdep,dep[x]);t[dep[x]]++;
	for(i=head[x];i;i=nxt[i]){
		j=to[i];if(j==fr || vs[j]) continue;
		dep[j]=dep[x]+1;dfs(j,x);
	}
}

void sol(int x)
{
	int i,j,sum=S;vs[x]=true;
	dep[x]=mxdep=0;dfs(x,0);
	for(i=2;i<=mxdep;++i) if(!pri[i]) ans+=t[i];
	t[0]=0;FFT::cal(mxdep,1);
	memset(t,0,(mxdep+1)<<2);
	for(i=head[x];i;i=nxt[i]){
		j=to[i];if(vs[j]) continue;
		dep[j]=mxdep=1;dfs(j,x);
		FFT::cal(mxdep,-1);
		memset(t,0,(mxdep+1)<<2);
	}
	for(i=head[x];i;i=nxt[i]){
		j=to[i];if(vs[j]) continue;
		S=(sz[j]>sz[x])?(sum-sz[x]):sz[j];MN=N;
		fdrt(j,x);sol(rt);
	}
}
}
using namespace DFZ;
int main(){
	int i,j,x,y;
	scanf("%d",&n);pre(n);
	for(i=1;i<n;++i){
		scanf("%d%d",&x,&y);
		lk(x,y);lk(y,x);
	}
	S=n;MN=n;fdrt(1,0);sol(rt);
	printf("%.8lf",ans/(db)((ll)n*(n-1)/2));
	return 0;
}