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}^na_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]*2]-sum[a[i]]$，前面的$a[n-1]*2$是前缀和最大的地方。

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

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

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

但是，博主tsb了！！！明知道n有1e5，然后计算$C_n^3$会爆int，所以我还用一个longlong的变量保存下来，即t=n*(n-1)*(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;
}

本文作者：jujimeizuo
本文地址： https://blog.jujimeizuo.cn/2020/10/09/hdu-4609/
本博客所有文章除特别声明外，均采用 CC BY-SA 3.0 协议。转载请注明出处！