解题报告（二）E、（BZOJ3513） [MUTC2013] idiots（生成函数 + FFT + 组合计数）

繁凡出品的全新系列：解题报告系列 —— 超高质量算法题单，配套我写的超高质量题解和代码，题目难度不一定按照题号排序，我会在每道题后面加上题目难度指数（$1\sim 5$），以模板题难度 $1$ 为基准。

---

这样大家在学习算法的时候就可以执行这样的流程：

阅读我的【学习笔记】 / 【算法全家桶】学习算法 $\Rightarrow$ 阅读我的相应算法的【解题报告】获得高质量题单 $\Rightarrow$ 根据我的一句话题解的提示尝试自己解决问题 $\Rightarrow$ 点开我的详细题解链接学习巩固（好耶）
%

题单链接：解题报告（二）多项式问题（多项式乘法及其各种运算）（ICPC / CCPC / NOIP / CF / AT / NC / P）超高质量题解

解题报告（二）E、（BZOJ3513） [MUTC2013] idiots（生成函数 + FFT + 组合计数）

Weblink

https://darkbzoj.tk/problem/3513

Problem

给定 $n$ 个长度分别为 $a_i$ 的木棒，问随机选择 $3$ 个木棒能够拼成三角形的概率。

input

第一行T( $T\le 100$ )，表示数据组数。
接下来若干行描述 T 组数据，每组数据第一行是 n ，接下来一行有 n 个数表示 $a_i$。

$3\le N\le 10^5,1\le a_i\le 10^5$

output

T 行，每行一个整数，四舍五入保留7位小数。

Solution

首先我们知道作为一个三角形，两边之和大于第三边。

答案要求的概率很明显就是能组成的三角形的方案数除以总方案数。

如果我们直接去统计一共有多少个符合要求的方案数的话无从下手（好吧，正着做好像也没什么区别），考虑经典正难则反。

我们首先计算总方案数，即从 $n$ 个木棒中选择三个，即

$tot=C_n^3=n\times (n-1)\times (n-2)/3/2$。

那么我们来尝试统计一下不符合要求的三个木棒的方案数。

不符合要求也就意味着是两个木棒 $a,b$ ，以及一个木棒 $c$ ，$c\ge a+b$,这样就不能组成一个三角形。

因为是统计方案数，且数据较大（$10^5$）不能暴力，只能 $O(nlogn)$ 来做，所以经典的 生成函数 + FFT。

所以我们考虑组合计数的思路：

即利用乘法原理，我们只需要求出左半部分 $A$ 的个数以及相对应的右半部分 $B$ 的个数，相乘并累加即可。

也就是：

我们用 $t[i]$ 统计所有大于等于 $i$ 的木棒的个数。

$g[i]$ 表示所有两个木棒的和为 $i$ 的木棒的个数，那么所有不符合要求的方案数也就是两个木棒 $a,b$ 以及一个 $c\ge a+b$ 的方案数显然就是 $illegal=\sum _{i=0}^{maxx}g[i]\times t[i]$ 。

最终的概率显然就是 $\frac{tot-illegal}{tot}$

那么也就意味着我们只需要计算 $g[i]$ 就可以得到答案。

我们考虑如何求得 $g[i]$。

$g[i]$ 的实际意义就是 $a+b=i$ 的方案数，可以理解为选择两个物品 $a,b$，权值和为 $i$ 的方案数，显然我们可以用生成函数求解。

我们设 $f[i]$ 表示长度为 $i$ 的木棒数

显然 $g[i]=\sum f[j]\times f[i-j]$。

一个标准的卷积形式可以直接用 FFT 求解。

