HDU 4609 3-idiots（FFT,计数）

题目链接
题意：给出n个长度，问任意选择出来3个能形成三角形的概率

FFT学习中，由卷积的定义,设一个多项式A(x)，其中x^k的系数ak 就是表示长度是k的个数，那么B(x)=A(x)*A(x)，B(x)的x^k系数bk就是A(x)系数对应的集合中，可以重复选择组合为k的种类数。利用这个 ，然后再去重处理。
详解HERE

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
#define cl(a,b) memset(a,b,sizeof(a))
#define ll long long
#define pb push_back
#define gcd __gcd

const double EPS = 1e-8;
const int maxn = 3e5+1000;
const int inf  = 0x3f3f3f3f;
const double PI = acos(-1.0);
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保证只交换一次
        k=len/2;//i做正常的加1，j做反转的加一，始终保持i，j是反转的
        while(j>=k){
            j-=k;k/=2;
        }
        if(j<k)j+=k;
    }
}
/*
    做FFT len必须是2的幂
    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 num[maxn<<1];
ll sum[maxn<<1];
int a[maxn];
complex x[maxn<<1];

int main(){
    int T;scanf("%d",&T);
    while(T--){
        int n;scanf("%d",&n);
        cl(num,0);
        for(int i=0;i<n;i++){
            scanf("%d",&a[i]);
            num[a[i]]++;
        }
        sort(a,a+n);
        int len1=a[n-1]+1;
        int len=1;
        while(len<len1*2)len<<=1;
        for(int i=0;i<len1;i++)x[i]=complex(num[i],0);
        for(int i=len1;i<len;i++)x[i]=complex(0,0);
        fft(x,len,1);
        for(int i=0;i<len;i++)x[i]=x[i]*x[i];
        fft(x,len,-1);
        for(int i=0;i<len;i++)num[i]=(ll)(x[i].r+0.5);
        len=2*a[n-1];
        for(int i=0;i<n;i++)num[a[i]<<1]--;
        for(int i=0;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-=(n-1);
            cnt-=(ll)(n-i-1)*(n-i-2)/2;
            cnt-=(ll)(n-1-i)*i;
        }
        ll tot=(ll)(n-2)*(n-1)*n/6;
        printf("%.7lf\n",(double)cnt/tot);
    }
    return 0;
}