hdu 4609 3-idiots（FFT）

Posted 2021-11-09 jpphy0

目录

• 问题
• 分析
• 代码

问题

hdu 4609 3-idiots - https://acm.hdu.edu.cn/showproblem.php?pid=4609

分析

• fft 方法计算两边之和
• 两边之和去重
  • 同一条边取两次，即 a + a = 2a
  • 两条边选取的次序交换，即 a + b = b + a，因此除2
• 两边之和大于第3边，按第3边最大枚举第3边去重
  • 将n条边排序
  • 升序枚举各边，前缀和方式计算大于当前边的两边可取数
  • 以下三种情况因重复不符去除
    • 两条边均比当前边长，设当前边是第i条，则去除$(n-i)\ast (n-i-1)/2$
    • 一条边不小于当前边，一条边不大于当前边，则去除$(i-1)\ast (n-i)$
    • 当前边已被使用，则去除$n-1$

代码

#include<bits/stdc++.h>
using namespace std;
#define LL long long
const double DFT = 2.0, IDFT = -2.0, PI = acos(-1);
const int MXN = 1e5+10, MXL = MXN << 2, MXV = MXN;
struct CP{
    double r, i;
    CP(double _r = 0, double _i = 0):r(_r), i(_i){};
    CP operator +(CP x){ return CP(r+x.r, i+x.i); }
    CP operator -(CP x){ return CP(r-x.r, i-x.i); }
    CP operator *(CP x){ return CP(r*x.r - i*x.i, r*x.i + i*x.r); }
    CP operator *=(CP x){ 
        double xr = x.r, xi = x.i, tr = r;
        r = r*xr - i*xi, i = tr*xi + i*xr;
        return *this;
    }
    void print(char pad){ printf("%d%c", (int)(r+0.5), pad); }
}p[MXL];
int n, r[MXL], d[MXN], cnt[MXV];
LL sum[MXV<<1], ans, tot;
void fft(CP p[], int len, double mode){
    for(int i = 0; i < len; ++i) if(i < r[i]) swap(p[i], p[r[i]]);
    for(int i = 2; i <= len; i <<= 1){
        CP wn(cos(mode*PI/i), sin(mode*PI/i));
        for(int j = 0; j < len; j += i){
            CP w(1.0, 0.0);
            for(int k = j; k < j+(i>>1); ++k, w *= wn){
                CP l = p[k], r = w*p[k+(i>>1)];
                p[k] = l + r, p[k+(i>>1)] = l - r;
            }
        }
    }
    if(mode == DFT) return;
    for(int i = 0; i <= len; ++i) p[i].r /= len;
}
int main(){
    int t, lim, bit, mx;
    scanf("%d", &t);
    while(t--){        
        memset(cnt, 0, sizeof cnt), memset(sum, 0, sizeof sum);
        p[0].r = 0, p[0].i = 0, mx = 0, ans = 0;
        scanf("%d", &n);
        for(int i = 1; i <= n; ++i) scanf("%d", d+i), mx = max(mx, d[i]), ++cnt[d[i]];
        for(int i = 1; i < MXV; ++i) p[i].r = cnt[i], p[i].i = 0;
        for(int i = MXV; i < MXL; ++i) p[i].r = 0, p[i].i = 0;
        lim = mx, bit = 2;
        while(lim >>= 1) ++bit;
        lim = 1 << bit;
        for(int i = 0; i < lim; ++i) r[i] = (r[i>>1]>>1)|((i&1)<<(bit-1));
        fft(p, lim, DFT);
        for(int i = 0; i <= lim; ++i) p[i] *= p[i];
        fft(p, lim, IDFT);
        for(int i = 1; i < MXV; ++i) if(cnt[i]) p[i<<1].r -= cnt[i]; // 同一个取两次去重
        for(int i = 0; i < (MXV<<1); ++i) p[i].r /= 2; // 不同的两个先后取去重
        for(int i = 1; i < (MXV<<1); ++i) sum[i] += sum[i-1] + (LL)(p[i].r+0.5);
        sort(d+1, d+n+1);
        for(int i = 1; i <= n; ++i){
            ans += sum[(mx<<1)+1] - sum[d[i]];
            ans -= 1LL*(i-1)*(n-i);
            ans -= n-1;
            ans -= 1LL*(n-i)*(n-i-1)/2;
        }
        tot = 1LL * n * (n - 1) * (n - 2) / 6;
        printf("%.7f\\n", 1.0*ans/tot);
    }
    return 0;
}

• 两边之和小于等于第三边不构成三角形
• 两边之差大于等于第三边不构成三角形
• 但第二种情形总可以被第一种情形包含

#include<bits/stdc++.h>
using namespace std;
typedef complex<double> cd;
typedef long long ll;
const double DFT = 2.0, IDFT = -2.0, PI = acos(-1);
const int MXN = 1e5+10, MXL = MXN << 2, MXV = MXN;
cd p[MXL];
int n, r[MXL], c[MXV];
ll tmp, pre, tot, sum[MXV<<1];
void fft(cd p[], int len, double mode){
    for(int i = 0; i < len; ++i) if(i < r[i]) swap(p[i], p[r[i]]);
    for(int i = 2; i <= len; i <<= 1){
        cd wn(cos(mode*PI/i), sin(mode*PI/i));
        for(int j = 0; j < len; j += i){
            cd w(1.0, 0.0);
            for(int k = j; k < j+(i>>1); ++k, w *= wn){
                cd l = p[k], r = w*p[k+(i>>1)];
                p[k] = l + r, p[k+(i>>1)] = l - r;
            }
        }
    }
    if(mode == DFT) return;
    for(int i = 0; i <= len; ++i) p[i] /= len;
}
int main(){
    int t, lim, bit;    
    scanf("%d", &t);
    while(t--){        
        memset(c, 0, sizeof c), p[0] = 0, lim = 0;
        scanf("%d", &n);
        for(int x, i = 1; i <= n; ++i) scanf("%d", &x), lim = max(lim, x), ++c[x];
        for(int i = 1; i < MXV; ++i) p[i] = c[i];
        for(int i = MXV; i < MXL; ++i) p[i] = 0;
        bit = 2;
        while(lim >>= 1) ++bit;
        lim = 1 << bit;
        for(int i = 0; i < lim; ++i) r[i] = (r[i>>1]>>1)|((i&1)<<(bit-1));
        fft(p, lim, DFT);
        for(int i = 0; i <= lim; ++i) p[i] *= p[i];
        fft(p, lim, IDFT);
        for(int i = 1; i < MXV; ++i) if(c[i]) p[i<<1] -= c[i];
        for(int i = 0; i < (MXV<<1); ++i) sum[i] = (ll)(p[i].real()/2 +0.5);
        tmp = tot = 1LL * n * (n - 1) * (n - 2) / 6, pre = 0;
        for(int i = 0; i < MXV; ++i) pre += sum[i], tmp -= 1LL*pre*c[i]; 
        printf("%.7f\\n", 1.0*tmp/tot);
    }
    return 0;
}