多校第一场——hdu4609——傅里叶变换

题意：随机从N个数中取3个数构成三角形的概率是多少。

注意到每个数的范围是1到10^5，可以借鉴生成函数的思想，把任意两个数的和看作多项式的系数用多项式乘法乘出来。可是一般的乘法的时间很作急，就要用到傅里叶变换了。

#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long lng;

const int maxn = (1 << 18) + 100;
const double pi = acos(-1.0);

struct fft
{
    int len, id;
    int sp[maxn];
    pair<double, double> p[3][maxn], tmp[maxn];
    
    void init(int l, int * s)
    {
        this->len = l;
        for(int i = 0; i < len; ++i)
            sp[i] = s[i];
    }
    
    void fill(int m, int d, int c)
    {
        if(m == len) p[c][d] = tmp[id++];
        else
        {
            fill(m << 1, d, c);
            fill(m << 1, d + m, c);
        }
    }
    
    void calc(double oper, int c)
    {
        for(int d = 0; (1 << d) < len; ++d)
        {
            int m = (1 << d);
            double x = pi / m * oper;
            double i0 = sin(x), r0 = cos(x);
            for(int k = 0; k < len; k += (m << 1))
            {
                double i = 0, r = 1;
                for(int j = 0; j < m; ++j)
                {
                    double tr = r * p[c][k + j + m].first - i * p[c][k + j + m].second;
                    double ti = r * p[c][k + j + m].second + i * p[c][k + j + m].first;
                    p[c][k + j + m].first = p[c][k + j].first - tr;
                    p[c][k + j + m].second = p[c][k + j].second - ti;
                    p[c][k + j].first += tr;
                    p[c][k + j].second += ti;
                    double tmpi = i;
                    i = i * r0 + r * i0; r = r * r0 - tmpi * i0;
                }
            }
        }
    }
    
    void multi()
    {
        for(int i = 0; i < len; ++i)
        {
            tmp[i].first = sp[i];
            tmp[i].second = 0.0;
        }
        id = 0; fill(1, 0, 0); calc(1.0, 0);
        id = 0; fill(1, 0, 1); calc(1.0, 1);
        for(int i = 0; i < len; ++i)
        {
            tmp[i].first = p[0][i].first * p[1][i].first - p[0][i].second * p[1][i].second;
            tmp[i].second = p[0][i].first * p[1][i].second + p[0][i].second * p[1][i].first;
        }
        id = 0; fill(1, 0, 2); calc(-1.0, 2);
        for(int i = 0; i < len; ++i) p[2][i].first /= len;
    }
    
}f;

int num[maxn], sum[maxn], n;
lng ans[maxn], len;

void prework()
{
    memset(sum, 0, sizeof(sum));
    scanf("%d", &n);
    len = 1;
    for(int i = 1; i <= n; ++i)
    {
        scanf("%d", num + i);
        sum[num[i]]++;
        while(len <= 2 * num[i]) len <<= 1;
    }
    f.init(len, sum);
    f.multi();
    for(int i = 0; i < len; ++i) ans[i] = f.p[2][i].first + 0.1;
    for(int i = 1; i <= n; ++i) ans[2 * num[i]]--;
    for(int i = 1; i < len; ++i) { ans[i] /= 2; ans[i] += ans[i - 1]; }
    for(int i = 1; i < len; ++i) sum[i] += sum[i - 1];
}

void solve()
{
    lng cnt = 0, tot = ((lng)n * (n - 1) * (n - 2)) / 6;
    for(int i = 1; i <= n; ++i)
    {
        lng lr = sum[num[i] - 1]; lng rr = n - lr - 1;
        lng tmp = ans[len - 1] - ans[num[i]];
        tmp -= n - 1; tmp -= (lr * rr); tmp -= (rr * (rr - 1)) / 2;
        cnt += tmp;
    }
    for(int i = 1; i < len; ++i)
    {
        lng tmp = sum[i] - sum[i - 1];
        if(tmp >= 2)
        {
            cnt += sum[i - 1] * ((tmp * (tmp - 1)) / 2);
            cnt += (tmp * (tmp - 1) * (tmp - 2)) / 6; 
        }
    }
    printf("%.7f\n", cnt * 1.0 / tot);
}

int main()
{
    int t; cin >> t;
    while(t--)
    {
        prework();
        solve();
    }
    return 0;
}