【HDU 4609】3-idiots（FFT）

这道题也是FFT的入门应用题之一了，求任意不同的三条边能形成三角形的概率。

首先我们可以先用FFT求出任意两条（可相同）的边能组成的长度的方案数。

那么现在就需要去重了，我们可以想到对于组成长度为奇数的边直接除以2即可（因为与x,y的先后顺序无关）

对于长度为偶数的，还要先减去自己与自己组成的方案数，除以2后再加回去与先后顺序无关的方案数。

算出方案数后，我们可以通过枚举第三条边来确定方案数。

可以想到求出不可行的方案数比较方便，因此直接前缀和搞一下就可以了。

具体见代码

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long LL;
const int N = 100005;
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 &a)const{
        return Complex(r + a.r,i + a.i);
    }
    Complex operator - (const Complex &a)const{
        return Complex(r - a.r,i - a.i);
    }
    Complex operator * (const Complex &a)const{
        return Complex(r * a.r - i * a.i,r * a.i + i * a.r);
    }
    Complex operator * (const double &a)const{
        return Complex(r * a,i * a);
    }
};
int cnt[N << 2],tn,val[N];
LL res[N << 2],sum[N << 2];
Complex x1[N << 2],tmp[N << 2];
int rev(int x){
    int ans = 0;
    for(int i = 0 ; i < tn ; i ++){
        if(x & 1) ans += 1 << tn - 1 - i;
        x >>= 1;
    }
    return ans;
}
void fft(Complex A[],int n,int op){
    for(int i = 0 ; i < n ; i ++) tmp[ rev(i) ] = A[i];
    for(int i = 0 ; i < n ; i ++) A[i] = tmp[i];
    for(int i = 1 ; (1 << i) <= n ; i ++){
        int m = 1 << i;
        Complex wn(cos(op * 2 * pi / m),sin(op * 2 * pi / m));
        for(int k = 0 ; k < n ; k += m){
            Complex w(1,0),u,t;
            for(int j = 0 ; j < m / 2 ; j ++){
                u = A[k + j];
                t = w * A[k + j + m / 2];
                A[k + j] = u + t;
                A[k + j + m / 2] = u - t;
                w = w * wn;
            }
        }
    }
    if(op == -1)
        for(int i = 0 ; i < n ; i ++)
            A[i] = A[i] * (1.0 / n);
}
void debug(int len){
    for(int i = 0 ; i < len ; i ++) printf("%d ",i);printf("\n");
    for(int i = 0 ; i < len ; i ++) printf("%I64d ",res[i]);printf("\n");
    for(int i = 0 ; i < len ; i ++) printf("%I64d ",sum[i]);printf("\n");
}
void solve(){
    int n,t,Max = 0;
    scanf("%d",&n);
    memset(cnt,0,sizeof(cnt));
    for(int i = 1 ; i <= n ; i ++){
        scanf("%d",&t);
        cnt[t] ++;
        Max = max(Max,t);
    }
    tn = ceil(log(Max + 1.0) / log(2.0)) + 1;
    int len = 1 << tn;
    sum[len] = 0;
    for(int i = len - 1 ; i >= 0 ; i --) sum[i] = sum[i + 1] + cnt[i];
    for(int i = 0 ; i < len ; i ++) x1[i] = Complex(cnt[i],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 ++) res[i] = LL(x1[i].r + 0.5);
    for(int i = 0 ; i < len ; i ++){
        if(i % 2 == 0) res[i] -= (LL)cnt[i / 2] * cnt[i / 2];
        res[i] /= 2;
        if(i % 2 == 0) res[i] += (LL)cnt[i / 2] * (cnt[i / 2] - 1) / 2;
    }
    double ans1 = 1.0 * n * (n - 1) * (n - 2) / 6;
    double ans2 = 0;
    for(int i = 0 ; i < len ; i ++) ans2 += 1.0 * res[i] * sum[i];
    printf("%.7f\n",1.0 - ans2 / ans1);
}
int main()
{
    int _;
    scanf("%d",&_);
    while(_--) solve();
    return 0;
}