[BZOJ3451][Tyvj1953]Normal(点分治+FFT)
https://www.cnblogs.com/GXZlegend/p/8611948.html
1 #include<cmath> 2 #include<cstdio> 3 #include<algorithm> 4 #define rep(i,l,r) for (int i=(l); i<=(r); i++) 5 #define For(i,x) for (int i=h[x],k; i; i=nxt[i]) 6 using namespace std; 7 8 const int N=100010; 9 const double pi=acos(-1.); 10 int n,u,v,S,rt,tot,cnt,mx,d[N],q[N],rev[N]; 11 int vis[N],sz[N],f[N],h[N],to[N<<1],nxt[N<<1]; 12 double ans,num[N]; 13 void add(int u,int v){ to[++cnt]=v; nxt[cnt]=h[u]; h[u]=cnt; } 14 15 struct C{ double x,y; }a[N]; 16 C operator +(const C &a,const C &b){ return (C){a.x+b.x,a.y+b.y}; } 17 C operator -(const C &a,const C &b){ return (C){a.x-b.x,a.y-b.y}; } 18 C operator *(const C &a,const C &b){ return (C){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x}; } 19 20 void find(int x,int fa){ 21 sz[x]=1; f[x]=0; 22 For(i,x) if (!vis[k=to[i]] && k!=fa) 23 find(k,x),sz[x]+=sz[k],f[x]=max(f[x],sz[k]); 24 f[x]=max(f[x],S-sz[x]); 25 if (f[x]<f[rt]) rt=x; 26 } 27 28 void get(int x,int fa){ 29 q[++tot]=d[x]; 30 For(i,x) if (!vis[k=to[i]] && k!=fa) d[k]=d[x]+1,get(k,x); 31 } 32 33 void DFT(C a[],int n,int f){ 34 for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]); 35 for (int i=1; i<n; i<<=1){ 36 C wn=(C){cos(pi/i),f*sin(pi/i)}; 37 for (int p=i<<1,j=0; j<n; j+=p){ 38 C w=(C){1,0}; 39 for (int k=0; k<i; k++,w=w*wn){ 40 C x=a[j+k],y=w*a[i+j+k]; a[j+k]=x+y; a[i+j+k]=x-y; 41 } 42 } 43 } 44 if (f==-1) for (int i=0; i<n; i++) a[i].x/=n; 45 } 46 47 void calc(int fl){ 48 int n=1,L=0; mx=0; 49 rep(i,1,tot) mx=max(mx,q[i]); 50 for (; n<=mx<<1; n<<=1) L++; 51 for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1)); 52 for (int i=0; i<n; i++) a[i].x=a[i].y=0; 53 rep(i,1,tot) a[q[i]].x++; 54 DFT(a,n,1); 55 for (int i=0; i<n; i++) a[i]=a[i]*a[i]; 56 DFT(a,n,-1); 57 rep(i,0,2*mx) num[i]+=fl*(int)(a[i].x+0.1); 58 } 59 60 void work(int x){ 61 vis[x]=1; d[x]=tot=0; get(x,0); calc(1); 62 For(i,x) if (!vis[k=to[i]]){ 63 d[k]=1; tot=0; get(k,0); calc(-1); 64 S=sz[k]; rt=0; find(k,0); work(rt); 65 } 66 } 67 68 int main(){ 69 freopen("bzoj3451.in","r",stdin); 70 freopen("bzoj3451.out","w",stdout); 71 scanf("%d",&n); 72 rep(i,2,n) scanf("%d%d",&u,&v),add(u+1,v+1),add(v+1,u+1); 73 S=f[0]=n; find(1,0); work(rt); 74 rep(i,1,n) ans+=num[i-1]/i; 75 printf("%.4lf\n",ans); 76 return 0; 77 }