bzoj3513

fft

考虑什么情况下不能组成三角形，那么就是a+b<=c,我们用fft求出所有a+b的情况，对于每个c统计就行了

#include<bits/stdc++.h>
using namespace std;
#define pi acos(-1)
const int N = 500005;
int n, m, k;
int A[N];
long long t[N], sum[N];
struct data {
    double a, b;
    data() { a = 0; b = 0; }
    data(double _, double __) : a(_), b(__) {}
    data friend operator + (const data &a, const data &b) { return data(a.a + b.a, a.b + b.b); }
    data friend operator - (const data &a, const data &b) { return data(a.a - b.a, a.b - b.b); }
    data friend operator * (const data &a, const data &b) { return data(a.a * b.a - a.b * b.b, a.a * b.b + a.b * b.a); }
} a[N];
void fft(data *a, int n, int f)
{
    for(int i = 0; i < n; ++i) 
    {
        int t = 0;
        for(int j = 0; j < k; ++j) if(i >> j & 1) t |= 1 << (k - j - 1);
        if(i < t) swap(a[i], a[t]); 
    }
    for(int l = 2; l <= n; l <<= 1) 
    {
        int m = l >> 1;
        data w = data(cos(pi / m), f * sin(pi / m));
        for(int i = 0; i < n; i += l) 
        {
            data t = data(1, 0);
            for(int k = 0; k < m; ++k, t = t * w)
            {
                data x = a[i + k], y = t * a[i + k + m];
                a[i + k] = x + y;
                a[i + k + m] = x - y;
            }
        }
    }
    if(f == -1) 
    {
        for(int i = 0; i < n; ++i) t[i] = (long long)(a[i].a / n + 0.1);
    }
}
int main()
{
    int T;
    scanf("%d", &T);
    while(T--)
    {
        scanf("%d", &n);
        int lim = 0;
        for(int i = 1; i <= n; ++i) 
        {
            scanf("%d", &A[i]);
            lim = max(lim, 2 * A[i]);
            ++sum[A[i]];
            a[A[i]].a += 1.0;
        }
        double ans = (double)n * (n - 1) * (n - 2) / 6.0, tmp = ans;
        m = 1;
        k = 0;
        while(m <= lim) m <<= 1, ++k;
        fft(a, m, 1);
        for(int i = 0; i < m; ++i) a[i] = a[i] * a[i];
        fft(a, m, -1);      
        for(int i = 1; i <= n; ++i) --t[A[i] * 2];
        for(int i = 1; i <= lim; ++i) t[i] = (t[i] >> 1) + t[i - 1];
        for(int i = 1; i <= lim; ++i) tmp -= (double)sum[i] * t[i]; 
        for(int i = 0; i < m; ++i) a[i] = data(0, 0), sum[i] = t[i] = 0;
        printf("%.7f
", tmp / ans);
    }
    return 0;
}