HDU-4609 3-idiots

3-idiots

大意

给你 \(n\) 根木根，任选 3 根，能组成三角形的概率是多少？

题目

King OMeGa catched three men who had been streaking in the street. Looking as idiots though, the three men insisted that it was a kind of performance art, and begged the king to free them. Out of hatred to the real idiots, the king wanted to check if they were lying. The three men were sent to the king's forest, and each of them was asked to pick a branch one after another. If the three branches they bring back can form a triangle, their math ability would save them. Otherwise, they would be sent into jail.
However, the three men were exactly idiots, and what they would do is only to pick the branches randomly. Certainly, they couldn't pick the same branch - but the one with the same length as another is available. Given the lengths of all branches in the forest, determine the probability that they would be saved.

Input:

An integer T(T≤100) will exist in the first line of input, indicating the number of test cases.
Each test case begins with the number of branches N(3≤N≤10^5).
The following line contains N integers a_i (1≤a_i≤10 5), which denotes the length of each branch, respectively.

Output:

Output the probability that their branches can form a triangle, in accuracy of 7 decimal places.

Sample Input:

2
4
1 3 3 4
4
2 3 3 4

Sample Output:

0.5000000
1.0000000

分析

组成三角形的充要条件任意两边和大于第三边。

不妨设三边为 \(a_i, a_j, a_k\) ，且 \(a_i \le a_j \le a_k\) 那么我们只需要知道组成 \(S = a_i + a_j\) 的方法数前缀和，就可以 \(O(n)\) 地扫过 \(a_k\) 得到总的方案数了。

那么现在就转化成如何计算两边和 \(S\) 的方案数的个数了。

简单粗暴地 \(n^2\) 计算肯定会超时的，我们可以先用生成函数相乘来计算方案数，然后用 FFT 来加速两个两个多项式相乘。总的时间复杂度为 \(O(n\log n)\) 。

设边长为 \(n\) 的棍子数为 \(a_n\) ，那么由生成函数有：

\[f(x) = a_0 + a_1 x^1 + a_2 x^2 + \cdots + a_n x^n \]

两边和的生成函数为：

\[g(x) = f(x) \cdot f(x) = b_0 + b_1 x^1 + b_2x^2 + \cdots +b_{2n}x^{2n} \]

其中 \(b_i\) 的含义为长度和为 \(i\) 的方案数。套 FFT 的板子就可以计算出 \(b_i\) 的各个值了。

样例 1 中的 \(f(x) = x + 2x^3 + x^4\) , \(g(x) = x^2+4x^4+2x^5+4x^6+4x^7+x^8\) 。

那 \(g(x)\) 中的 \(x^2\) 含义就是组成长度和为 \(2\) 的方案个数为 \(1\) 。

由生成函数的意义我们可以知道，这里其实算重了一部分。这里首先对于每一根棍子，它都贡献了一个第一条边选自己，第二条边也选自己的方案。所以我们要减去。也就是：

for (int i = 0; i < n; i++) num[arr[i] * 2]--;

然后，这里选的每根棍子都应该与顺序无关。 \(f(x) \cdot f(x)\) 表示选取的时候有顺序了，我们要把它消去顺序，也就是除以 2 。

for (int i = 0; i < len; i++) num[i] /= 2;

这样算出来的的就是真正的两边和的方案数了。先算一波前缀和。

再然后，对于每根棍子 \(a_k\) 对于答案的贡献就是组成比 \(a_k\) 大的两边和的方案数。也就是：

ans += sum[len - 1] - sum[arr[i]];

但这样选取的两边和可能会造成三种情况：

1. 选了一根比 \(a_k\) 小的，一根比 \(a_k\) 大的
2. 选了 \(a_k\) 自己，然后取其他的
3. 选了两根比 \(a_k\) 大的

那么我们就要相应地减去对应的方案数：

