【hdu4609】 3-idiots FFT

题外话：好久没写blog了啊～～

题目传送门

题目大意：给你m条长度为ai的线段，求在其中任选三条出来，能构成三角形的概率。即求在这n条线段中找出三条线段所能拼出的三角形数量除以$\binom{m}{3}$。

假设我们手中有3条长度分别为$x,y,z$的边(为了简化问题我们假设$x<y<z$，$x,y,z$相等的情况另行讨论)，如果他们能拼成三角形，必然满足$x+y>z$且$z-y<x$。

该题的$O(m^3)$做法：枚举其中的3条边，套用上面的判断公式，进行累计。

但通过简单的变式，我们假设我们已经确定了$x$和$y$，那么$z$的范围即为$[y,x+y)$，我们维护一个数组$num$，$num_i$表示长度为i的线段数量。再维护一个num的前缀和sum。则$ans=\sum^{n-1}_{x=1}  \sum^{n}_{y=x+1} num_x * num_y*(sum[x+y-1]-sum[y])$，运算该和式的时间复杂度为$O(n^2)$，此处的n表示最长线段的长度。

然而还是会TLE。。。。

我们考虑对其做一些变式。

$ans=\sum^{n-1}_{x=1}  \sum^{n}_{y=x+1} num_x * num_y*(sum[x+y-1]-sum[y])$

 $=\sum^{n-1}_{x=1}  \sum^{n}_{y=x+1} num_x * num_y*sum[x+y-1] -\sum^{n-1}_{x=1}  \sum^{n}_{y=x+1} num_x * num_y*sum[y]$

该式子的后半部分，我们可以通过维护$num_y*sum[y]$的后缀和，实现O(n)的计算。

下面我们继续对式子的前半部分变式。令t=x+y。则有

$=\frac{1}{2}\sum^{2n}_{t=2}  \sum^{t-1}_{p=1} num_p*num_{t-p}*sum[t-1] -\frac{1}{2}\sum^{n}_{i=1}num_i^2*sum[2*i-1]$

我们发现，该式子的$\sum_{p=1}^{t-1} num_{p}\times num_{t-p}$，可以用FFT求出。则时间复杂度成功降低至O(n log n)，而式子的后半部分可以O(n)求出，求和时间复杂度降低至O(n log n)。

下面说下x=y=z，x=y<z，x<y=z的处理方法。

x=y=z:$\sum _{\forall num_x>2} \binom{num_x}{3}$

x=y<z: $\sum _{x=1}^{n} \binom{num_x}{2}*(sum[2x-1]-sum[x])$

x<y=z:$\sum_{x=1}^{n}(sumc2_n-sumc2_x)$ ，其中$sumc2_x$ 表示$\sum_{i=1}^{x}\binom{num_x}{2}$。

把这三种情况和最初描述的情况相加即可。

#include<bits/stdc++.h>
#define M 270000
#define cp complex<double>
#define PI acos(-1)
#define L long long
using namespace std;
cp a[M];
void change(cp a[],int len){
    for(int i=0,j=0;i<len-1;i++){
        if(i<j) swap(a[i],a[j]);
        int k=len>>1;
        while(j>=k) j-=k,k>>=1;
        j+=k;
    }
}
void fft(cp a[],int n,int on){
    change(a,n); cp t,u;
    for(int h=2;h<=n;h<<=1){
        cp wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
        for(int j=0;j<n;j+=h){
            cp w(1,0);
            for(int k=j;k<j+(h>>1);k++){
                u=a[k]; t=w*a[k+(h>>1)];
                a[k]=u+t; a[k+(h>>1)]=u-t;
                w=w*wn;
            }
        }
    }
}
L num[M]={0},sum[M]={0},sumhh[M]={0},sumc2[M]={0},sc[M]={0},ans=0;

int Main(){
    int m,n,maxn=0; scanf("%d",&n);
    for(int i=1;i<=n;i++){
        int x; scanf("%d",&x);
        num[x]++; maxn=max(maxn,x);
    }
    for(m=1;m<(maxn*2+1);m<<=1);
    for(int i=0;i<m;i++) a[i]=cp(num[i],0);
    fft(a,m,1);
    for(int i=0;i<m;i++) a[i]=a[i]*a[i];
    fft(a,m,-1);
    for(int i=0;i<m;i++) sc[i]=(a[i+1].real()+0.5)/m;
    for(int i=1;i<m;i++){
        sum[i]=sum[i-1]+num[i];
        sumc2[i]=sumc2[i-1]+num[i]*(num[i]-1)/2;
    }
    
    L ans1=0,ans2=0,ans3=0,ans4=0;
    for(int i=1;i<=maxn;i++) if(num[i]>2) 
    ans1+=(L)num[i]*(num[i]-1)*(num[i]-2);
    ans1/=6;//x=y=z
    for(int i=1;i<=maxn;i++) 
    ans2+=(L)num[i]*(num[i]-1)*(sum[2*i-1]-sum[i]);
    ans2/=2;//x=y<z
    for(int i=1;i<=maxn;i++)
    ans3+=num[i]*(sumc2[maxn]-sumc2[i]);
    //x<y=z
    for(int i=1;i<2*maxn;i++) ans4+=sum[i]*sc[i];
    for(int i=1;i<=maxn;i++) ans4-=num[i]*num[i]*sum[2*i-1];
    ans4/=2;//x<y<z卷积部分
    
    for(int i=maxn;i;i--) sumhh[i]=sumhh[i+1]+num[i]*sum[i];
    for(int i=1;i<maxn;i++) ans4-=num[i]*sumhh[i+1];
    ans=ans1+ans2+ans3+ans4;
    
    double fenmu=(L)n*(n-1)*(n-2)/6;
    double hh=ans/fenmu;
    printf("%.7lf\n",hh);
}
int main(){
    //freopen("in.txt","r",stdin);
    //freopen("out.txt","w",stdout);
    int cas; scanf("%d",&cas);
    while(cas--){
        memset(num,0,sizeof(num)); memset(sum,0,sizeof(sum));
        memset(sumhh,0,sizeof(sumhh)); memset(sumc2,0,sizeof(sumc2));
        memset(sc,0,sizeof(sc)); memset(a,0,sizeof(a)); ans=0;
        Main();
    }
}