HDU 4609 FFT

题意

传送门 HDU 4609 3-idiots

题解

$N$ 个数中任取三个，求选择的数可以作为三角形边长的概率。对于三角形的判定条件，即 $a_i+a_j>a_k,a_k\ge a_j,a_i$。计数时还需考虑所选择的数互不相同。

令 $A(x)=\sum x^{a_i}$，简单起见写为 $\sum x$，根据容斥原理，考虑相同项的数量，可以得到 $(\sum x{)}^2=\sum x^2+2\sum xy$。$[x^k]\sum xy$ 即选择两个不同数使之满足求和为 $k$ 的方案数。最终减去不满足 $i<j<k$ 的方案数即可，这样的方案包含了使 $a$ 有序后索引大于等于 $k$ 的方案。

#include <bits/stdc++.h>
using namespace std;
using cp = complex<double>;
constexpr double PI = acos(-1);
vector<int> rev;
vector<cp> omg, omgInv;
struct Poly : vector<cp>
{
    Poly() {}
    Poly(int n) : vector<cp>(n) {}
    Poly(const initializer_list<cp> &list) : vector<cp>(list) {}
    void fft(int n, bool inverse)
    {
        if ((int)rev.size() != n)
        {
            rev.resize(n);
            omg.resize(n);
            omgInv.resize(n);
            for (int i = 0; i < n; ++i)
            {
                rev[i] = rev[i >> 1] >> 1 | (i & 1 ? n >> 1 : 0);
                omg[i] = cp(cos(2 * PI * i / n), sin(2 * PI * i / n));
                omgInv[i] = conj(omg[i]);
            }
        }
        resize(n);
        for (int i = 0; i < n; ++i)
            if (i < rev[i])
                std::swap(at(i), at(rev[i]));

        for (int m = 1; m < n; m <<= 1)
        {
            int m2 = m << 1;
            for (int i = 0; i < n; i += m2)
                for (int j = 0; j < m; ++j)
                {
                    cp &x = at(i + j), &y = at(i + j + m);
                    cp t = (inverse ? omgInv[n / m2 * j] : omg[n / m2 * j]) * y;
                    y = x - t;
                    x += t;
                }
        }
    }
    void dft(int n) { fft(n, 0); };
    void idft(int n)
    {
        fft(n, 1);
        for (int i = 0; i < n; ++i)
            at(i) /= n;
    }
    Poly operator*(const Poly &p) const
    {
        auto a = *this, b = p;
        int k = 1, n = a.size() + b.size() - 1;
        while (k < n)
            k <<= 1;
        a.dft(k), b.dft(k);
        for (int i = 0; i < k; ++i)
            a[i] *= b[i];
        a.idft(k);
        a.resize(n);
        return a;
    }
};
using ll = long long;
constexpr int MAXN = 1E5 + 5;
int T, N, A[MAXN];
ll C[MAXN][4];

int main()
{
    ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    for (int i = 0; i < MAXN; ++i)
        C[i][0] = 1;
    for (int i = 1; i < MAXN; ++i)
        for (int j = 1; j <= 3; ++j)
            C[i][j] = C[i - 1][j - 1] + C[i - 1][j];
    cin >> T;
    while (T--)
    {
        cin >> N;
        int amax = 0;
        for (int i = 0; i < N; ++i)
            cin >> A[i], amax = max(amax, A[i]);
        sort(A, A + N);
        Poly f(amax + 1);
        vector<int> a(amax + 1);
        for (int i = 0; i < N; ++i)
            ++a[A[i]];
        for (int i = 0; i <= amax; ++i)
            f[i].real(a[i]);
        f = f * f;
        int n = f.size();
        vector<ll> b(n);
        for (int i = 0; i < n; ++i)
            b[i] = floor(f[i].real() + 0.5);
        for (int i = 0; i < N; ++i)
            b[A[i] * 2]--;
        for (int i = 0; i < n; ++i)
            b[i] /= 2;
        vector<ll> sum(n);
        for (int i = 1; i < n; ++i)
            sum[i] = sum[i - 1] + b[i];
        ll res = 0;
        for (int i = N - 1; i >= 0; --i)
        {
            ll t = sum[n - 1] - sum[A[i]];
            t -= (ll)(i + N - 1) * (N - i) / 2;
            res += t;
        }
        cout << fixed << setprecision(7) << (double)res / C[N][3] << '\n';
    }
    return 0;
}