HDU 4609 FFT模板

http://acm.hdu.edu.cn/showproblem.php?pid=4609

题意：给你n个数，问任意取三边能够，构成三角形的概率为多少。

思路：使用FFT对所有长度的个数进行卷积（\(C = \{x + y| x \in A, y \in B \} \)），得到所有两边和情况数，再考虑去掉重复的情况。找边并计数的时候再去掉不可能的情况。具体操作看bin神的博客   

另FFT还可以用来进行多项式和高精度乘法，又难懂又难用的东西=x=

/** @Date    : 2016-12-04-16.31
  * @Author  : Lweleth ([e-mail])
  * @Link    : https://github.com/
  * @Version :
  */

#include<bits/stdc++.h>
#define LL long long
#define PII pair
#define MP(x, y) make_pair((x),(y))
#define fi first
#define se second
#define PB(x) push_back((x))
#define MMG(x) memset((x), -1,sizeof(x))
#define MMF(x) memset((x),0,sizeof(x))
#define MMI(x) memset((x), INF, sizeof(x))
using namespace std;

const int INF = 0x3f3f3f3f;
const int N = 1e5+20;
const double eps = 1e-8;
const double PI = acos(-1.0);

struct Complex
{
    double a, i;//实部，虚部
    Complex(double aa = 0, double ii = 0)
    {
        a = aa, i = ii;
    }
    Complex operator +(const Complex &y)
    {
        return Complex(a + y.a, i + y.i);
    }
    Complex operator -(const Complex &y)
    {
        return Complex(a - y.a, i - y.i);
    }
    Complex operator *(const  Complex &y)
    {
        return Complex(a * y.a - i * y.i, a * y.i + i * y.a);
    }
};

void change(Complex y[], int len)//len 必须为2的幂；位置i和二进制反转的位置互换
{
    int i, j, k;
    for(i = 1, j = len / 2; i < len - 1; i++)
    {
        if(i < j)
            swap(y[i], y[j]);
        //交换互为小标反转的元素，i<j交换一次
        //i正常++ j左反转类型的+1 始终保持二者反转
        k = len / 2;
        while(j >= k)
        {
            j -= k;
            k /= 2;
        }
        if(j < k)
            j += k;
    }
}

void fft(Complex y[], int len, int on)//on==1时DFT -1时IDFT len必须为2的幂
{
    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].a /= len;
}

Complex x1[4 * N];
int a[N];
LL num[4 * N];
LL sum[4 * N];

int main()
{
    int T;
    cin >> T;
    while(T--)
    {
        int n;
        scanf("%d", &n);
        MMF(num);
        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 < 2 * len1)
            len <<= 1;
        for(int i = 0; i < len1; i++)
            x1[i] = Complex(num[i], 0);
        for(int i = len1; i < len; i++)
            x1[i] = Complex(0 , 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++)
            num[i] = (LL)(x1[i].a + 0.5);
        len = 2 * a[n - 1];

        /////
        for(int i = 0; i < n; i++)//卷积中同样的数字选取两次
            num[a[i] + a[i]]--;
        for(int i = 1; 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 -= (LL)(n - i - 1) * i;
            cnt -= (LL)(n - 1) * 1;
            cnt -= (LL)(n - i - 1) * (n - i - 2) / 2;
        }
        printf("%.7lf\n", (double)cnt * 6.0 / n / (n - 1.0) / (n - 2.0));
    }
    return 0;
}