【XSY2332】Randomized Binary Search Tree 概率DP FFT

题目描述

　　(forall 0leq i<n)，求有多少棵(n)个点，权值和优先级完全随机的treap的树高为(i)。

　　(nleq 30000)

题解

　　设(f_{i,j})为(j)个点的树，树高不超过为(i)的概率

[f_{i,j}=frac{1}{j}sum_{k=1}^{j}f_{i-1,j-1} imes f_{i-1,j-k} ]

　　枚举一个点左子树大小(k-1)，那么右子树大小为(j-k)。且这个点的优先级为这(j)个点最小的概率是(frac{1}{j})。

　　这个东西是个卷积，可以用FFT加速。

　　其实期望树高是(O(log n))的。实际上只有前面一部分的答案不为(0)。所以我们只用计算树高(leq 100)的答案。

　　时间复杂度：(O(nlog^2n))

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
#include<utility>
#include<cmath>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
double pi=acos(-1);
struct cp
{
	double x,y;
	cp(double a=0,double b=0)
	{
		x=a;
		y=b;
	}
};
cp operator +(cp a,cp b)
{
	return cp(a.x+b.x,a.y+b.y);
}
cp operator -(cp a,cp b)
{
	return cp(a.x-b.x,a.y-b.y);
}
cp operator *(cp a,cp b)
{
	return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
cp operator /(cp a,double b)
{
	return cp(a.x/b,a.y/b);
}
namespace fft
{
	cp w1[100010];
	cp w2[100010];
	int rev[100010];
	int n;
	void init(int m)
	{
		n=1;
		while(n<=m)
			n<<=1;
		int i;
		for(i=2;i<=n;i<<=1)
		{
			w1[i]=cp(cos(2*pi/i),sin(2*pi/i));
			w2[i]=cp(cos(2*pi/i),-sin(2*pi/i));
		}
		rev[0]=0;
		for(i=1;i<n;i++)
			rev[i]=(rev[i>>1]>>1)|(i&1?n>>1:0);
	}
	void fft(cp *a,int t)
	{
		int i,j,k;
		cp u,v,w,wn;
		for(i=0;i<n;i++)
			if(rev[i]<i)
				swap(a[i],a[rev[i]]);
		for(i=2;i<=n;i<<=1)
		{
			wn=(~t?w1[i]:w2[i]);
			for(j=0;j<n;j+=i)
			{
				w=1;
				for(k=j;k<j+i/2;k++)
				{
					u=a[k];
					v=a[k+i/2]*w;
					a[k]=u+v;
					a[k+i/2]=u-v;
					w=w*wn;
				}
			}
		}
		if(t==-1)
			for(i=0;i<n;i++)
				a[i]=a[i]/n;
	}
}
cp a[100010];
double f[110][30010];
double ans[30010];
int main()
{
	freopen("b.in","r",stdin);
	freopen("b.out","w",stdout);
	int n;
	scanf("%d",&n);
	int m=100;
	int i,j;
	fft::init(2*n);
	a[0]=a[1]=1;
	double last;
	for(i=0;i<n;i++)
		ans[i]=0;
	last=ans[0]=a[n].x;
	for(i=1;i<=m;i++)
	{
		fft::fft(a,1);
		for(j=0;j<fft::n;j++)
			a[j]=a[j]*a[j];
		fft::fft(a,-1);
		for(j=fft::n-1;j>=1;j--)
			a[j]=a[j-1];
		a[0]=a[1]=1;
		for(j=2;j<=n;j++)
			a[j]=a[j]/j;
		for(j=n+1;j<fft::n;j++)
			a[j]=0;
		ans[i]=a[n].x-last;
		last=a[n].x;
	}
//	for(i=0;i<=m;i++)
//		f[i][1]=f[i][0]=1;
//	for(i=1;i<=m;i++)
//		for(j=2;j<=n;j++)
//			for(k=1;k<=j;k++)
//				f[i][j]+=f[i-1][k-1]*f[i-1][j-k]/j;
//	for(i=0;i<=m;i++)
//	{
//		ans[i]=f[i][n];
//		if(i>=1)
//			ans[i]-=f[i-1][n];
//	}
	for(i=0;i<=n-1;i++)
		printf("%.10lf
",ans[i]);
	return 0;
}