【XSY2332】Randomized Binary Search Tree 概率DP FFT

题目描述

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

　　$n\leq 30000$

题解

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

$$f_{i,j}=\frac{1}{j}\sum_{k=1}^{j}f_{i-1,j-1}\times f_{i-1,j-k}$$

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

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

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

　　时间复杂度：$O(n{log}^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\n",ans[i]);
    return 0;
}