HDU 4609 3-idiots FFT

http://acm.hdu.edu.cn/showproblem.php?pid=4609
King OMeGa catched three men who had been streaking in the street. Looking as idiots though, the three men insisted that it was a kind of performance art, and begged the king to free them. Out of hatred to the real idiots, the king wanted to check if they were lying. The three men were sent to the king’s forest, and each of them was asked to pick a branch one after another. If the three branches they bring back can form a triangle, their math ability would save them. Otherwise, they would be sent into jail.
However, the three men were exactly idiots, and what they would do is only to pick the branches randomly. Certainly, they couldn’t pick the same branch - but the one with the same length as another is available. Given the lengths of all branches in the forest, determine the probability that they would be saved.
Input
An integer T(T≤100) will exist in the first line of input, indicating the number of test cases.
Each test case begins with the number of branches N(3≤N≤10 5).
The following line contains N integers a_i (1≤a_i≤10 5), which denotes the length of each branch, respectively.
Output
Output the probability that their branches can form a triangle, in accuracy of 7 decimal places.
Sample Input
2
4
1 3 3 4
4
2 3 3 4
Sample Output
0.5000000
1.0000000

题目大意：给出$n$条线段的长度，问从中任取三条可以组成三角形的概率(每条线段只能取一次)。
思路：$FFT$，需要用到很巧妙的转换，推荐看$kuangbin$大佬的博客，讲的超级棒。https://www.cnblogs.com/kuangbin/archive/2013/07/24/3210565.html

#include<bits/stdc++.h>
using namespace std;        //FFT模板
typedef long long ll;

const int maxn=1e5+5;

struct Complex
{
    double x,y;
    Complex(double dx=0,double dy=0)
    {
        x=dx;
        y=dy;
    }
};

Complex operator +(Complex a,Complex b)
{
    return Complex(a.x+b.x,a.y+b.y);
}
Complex operator -(Complex a,Complex b)
{
    return Complex(a.x-b.x,a.y-b.y);
}
Complex operator *(Complex a,Complex b)
{
    return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}

const double pi=acos(-1.0); //PI
int limit=1,bit=0;
int wz[maxn<<2];
Complex a[maxn<<2];
int v[maxn],cnt[maxn];
ll num[maxn<<2],sum[maxn<<2];

void FFT(Complex *A,int inv)
{
    for(int i=0;i<limit;i++)
        if(i<wz[i])
            swap(A[i],A[wz[i]]);
    for(int mid=1;mid<limit;mid<<=1)
    {
        Complex wn(cos(pi/mid),inv*sin(pi/mid));
        for(int i=0;i<limit;i+=mid<<1)
        {
            Complex w(1,0);
            for(int j=0;j<mid;j++,w=w*wn)
            {
                Complex t1=A[i+j];
                Complex t2=w*A[i+mid+j];
                A[i+j]=t1+t2;
                A[i+mid+j]=t1-t2;
            }
        }
    }
}

int main()
{
    int t,n;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d",&n);
        memset(a,0,sizeof(a));
        memset(cnt,0,sizeof(cnt));
        limit=1,bit=0;
        for(int i=0;i<n;i++)
        {
            scanf("%d",&v[i]);
            cnt[v[i]]++;
        }
        sort(v,v+n);
        int len=v[n-1];
        int tmp=len;
        len<<=1;
        while(limit<=len)
        {
            limit<<=1;
            bit++;
        }
        for(int i=0;i<limit;i++) //从低位到高位
        {
            if(i<=tmp)
                a[i].x=cnt[i];
            else
                a[i].x=0;
        }
        for(int i=0;i<limit;i++)
            wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
        FFT(a,1);
        for(int i=0;i<limit;i++)
            a[i]=a[i]*a[i];
        FFT(a,-1);
        for(int i=0;i<limit;i++)
            num[i]=(ll)(a[i].x/limit+0.5);
        for(int i=0;i<n;i++)
            num[v[i]<<1]--; //选取同一根木棒两次 不合题意
        for(int i=0;i<limit;i++)
            num[i]>>=1; //选择是无序的 所以要除2
        for(int i=1;i<limit;i++)
            sum[i]=sum[i-1]+num[i];
        ll ans=0;
        for(int i=0;i<n;i++)//枚举三角形的最长边 v_k
        {
            ans+=sum[limit-1]-sum[v[i]]; // v_i + v_j > v_k 满足三角形性质
            ans-=(ll)(n-i-1)*i; // 一大一小 不合题意
            ans-=n-1;   // v_i || v_j 与 v_k 是同一根 不合题意
            ans-=(ll)(n-i-1)*(n-i-2)/2; //v_i >= v_k && v_j >= v_k 不合题意
        }
        ll tot=(ll)n*(n-1)*(n-2)/6;
        printf("%.7f\n",ans*1.0/tot);
    }
    return 0;
}