ans -= 1ll * (n - 1 - i) * i;               //情况 1 
ans -= 1ll * (n - 1);                       //情况 2 
ans -= 1ll * (n - 1 - i) * (n - 2 - i) / 2; //情况 3 

最后得到的就是能选取三角形的方案数。除以总的方案数就能得到最后的答案。

参考代码

#include <bits/stdc++.h>
#define pii pair<int, int>
#define pdd pair<double, double>
#define pdi pair<pdd, int>
#define pil pair<int, long long>
#define pll pair<long long, long long>
#define lowbit(x) ((x) & (-x))
#define mem(i, a) memset(i, a, sizeof(i))
#define sqr(x) ((x) * (x))
#define ls(k) (k << 1)
#define rs(k) (k << 1 | 1)
typedef long long ll;
const double eps = 1e-8;
const int inf = 0x3f3f3f3f;
const int maxn = 5e5 + 7;
using namespace std;
const double pi = acos(-1);
struct Complex {
    double x, y;
    Complex(double _x = 0.0, double _y = 0.0) {
        x = _x;
        y = _y;
    }
    Complex operator-(const Complex &b) const {
        return Complex(x - b.x, y - b.y);
    }
    Complex operator+(const Complex &b) const {
        return Complex(x + b.x, y + b.y);
    }
    Complex operator*(const Complex &b) const {
        return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
    }
};
void change(Complex y[], int len) {
    for (int i = 1, j = len / 2, k; i < len - 1; i++) {
        if (i < j) swap(y[i], y[j]);
        k = len / 2;
        while (j >= k) {
            j -= k;
            k /= 2;
        }
        if (j < k) j += k;
    }
}
void fft(Complex y[], int len, int on) {
    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].x /= len;
}
ll f(ll n) {
    bool f2 = 0, f3 = 0;
    ll arr[3], ans = 1;
    for (int i = 0; i < 3; i++) {
        arr[i] = n - i;
        if (!(arr[i] % 2) && !f2) arr[i] /= 2, f2 = 1;
        if (!(arr[i] % 3) && !f3) arr[i] /= 3, f3 = 1;
        ans *= arr[i];
    }
    return ans;
}
ll sum[maxn], num[maxn];
ll arr[maxn];
Complex a[maxn], b[maxn];
int main(void) {
#ifdef ljxtt
    freopen("data.in", "r", stdin);
// freopen("data.out", "w", stdout);
#endif

    int T;
    scanf("%d", &T);
    while (T--) {
        int n;
        scanf("%d", &n);
        for (int i = 0; i < n; i++) {
            scanf("%lld", &arr[i]);
        }
        sort(arr, arr + n);
        int m = arr[n - 1], len = 1;
        while (len <= 2 * m) len <<= 1;
        for (int i = 0; i < len; i++) a[i] = b[i] = Complex(0, 0);
        for (int i = 0; i < n; i++) a[arr[i]].x++, b[arr[i]].x++;
        fft(a, len, 1);
        fft(b, len, 1);
        for (int i = 0; i < len; i++) a[i] = a[i] * b[i];
        fft(a, len, -1);
        mem(num, 0);
        for (int i = 0; i < len; i++) num[i] = ll(a[i].x + 0.5);
        for (int i = 0; i < n; i++) num[arr[i] * 2]--;
        for (int i = 0; i < len; i++) num[i] /= 2;

        sum[0] = num[0];
        for (int i = 1; i < len; i++) {
            sum[i] = sum[i - 1] + num[i];
        }
        ll ans = 0;
        for (int i = 0; i < n; i++) {
            ans += sum[len - 1] - sum[arr[i]];
            ans -= 1ll * (n - 1 - i) * i;
            ans -= 1ll * (n - 1);
            ans -= 1ll * (n - 1 - i) * (n - 2 - i) / 2;
        }
        printf("%.7lf\n", (1.0 * ans) / f(n));
    }
    return 0;
}