但是很明显我们将两个 $f$ 卷起来会有重复，因为我们只有一个 $f$ 序列，自己凑，但是生成函数在卷的时候是当成两个相同的 $f$ 序列卷起来的，例如 $1+9=10$，两个 $f$ ，所以会有两个 $1$ 和两个 $9$ ，但实际上只有一个 $1$ 和 1个 $9$ ，方案数应该是 $[1,9],[9,1]$ 两种，但是生成函数 $f$ 卷的时候会得到 $[1,9],[9,1],[1,9],[9,1]$ 四种方案数，所以最后的方案数要除以 $2$。并且注意对于序列中仅有一个 $5$ ，卷的时候会出现 $5+5=10$ 这种本来是没有的方案数，如果我们除以二的话也不一定能将其消掉（奇数个除以二下取整会消掉，但万一有偶数个这样的呢）所以我们可以在输入的时候，每次输入 $x$ ，就将 $g[2\times x]-1$ ，这样就可以将本身不合法的自己加自己（$5+5=10$） 方案数减去。正因如此，我们需要先将 $f$ 与 $g$ 加起来消掉自己加自己之后再除以二消掉冗余的双倍方案。

AC Code

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <map>
#include <queue>
using namespace std;
typedef long long ll;
typedef int itn;
#define int long long
typedef pair<int, int>PII;
const int N = 5e5 + 7, mod = 1e6;
const double PI = acos(-1.0);

int n, m, k, L, limit = 1;
int g[N], t[N], RR[N];

struct Complex
{
    double x, y;
    Complex (double x = 0, double y = 0) : x(x), y(y) { }
}f[N];

Complex operator + (Complex a, Complex b) {return Complex (a.x + b.x, a.y + b.y);}
Complex operator - (Complex a, Complex b) {return Complex (a.x - b.x, a.y - b.y);}
Complex operator * (Complex a, Complex b) {return Complex (a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}

void FFT(Complex *A, double type)
{
    for(int i = 0; i < limit; ++ i) {
        if(i < RR[i])
            swap(A[i], A[RR[i]]);
    }
    for(int mid = 1; mid < limit; mid <<= 1) {
        Complex wn(cos(PI / mid), type * sin(PI / mid));

        for(int len = mid << 1, pos = 0; pos < limit; pos += len) {
            Complex w(1, 0);

            for(int k = 0; k < mid; ++ k, w = w * wn) {
                Complex x = A[pos + k];
                Complex y = w * A[pos + mid + k];
                A[pos + k] = x + y;
                A[pos + mid + k] = x - y;
            }
        }
    }
    if(type == -1) {
        for(int i = 0; i <= limit; ++ i)
            A[i].x /= limit;
    }
}

void init()
{
    memset(g, 0, sizeof g);
    memset(f, 0, sizeof f);
    memset(t, 0, sizeof t);

}

void solve()
{
    init();
    scanf("%lld", &n);
    int maxx = 0;
    ll tot = (n * (n - 1) * (n - 2)) / 6;
    for(int i = 1; i <= n; ++ i) {
        ll x;
        scanf("%lld", &x);
        g[x << 1] -- ;
        t[x] ++ ;
        f[x].x ++ ;
        maxx = max(maxx, x);
    }
    for(int i = maxx - 1; i >= 1; -- i)
        t[i] += t[i + 1];//大于等于 i 的木棒数量
    maxx = maxx - 1 << 1;
    L = 0, limit = 1;
    while(limit <= maxx) limit <<= 1, L ++ ;
    for(int i = 0; i <= limit; ++ i) {
        RR[i] = (RR[i >> 1] >> 1) | ((i & 1) << (L - 1));
    }
    FFT(f, 1);
    for(int i = 0; i < limit; ++ i) {
        f[i] = f[i] * f[i];
    }
    FFT(f, -1);
    for(int i = 0; i < limit; ++ i) {
        g[i] += (int)(f[i].x + 0.5);
    }
    for(int i = 0; i < limit; ++ i)
        g[i] >>= 1;
    ll illegal = 0;
    for(int i = 0; i < limit; ++ i) {
        illegal += g[i] * t[i];
    }

    double res = (double)(tot - illegal) / tot;
    printf("%.7f\n", res);
    return ;
}

signed main()
{
    itn t;
    scanf("%lld", &t);
    while(t -- ) {
        solve();
    }
    return 0;
}