hdu 4609 3-idiots（快速傅里叶FFT）

比较裸的FFT（快速傅里叶变换），也是为了这道题而去学的，厚的白书上有简单提到，不过还是推荐看算法导论，讲的很详细。

代码的话是照着别人敲的，推荐：http://www.cnblogs.com/kuangbin/archive/2013/07/24/3210565.html写的很详细。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
#define LL __int64

const double PI=acos(-1.0);

struct complex{     //实数：r实部，i虚部
    double r,i;
    complex(double rr=0,double ii=0)
    {
        r=rr;
        i=ii;
    }
    complex operator +(const complex &b)
    {
        return complex(r+b.r,i+b.i);
    }
    complex operator -(const complex &b)
    {
        return complex(r-b.r,i-b.i);
    }
    complex operator *(const complex &b)
    {
        return complex(r*b.r-i*b.i,r*b.i+i*b.r);
    }
};

void change(complex y[],int len)  //位逆序置换
{
    int i,j,k;
    for(i=1,j=len/2;i<len-1;i++)
    {
        if(i<j)
            swap(y[i],y[j]);
        k=len/2;
        while(j>=k)
        {
            j-=k;
            k/=2;
        }
        if(j<k)
            j+=k;
    }
}

void fft(complex y[],int len,int on)
{
    change(y,len);
    for(int h=2;h<=len;h<<=1)
    {
        complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
        for(int j=0;j<len;j+=h)
        {
            complex w(1,0);
            for(int k=j;k<j+h/2;k++)
            {
                complex u=y[k];
                complex t=w*y[k+h/2];
                y[k]=u+t;
                y[k+h/2]=u-t;
                w=w*wn;
            }
        }
    }
    if(on==-1)
        for(int i=0;i<len;i++)
            y[i].r/=len;
}

const int MAXN=400040;
complex x1[MAXN];
int a[MAXN/4];
LL num[MAXN];
LL sum[MAXN];

int main()
{
    int T,n;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d",&n);
        memset(num,0,sizeof(num));
        for(int i=0;i<n;i++)
        {
            scanf("%d",&a[i]);
            num[a[i]]++;
        }
//FFT
        sort(a,a+n);
        int len1=a[n-1]+1;
        int len=1;
        while(len<2*len1)
            len<<=1;
        for(int i=0;i<len1;i++)
            x1[i]=complex(num[i],0);
        for(int i=len1;i<len;i++)  //补0
            x1[i]=complex(0,0);
        fft(x1,len,1);              //求值
        for(int i=0;i<len;i++)     //乘法
            x1[i]=x1[i]*x1[i];
        fft(x1,len,-1);             //插值
//
        for(int i=0;i<len;i++)
            num[i]=(LL)(x1[i].r+0.5);
        len=2*a[n-1];
        for(int i=0;i<n;i++)
            num[a[i]+a[i]]--;
        for(int i=1;i<=len;i++)
            num[i]/=2;
        sum[0]=0;
        for(int i=1;i<=len;i++)
            sum[i]=sum[i-1]+num[i];
        LL cnt=0;
        for(int i=0;i<n;i++)
        {
            cnt+=sum[len]-sum[a[i]];
            cnt-=(LL)(n-1-i)*(n-i-2)/2;
            cnt-=(n-1);
            cnt-=(LL)(n-1-i)*(n-i-2)/2;
        }
        LL tot=(LL)n*(n-1)*(n-2)/6;
        printf("%.7f
",(double)cnt/tot);
    }
    return 0;
}