HDU - 4609 FFT

题意：

给出n个长度，要求任意选择3个长度，能组成三角形的概率。

思路：

这也是一种常见的FFT题型，初看这个可能很难和FFT能解决的多项式联系到一起。

首先既然要组成三角形，一定要保证两边之和大于第三边，那么我们可以枚举一条边为x，然后看不超过x的两边之和一共有多少情况，这个很容易，只要求出来两边之和==0,1,..Max的种类数，然后用一个前缀和保存就行了。

现在关键是怎么对于任意一个x，求出有多少对长度的和等于x，令cnt[x]保存长度为x的边有多少个

这样样例给出四条边 { 1，3，3，4 }就可以表示成一个cnt 向量 ：{ 0，1，0，2，1 }

如果不考虑重复或者选择自身两次，可以发现任意两个和的向量为： { 0，0，1，0，4，2，4，4，1 }

对于所求的结果可以这样理解：

取两个数和为 2 的取法是一种：1+1
              和为 4 的取法有四种：1+3，1+3，3+1，3+1
              和为 5 的取法有两种：1+4，4+1；
              和为 6 的取法有四种： 3+3，3+3，3+3，3+3，3+3
              和为 7 的取法有四种： 3+4，3+4，4+3，4+3
              和为 8 的取法有一种：4+4

这里其实可以这样考虑如果要求和为x的取法是f[x]，则：

如果把cnt向量看作一个多项式的系数，其实求解f(x)的过程就可以看作多项式相乘：

这样，就可以用FFT来加速计算了，但是要注意上面求出的种类数包括了选择自身两次以及(a,b)和(b,a)的重复，这些都要去掉。

代码：

#include <bits/stdc++.h>
using namespace std;
#define pi acos (-1)
const int MAXN = 555555;

struct plex {       // 定义复数类
    double x, y;
    plex (double _x = 0.0, double _y = 0.0) : x (_x), y (_y) {}
    plex operator + (const plex &a) const {
        return plex (x + a.x, y + a.y);
    }
    plex operator - (const plex &a) const {
        return plex (x - a.x, y - a.y);
    }
    plex operator * (const plex &a) const {
        return plex (x * a.x - y * a.y, x * a.y + y * a.x);
    }
};

void change (plex *y, int len) {
    for (int i = 1, j = len / 2; i < len - 1; i++) {
        if (i < j) swap(y[i], y[j]);
        int k = len / 2;
        while (j >= k) {
            j -= k;
            k /= 2;
        }
        if (j < k) j += k;
    }
}

void fft(plex y[], int len, int on) {   // FFT过程，on==1时，将系数表达转换成点值表达，on==-1时，将点值表达转换成系数表达
    change(y, len);
    for(int h = 2; h <= len; h <<= 1) {
        plex wn(cos(-on * 2 * pi / h), sin(-on * 2 * pi / h));
        for(int j = 0; j < len; j += h) {
            plex w(1, 0);
            for(int k = j; k < j + h / 2; k++) {
                plex u = y[k];
                plex 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].x /= len;
    }
}

long long a[MAXN], cnt[MAXN], sum[MAXN];
plex x[MAXN];

int main () {
    int T;
    scanf("%d", &T);
    while (T--) {
        int n;
        long long Max = 0;
        scanf("%d", &n);
        memset(cnt, 0, sizeof(cnt));
        for (int i = 1; i <= n; i++) {
            scanf("%I64d", &a[i]);
            cnt[a[i]]++;
            Max = max(Max, a[i]);
        }
        ++Max;
        int len = 2;
        while (len < Max * 2) len <<= 1;
        for (int i = 0; i <= Max; i++)
            x[i] = plex(cnt[i], 0);
        for (int i = Max + 1; i < len; i++)
            x[i] = plex(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++) {
            cnt[i] = (long long) (x[i].x + 0.5);
        }
        for (int i = 1; i <= n; i++)
            cnt[a[i] + a[i]]--;
        for (int i = 0; i < len; i++)
            cnt[i] /= 2;
        sum[0] = 0;
        for (int i = 1; i < len; i++)
            sum[i] = sum[i - 1] + cnt[i];
        //sort (a + 1, a + 1 + n);
        long long tot = (long long)n * (n - 1) * (n - 2) / 6LL, ans = tot;
        for (int i = 1; i <= n; i++) {
            ans -= sum[a[i]];
        }
        printf("%.7f\n", ans * 1.0 / tot);
    }
    return 0;
}