HDU-4609 3-idiots FFT

　　题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=4609　

　　题意：给n个棒子，求任意组成3根能形成三角形的概率。

　　用FFT求出任意两根棒子组合成新的长度，这种长度的组合有多少种。

　　用cnt[i]表示和的长度为i的棒子组合有多少种。

　　(1) 首先去掉cnt里重复的部分，一根棒子不能与自己组合，所以有cnt[a[i] * 2] --

　　(2) 任意两个棒子组合的顺序我们不需要考虑，而实际上他们算了两次，所以有cnt[i] /=2

　　(3) 枚举所有棒子，假设这根棒子a[i]是组成三角形其中最长的那根，我们首先求出另外两根的长度和比它的长的总种数，即sigma(cnt[j])   (a[i]<j<=crest*2，crest是最长的棒子长度)

　　(4) 另外两根中，这根被枚举的棒子不能再出现了，所以要减去n-1

　　(5) 另外两根中，可能有一根长度大于a[i]，另一根小于a[i]，所以要减去(n-i) * (i-1)

　　(6) 另外两根中，可能两根长度都大于a[i]，再减去(n-i) * (n-i-1) / 2

//STATUS:C++_AC_2140MS_7128KB
#include <functional>
#include <algorithm>
#include <iostream>
//#include <ext/rope>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <numeric>
#include <cstring>
#include <cassert>
#include <cstdio>
#include <string>
#include <vector>
#include <bitset>
#include <queue>
#include <stack>
#include <cmath>
#include <ctime>
#include <list>
#include <set>
#include <map>
using namespace std;
//#pragma comment(linker,"/STACK:102400000,102400000")
//using namespace __gnu_cxx;
//define
#define pii pair<int,int>
#define mem(a,b) memset(a,b,sizeof(a))
#define lson l,mid,rt<<1
#define rson mid+1,r,rt<<1|1
#define PI acos(-1.0)
//typedef
typedef __int64 LL;
typedef unsigned __int64 ULL;
//const
const int N=400010;
const int INF=0x3f3f3f3f;
const int MOD=10007,STA=8000010;
const LL LNF=1LL<<55;
const double EPS=1e-4;
const double OO=1e30;
const int dx[4]={-1,0,1,0};
const int dy[4]={0,1,0,-1};
const int day[13]={0,31,28,31,30,31,30,31,31,30,31,30,31};
//Daily Use ...
inline int sign(double x){return (x>EPS)-(x<-EPS);}
template<class T> T gcd(T a,T b){return b?gcd(b,a%b):a;}
template<class T> T lcm(T a,T b){return a/gcd(a,b)*b;}
template<class T> inline T lcm(T a,T b,T d){return a/d*b;}
template<class T> inline T Min(T a,T b){return a<b?a:b;}
template<class T> inline T Max(T a,T b){return a>b?a:b;}
template<class T> inline T Min(T a,T b,T c){return min(min(a, b),c);}
template<class T> inline T Max(T a,T b,T c){return max(max(a, b),c);}
template<class T> inline T Min(T a,T b,T c,T d){return min(min(a, b),min(c,d));}
template<class T> inline T Max(T a,T b,T c,T d){return max(max(a, b),max(c,d));}
//End
//复数结构体
struct complex
{
    double r,i;
    complex(double _r = 0.0,double _i = 0.0)
    {
        r = _r; i = _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);
    }
    complex operator *(const complex &b)
    {
        return complex(r*b.r-i*b.i,r*b.i+i*b.r);
    }
};
/*
 * 进行FFT和IFFT前的反转变换。
 * 位置i和 （i二进制反转后位置）互换
 * len必须去2的幂
 */
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]);
        //交换互为小标反转的元素，i<j保证交换一次
        //i做正常的+1，j左反转类型的+1,始终保持i和j是反转的
        k = len/2;
        while( j >= k)
        {
            j -= k;
            k /= 2;
        }
        if(j < k) j += k;
    }
}
/*
 * 做FFT
 * len必须为2^k形式，
 * on==1时是DFT，on==-1时是IDFT
 */
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;
}

LL sum[N/2],tot[N/2];
int cnt[N/4],l[N/4];
complex a[N],b[N];
int T,n;

int main(){
 //   freopen("in.txt","r",stdin);
    int i,j,len,hig;
    LL ans;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d",&n);
        mem(cnt,0);hig=0;
        for(i=1;i<=n;i++){
            scanf("%d",&l[i]);
            cnt[l[i]]++;
            hig=Max(hig,l[i]);
        }
        sort(l+1,l+n+1);
        hig++;
        for(len=1;len<(hig<<1);len<<=1);
        for(i=0;i<hig;i++)
            a[i]=b[i]=complex(cnt[i],0);
        for(;i<=len;i++)
            a[i]=b[i]=complex(0,0);

        FFT(a,len,1);
        FFT(b,len,1);
        for(i=0;i<=len;i++)
            a[i]=a[i]*b[i];
        FFT(a,len,-1);
        len=(hig<<1)-1;
        for(i=0;i<=len;i++)
            tot[i]=(LL)(a[i].r+0.5);
        for(i=1;i<=n;i++)
            tot[l[i]<<1]--;
        for(i=1;i<=len;i++)tot[i]>>=1;
        sum[0]=0;
        for(i=1;i<=len;i++)
            sum[i]=tot[i]+sum[i-1];
        ans=0;
        for(i=1;i<=n;i++){
            ans+=sum[len]-sum[l[i]];
            ans-=(LL)(n-1)+(LL)(i-1)*(n-i)+(LL)(n-i)*(n-i-1)/2;
        }

        printf("%.7lf
",(double)ans/n/(n-1)/(n-2)*6);
    }
    return 0;
}