HDU 4609 3-idiots (FFT)

题意：

给你n (n <= 10^5)个边长ai (ai <= 10^5)，问随机取出三条边可以构成三角形的概率是多少。

解题思路：

首先要学会FFT，然后就很好做了。cnt[i]表示边长为i的有多少条，如果要算出所有两条边长的和分别是多少，普通算法是O(n^2)，利用FFT卷积可以O(nlogn)得到所有的，然后去重一下O(n)处理总数即可。具体见代码

/* **********************************************
Author      : JayYe
Created Time: 2013-9-26 20:32:45
File Name   : JayYe.cpp
*********************************************** */

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <algorithm>
using namespace std;
typedef long long ll;
const double pi = acos(-1.0);
const int maxn = 100000 + 5;
const double eps = 1e-6;

struct Complex {
    double a, b;
    Complex() {}
    Complex(double a, double b) : a(a), b(b) {}

    Complex operator + (const Complex& t) const {
        return Complex(a + t.a, b + t.b);
    }

    Complex operator - (const Complex& t) const {
        return Complex(a - t.a, b - t.b);
    }
    Complex operator * (const Complex& t) const {
        return Complex(a*t.a - b*t.b, a*t.b + b*t.a);
    }
};

void brc(Complex *a, int n) {
    int i, j, k;
    for(i = 1, j = n>>1;i < n - 1; i++) {
        if(i < j)   swap(a[i], a[j]);
        k = n>>1;
        while(j >= k) {
            j -= k;
            k >>= 1;
        }
        if(j < k)
            j += k;
    }
}

void FFT(Complex *a, int n, int on) {
    int h, i, j, k, p;
    double r;
    Complex u, t;
    brc(a, n);
    for(h = 2;h <= n; h <<= 1) {
        r = on*2.0*pi / h;
        Complex wn(cos(r), sin(r));
        p = h>>1;
        for(j = 0;j < n;j += h) {
            Complex w(1, 0);
            for(k = j;k < j + p; k++) {
                u = a[k];
                t = w*a[k + p];
                a[k] = u + t;
                a[k+p] = u - t;
                w = w*wn;
            }
        }
    }
    if(on == -1) {
        for(i = 0;i < n; i++)
            a[i].a = a[i].a / n + eps;
    }
}

Complex x1[maxn<<2];
int sa[maxn];
ll cnt[maxn<<1];

void solve() {
    int n, mx = 0;
    memset(cnt, 0, sizeof(cnt));
    scanf("%d", &n);
    for(int i = 0;i < n; i++) {
        scanf("%d", &sa[i]);
        cnt[sa[i]] ++;
        mx = max(mx, sa[i]);
    }
    int N = 1, tmpn = 2*mx+1;
    while(N < tmpn) N <<= 1;
    for(int i = 0;i < N; i++)
        x1[i].a = x1[i].b = 0;
    for(int i = 0;i <= mx; i++)
        x1[i].a = cnt[i];
    FFT(x1, N, 1);
    for(int i = 0;i < N; i++)
        x1[i] = x1[i]*x1[i];
    FFT(x1, N, -1);
    mx <<= 1;
    for(int i = 0;i <= mx; i++)
        cnt[i] = (ll) (x1[i].a + eps);
    // FFT后得到的是cnt[i]表示任取两条边和为i的有多少种
    for(int i = 0;i < n; i++) {
        cnt[sa[i]*2]--;　       // 减去自己边与自己边的和
    }
    // 由于是取出两条边，所以要除二
    for(int i = 0;i <= mx; i++)
        cnt[i] /= 2;
    for(int i = 1;i <= mx; i++)
        cnt[i] += cnt[i-1];     // 算出前缀和
    sort(sa, sa + n);
    ll ans = 0;
    for(int i = 0;i < n; i++) {
        ans += cnt[mx] - cnt[sa[i]];
        ans -= (ll)(n - i - 1)*(n - i - 2)/2 + n-1 + (ll)(n - i - 1)*i;
    }
    printf("%.7lf
", (double)ans / ((ll)n*(n-1)*(n-2)/6));
}

int main() {
    int t;
    scanf("%d", &t);
    while(t--) {
        solve();
    }
    return 0;
}