HDU - 4609 ,3-idiots,FFT

HDU - 4609 ,3-idiots

https://vjudge.net/problem/HDU-4609/origin
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.
题意：有n个木棍，让你从中随机选取3个问能够组成三角形的概率
思路：首先能够组成三角形表示两根长度小的和大于长度最大的，现先求方案数，枚举每根木棍作为最大长度u[i]，则需要求出其他任意两根长度小于等于u[i]的木棍长度和大于u[i]的方案数，设大于u[i]的木棍有x个，小于等于u[i]的木棍有y个
则枚举u[i]时的贡献 = 任意两根不同木棍和大于u[i]的方案数 - $C_x^2$- x*y - (n-1)
即分别排除选两大，一小一大，选u[i]以及另外一个 这三种情况（这三种情况的选的两木棍和都会大于u[i]）
最后即求任意两根不同木棍和大于u[i]的方案数
很明显所有木棍的长度都小于$10^5$，直接记录每种长度的个数，自卷积一下然后排除掉同一木棍选两次的情况，弄个前缀和即可

#include<bits/stdc++.h>
#define ll long long
#define ri register int
#define MAXN 100005
using namespace std;
const double pi=acos(-1.0);
struct cpx
{
    double r, i;
    inline void operator +=(const cpx &b){ r+=b.r, i+=b.i;}
    inline cpx operator +(const cpx &b)const{ return (cpx){r+b.r, i+b.i};}
    inline cpx operator -(const cpx &b)const{ return (cpx){r-b.r, i-b.i};}
    inline cpx operator *(const cpx &b)const{ return (cpx){r*b.r-i*b.i, r*b.i+i*b.r};}
    inline cpx operator *(const double b)const{ return (cpx){r*b, i*b};}
    inline cpx operator ~()const{return (cpx){r, -i};}
}a[MAXN<<2], w[MAXN<<2];

inline void DFT_(cpx *f, int n)
{
    for(ri i=0, j=0; i<n; ++i)
    {
        if(i>j) swap(f[i], f[j]);
        for(ri k=n>>1; (j^=k)<k; k>>=1);
    }
    for(ri i=1; i<n; i<<=1) for(ri j=0; j<n; j+=i<<1)
            for(ri k=j; k<j+i; ++k)
            {
                cpx t=w[i+k-j]*f[k+i];
                f[k+i]=f[k]-t, f[k]+=t;
            }
}
inline void DFT(cpx *f, int n)
{
    if(n==1) return;
    n>>=1;
    static cpx a[MAXN<<1];
    for(ri i=0; i<n; ++i) a[i]=(cpx){f[i<<1].r, f[i<<1|1].r};
    DFT_(a, n);
    for(ri i=0; i<n; ++i)
    {
        cpx q=~a[(n-i)&(n-1)], x=(a[i]+q)*0.5, y=(a[i]-q)*(cpx){0, -0.5}, t=y*w[n+i];
        f[i]=x+t, f[n+i]=x-t;
    }
}
inline void IDFT(cpx *f, int n)
{
    if(n==1) return;
    reverse(f+1, f+n), n>>=1;
    static cpx a[MAXN<<1];
    for(ri i=0; i<n; ++i)
        a[i]=(f[i]+f[i+n])*0.5 + (f[i]-f[i+n])*(cpx){0, 0.5}*w[n+i];
    DFT_(a, n);
    double k=1.0/n;
    for(ri i=0; i<n; ++i) f[i<<1]=(cpx){a[i].r*k, 0}, f[i<<1|1]=(cpx){a[i].i*k, 0};
}

int n,x[MAXN],len1,b[MAXN],u[MAXN];
ll ans[MAXN << 2],sum[MAXN];
int main()
{
    int t;
    scanf("%d",&t);
    while(t--)//输入两个多项式的长度及模数,注意这里的长度比实际长度小1
    {
        scanf("%d",&n);
        len1 = 0;
        memset(x,0,sizeof(x));
        for(int i = 1;i <= n;++i)
        {
            scanf("%d",&u[i]);
            ++x[u[i]];
            len1 = max(len1,u[i]);
        }
        ++len1;
        sort(u+1,u+n+1);
        int len = 1;
        while(len <(len1<<1)) len<<=1;
        for(int i = 0;i < len1;i++)
            a[i] = (cpx){(double)x[i],0};
        for(int i = len1;i < len;i++)
            a[i] = (cpx){0,0};
        for(ri i=1; i<len; i<<=1)
        {
            w[i]=(cpx){1, 0};
            for(ri j=1; j<i; ++j)
                w[i+j]=((j&31)==1?(cpx){cos(pi*j/i), sin(pi*j/i)}:w[i+j-1]*w[i+1]);
        }

        DFT(a,len);
        for(int i = 0;i < len;++i)
            a[i] = a[i]*a[i];
        IDFT(a,len);

        ll res = 0;
        for(int i = 0;i < len;++i)
            ans[i] = (ll)(0.5+a[i].r);
        for(int i = 1;i <= n;++i)
            --ans[u[i]*2];
        for(int i = 1;i < len;++i)
            ans[i] = ans[i-1] + ans[i];
        for(int i = 1;i <= n;++i)
        {
            res += (ans[len-1]-ans[u[i]])/2 - 1ll*(n - i)*(i-1)-1ll*(n-i)*(n-i-1)/2-n+1;
        }
        ll res1 = 1ll*n*(n-1)*(n-2)/6;
        printf("%.7f\n",1.0*res/res1);
    }
    return 0;
}