您的位置:首页 > 其它

Tyvj1953:Normal (点分治+FFT)

2017-12-14 12:57 281 查看
题目传送门:http://tyvj.joyoi.cn/p/1953

题目分析:好神的一道题,%了一波dalao的题解才会做,然后我发现我对期望一无所知QAQ……

这题的难点其实既不在于点分,也不在于FFT,而是对期望的转化。首先由于期望的线性性,我们可以单独计算每一个点x对期望时间的贡献。而点x每被计算一次,就意味着存在一个分治中心y(x可以等于y),使得x到y路径上的其它点都还没有被作为分治中心。由于x到y路径上的每个点,作为这条路径第一个分治中心的概率是相等的,所以如果设dis(x,y)表示x到y路径上的点数,那么有:

ans=∑x=1n∑y=1n1dis(x,y)

又由于dis(x,y)在分母处,所以不能直接用裸的点分治求答案。不妨记num[x]为长度是x的路径条数,然后用点分治+FFT算num数组,最后再统计答案。注意点分治要用容斥的写法,即算完当前连通块的答案再递归子连通块减去非法答案。如果直接将所有子树逐个用FFT卷起来,是没有办法保证时间复杂度的。我不会告诉你tututu用先卷深度小的子树,再卷深度大的子树的方法卡过了

CODE:

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
using namespace std;

const int maxn=300100;
const double pi=acos(-1.0);

struct Complex
{
double X,Y;
Complex (double a=0.0,double b=0.0) : X(a),Y(b) {}
} ;

Complex operator+(Complex a,Complex b){return Complex(a.X+b.X,a.Y+b.Y);}
Complex operator-(Complex a,Complex b){return Complex(a.X-b.X,a.Y-b.Y);}
Complex operator*(Complex a,Complex b){return Complex(a.X*b.X-a.Y*b.Y,a.X*b.Y+a.Y*b.X);}

int Rev[maxn];
Complex A[maxn];

struct edge
{
int obj;
edge *Next;
} e[maxn<<1];
edge *head[maxn];
int cur=-1;

int Size[maxn];
int Max[maxn];
int now[maxn];
int tail;

int val[maxn];
bool vis[maxn];
int num[maxn];
int n,N,Lg;

void Add(int x,int y)
{
cur++;
e[cur].obj=y;
e[cur].Next=head[x];
head[x]=e+cur;
}

void Dfs1(int node,int fa)
{
now[++tail]=node;
Size[node]=1;
Max[node]=0;
for (edge *p=head[node]; p; p=p->Next)
{
int to=p->obj;
if ( !vis[to] && to!=fa )
{
Dfs1(to,node);
Size[node]+=Size[to];
Max[node]=max(Max[node],Size[to]);
}
}
}

void Dfs2(int node,int fa,int dep)
{
val[dep]++;
for (edge *p=head[node]; p; p=p->Next)
{
int to=p->obj;
if ( !vis[to] && to!=fa ) Dfs2(to,node,dep+1);
}
}

void DFT(Complex *a,double f)
{
for (int i=0; i<N; i++)
if (i<Rev[i]) swap(a[i],a[ Rev[i] ]);
for (int len=2; len<=N; len<<=1)
{
int mid=len>>1;
double ang=2.0*pi/((double)len);
Complex e( cos(ang) , f*sin(ang) );
for (Complex *p=a; p!=a+N; p+=len)
{
Complex wn(1.0,0.0);
for (int i=0; i<mid; i++)
{
Complex temp=wn*p[mid+i];
p[mid+i]=p[i]-temp;
p[i]=p[i]+temp;
wn=wn*e;
}
}
}
}

void FFT(int m,int f,int x)
{
N=1,Lg=0;
while (N<2*m) N<<=1,Lg++;
for (int i=0; i<N; i++)
{
Rev[i]=0;
for (int j=0; j<Lg; j++)
if ( i&(1<<j) ) Rev[i]|=( 1<<(Lg-j-1) );
}
for (int i=0; i<m; i++) A[i]=Complex((double)val[i+1],0.0);
for (int i=m; i<N; i++) A[i]=Complex(0.0,0.0);
DFT(A,1.0);
for (int i=0; i<N; i++) A[i]=A[i]*A[i];
DFT(A,-1.0);
for (int i=0; i<N; i++) A[i].X/=((double)N);
for (int i=0; i<N; i++) num[i+1+x]+=( f*(int)floor( A[i].X+0.5 ) );
}

void Dec(int node)
{
tail=0;
Dfs1(node,node);
for (int i=1; i<=tail; i++) val[i]=0;
Dfs2(node,node,1);
FFT(tail,-1,2);
}

void Solve(int node)
{
tail=0;
Dfs1(node,node);
if (tail==1)
{
num[1]++;
return;
}
int root=node,sum=Size[node];
for (int i=2; i<=tail; i++)
{
int x=now[i];
if ( max(Max[x],sum-Size[x])<max(Max[root],sum-Size[root]) ) root=x;
}
for (int i=1; i<=sum; i++) val[i]=0;
Dfs2(root,root,1);
FFT(sum,1,0);
vis[root]=true;
for (edge *p=head[root]; p; p=p->Next)
{
int to=p->obj;
if (!vis[to]) Dec(to),Solve(to);
}
}

int main()
{
freopen("1953.in","r",stdin);
freopen("1953.out","w",stdout);

scanf("%d",&n);
for (int i=1; i<=n; i++) head[i]=NULL;
for (int i=1; i<n; i++)
{
int x,y;
scanf("%d%d",&x,&y);
x++,y++;
Add(x,y);
Add(y,x);
}

Solve(1);

double ans=0.0;
for (int i=1; i<=n; i++) ans+=( (double)num[i]/(double)i );
printf("%.4lf\n",ans);

return 0;
}

//FFT重载运算符的时候要仔细检查!!!
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: