HDU 4609 3-idiots

传送门： http://acm.hdu.edu.cn/showproblem.php?pid=4609

题意

给n条边，问任取三条边是否能构成三角形的概率。

思路

$首先把题意转化一下，即找下列式子的个数：$
$$\frac{{\sum }_{i=1}^n{\sum }_{j=1}^n{\sum }_{k=1}^n{a}_i+a_j>a_k}{C_n^3}(i\ne j\ne k且a_i和a_j是较小的两条边)$$

$看到三角形，首先想到就是把所有的a_i+a_j的个数存在一个数组里，然后再逐一判断每一个a_k是否符合。$

$但是，O(n^2)的复杂度是会T的，但是O(nlogn)的复杂度不会T，那么就需要用到FFT的知识了。$

$FFT是两个多项式相乘，但是上面是要让我们两个相加，没关系，我们让a_i为x的指数，a_i出现的次数为该x的系数，用一个cnt数组保存，即：$

$$A(x)=cnt[a[0]]x^{a[0]}+cnt[a[1]]x^{a[1]}+...+cnt[a[n-1]]x^{a[n-1]}$$

$所以我们让A(x)\times A(x)，就可以得到任意两条边之和，并且得到该和的个数，并存在num数组里。当然，这里会出现很多问题。$

$1、i\ne j，所以要让a_i+a_i的个数减1,即num[a[i]+a[i]]--;$
$2、先取a_1后取a_3和先取a_3后取a_1是一样的，所以我们需要让每一个和都要除2，即num[i]/=2;$

$然后我们用一个sum数组存下num数组的前缀和，然后遍历每条边，取两条边之和比该条边大的个数，即ans+=sum[a[n-1]\ast 2]-sum[a[i]]，前面的a[n-1]\ast 2是前缀和最大的地方。$

$但是，还是有很多问题！！！在刚开始的时候说过，两条边之和必须是两条较小的边，所以还要处理下面几种情况：$

$1、两条边都大于a[i]，即ans-=(n-i-1)\ast (n-i-1)/2;$
$2、一条边大于a[i]，一条边小于a[i]，即ans-=(n-i-1)\ast i;$
$3、两条边中有一条边等于a[i],即ans-=(n-1);$

$最后的ans就是满足条件的所有个数，然后计算C_n^3，最后算概率，保留7位小数。$

$但是，博主tsb了！！！明知道n有1e5，然后计算C_n^3会爆int，所以我还用一个longlong的变量保存下来，即t=n\ast (n-1)\ast (n-2)/6，是不是有一点问题，那就是要在n之前加一个(longlong)强制转换一下，我不过在这个地方错了10多发吧。。。$

Code（826MS）

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;
typedef long double ld;
typedef pair<int, int> pdd;

#define INF 0x3f3f3f3f
#define lowbit(x) x & (-x)
#define mem(a, b) memset(a , b , sizeof(a))
#define FOR(i, x, n) for(int i = x;i <= n; i++)

// const ll mod = 998244353;
// const ll mod = 1e9 + 7;
// const double eps = 1e-6;
const double PI = acos(-1);

const int N = 4e5 + 10;

struct Complex {
    double r, i;
    Complex() {r == 0; i = 0;};
    Complex(double real, double imag) : r(real), i(imag) {};
}A[N], B[N];

Complex operator + (Complex a, Complex b) {
    return Complex(a.r + b.r, a.i + b.i);
}

Complex operator - (Complex a, Complex b) {
    return Complex(a.r - b.r, a.i - b.i);
}

Complex operator * (Complex a, Complex b) {
    return Complex(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r);
}

int rev[N], len, lim = 1;

void FFT(Complex *a, int opt) {
    for(int i = 0;i < lim; i++) {
        if(i < rev[i])
            swap(a[i], a[rev[i]]);
    }
    for(int dep = 1;dep <= log2(lim); dep++) {
        int m = 1 << dep;
        Complex wn = Complex(cos(2.0 * PI / m), opt * sin(2.0 * PI / m));
        for(int k = 0;k < lim; k += m) {
            Complex w = Complex(1, 0);
            for(int j = 0;j < m / 2; j++) {
                Complex t = w * a[k + j + m / 2];
                Complex u = a[k + j];
                a[k + j] = u + t;
                a[k + j + m / 2] = u - t;
                w = w * wn;
            }
        }
    }
    if(opt == -1) { // 逆变换
        for(int i = 0;i < lim; i++) a[i].r /= lim;
    }
}

int getBit(int n) {
    while(lim < (n << 1)) {
        lim <<= 1;
        len++;
    }
    for(int i = 0;i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
}

ll a[N], cnt[N];
ll sum[N], num[N];

void solve() {
    int T;
    scanf("%d",&T);
    while(T--) {
        mem(rev, 0);
        mem(A, 0);
        mem(num, 0);
        mem(cnt, 0);
        int n;
        len = 0, lim = 1;
        scanf("%d",&n);
        for(int i = 0;i < n; i++) {
            scanf("%lld",&a[i]); // 把a[i]看成x的指数，使乘法变加法
            cnt[a[i]]++; // 把a[i]看成每个x^{a_i}的系数，然后再FFT
        }
        sort(a, a + n);
        getBit(a[n - 1] + 1);
        int Max = a[n - 1] * 2;
        for(int i = 0;i <= a[n - 1]; i++) {
            A[i].r = cnt[i];
        }
        FFT(A, 1);
        for(int i = 0;i < lim; i++) {
            A[i] = A[i] * A[i];
        }
        FFT(A, -1);
        for(int i = 0;i < lim; i++) {
            num[i] = A[i].r = (ll)(A[i].r + 0.5); // FFT之后每个a_i+a_j的系数存在num里
        }
        for(int i = 0;i < n; i++) {
            num[a[i] + a[i]]--; // 把i=j的部分减掉
        }
        sum[0] = 0;
        for(int i = 0;i <= Max; i++) {
            num[i] /= 2; // 重复取，除2即可
            sum[i] = sum[i - 1] + num[i]; // 取前缀和，最后判断的时候方便一点
        }
        ll ans = 0;
        for(int i = 0;i < n; i++) {
            ans += sum[Max] - sum[a[i]]; // 统计大于a[i]
        }
        for(int i = 0;i < n; i++) { // 处理符合两边大于a_k但是不正确的个数(因为判断是否为三角形必须是两条较小的边大于最长边)
            ans -= ll(n - i - 1) * (n - i - 2) / 2; // 两边都大于a[i]
            ans -= ll(n - i - 1) * i; // 一边大于a[i]，一边小于a[i]
            ans -= ll(n - 1); // 两边边中有一条边等于a[i]
        }
        ll t = 1ll * n * (n - 1) * (n - 2) / 6;

        printf("%.7lf\n",(double)ans / t);
    }
}

signed main() {
    ios_base::sync_with_stdio(false);
    //cin.tie(nullptr);
    //cout.tie(nullptr);
#ifdef FZT_ACM_LOCAL
    freopen("in.txt", "r", stdin);
    freopen("out.txt", "w", stdout);
    signed test_index_for_debug = 1;
    char acm_local_for_debug = 0;
    do {
        if (acm_local_for_debug == '$') exit(0);
        if (test_index_for_debug > 20)
            throw runtime_error("Check the stdin!!!");
        auto start_clock_for_debug = clock();
        solve();
        auto end_clock_for_debug = clock();
        cout << "Test " << test_index_for_debug << " successful" << endl;
        cerr << "Test " << test_index_for_debug++ << " Run Time: "
             << double(end_clock_for_debug - start_clock_for_debug) / CLOCKS_PER_SEC << "s" << endl;
        cout << "--------------------------------------------------" << endl;
    } while (cin >> acm_local_for_debug && cin.putback(acm_local_for_debug));
#else
    solve();
#endif
    return 0;
}