hdu 4609 3-idiots FFT 计数

最新推荐文章于 2019-01-24 22:41:05 发布

sunshine47

最新推荐文章于 2019-01-24 22:41:05 发布

阅读量344

点赞数

分类专栏： FFT-NTT 计数

本文链接：https://blog.csdn.net/sunshine47/article/details/48317163

版权

FFT-NTT 同时被 2 个专栏收录

1 篇文章 0 订阅

订阅专栏

计数

1 篇文章 0 订阅

订阅专栏

Source

hdu 4609

题意

N根棍子，任取三根，问能组成三角形的概率。

分析

考虑Ａ+B > C，所以枚举最长边C，看有多少A+B，就可以FFT了。
后面处理的部分真是想哭=。=我把长度相同的都放起来了，然后不好统计，枚举当前长度 i 作为最长，需要考虑有一根、两根、三根长度为 i 的。

如果不放到一起，直接for的话，会比较好想。。
网上的粗略看了看，大概有这么些方法
容易想到的：
http://www.cnblogs.com/kuangbin/archive/2013/07/24/3210565.html
实际上这个做法不排序也可以：
http://www.cnblogs.com/jklongint/p/4693536.html

求补集，从反面考虑，枚举较小的两边A+B，去掉C >= A+B：
http://blog.csdn.net/jackyguo1992/article/details/12622139
求补集，也是从反面考虑，枚举最长边C，去掉A+B<=C：
http://blog.csdn.net/zz_1215/article/details/38341381

代码

#include<cstdio>
#include<complex>
#include<cmath>
#include<algorithm>
#include<cstring>
using namespace std;

const int N = 1<<18;
const double PI = acos(-1.0);

struct Complex{
    double x, y;
    Complex(double xx = 0, double yy = 0){ x = xx, y = yy; }
    friend Complex operator + (const Complex &a, const Complex &b){
        return Complex(a.x+b.x, a.y+b.y);
    }
    friend Complex operator - (const Complex &a, const Complex &b){
        return Complex(a.x-b.x, a.y-b.y);
    }
    friend Complex operator * (const Complex &a, const Complex &b){
        return Complex(a.x*b.x-a.y*b.y, a.y*b.x+b.y*a.x);
    }
};
typedef Complex CD;
//typedef complex<double> CD;
//注释部分采用stl的complex，会慢些

//oper==1 DFT oper==-1 IDFT
void FFT(CD* a, int n, int oper){        // n must be 2^exp
    for (int i = 1, j = 0; i < n; i++){
        for (int s = n; j^=s>>=1, ~j&s;);
        if (i < j) swap(a[i], a[j]);
    }
    for (int m = 1; m < n; m <<= 1){
        int m2 = m<<1;
        double p = PI/m*oper;
        CD w = CD(cos(p), sin(p));
        for (int i = 0; i < n; i += m2){
            CD unit = 1;
            for (int j = 0; j < m; j++){
                CD &x = a[i+j+m], &y = a[i+j], t = unit*x;
                x = y-t;
                y = y+t;
                //unit *= w;
                unit = unit * w;
            }
        }
    }
    //if (oper == -1) for (int i = 0; i < n; i++) a[i] /= n;
    if (oper == -1) for (int i = 0; i < n; i++) a[i].x /= n;
}

void Square(int* A, int la, long long* C, int& lc){
    int n = 1; while(n < la*2) n <<= 1;
    static CD a[N], B[N];
    for (int i = 0; i < n; i++) a[i] = i<la? A[i]: 0;
    FFT(a, n, 1);
    for (int i = 0; i < n; i++) a[i] = a[i] * a[i];
    FFT(a, n, -1);
    lc = la*2-1;
    for (int i = 0; i < lc; i++) C[i] = (long long)(a[i].x + 0.5);
}

const int maxn = 1e5+10;
int T, n;
int a[maxn], b[maxn];
long long c[maxn<<1];
int main()
{
    scanf("%d", &T);
    while(T--){
        scanf("%d", &n);
        int mx = 0, mm;
        for (int x, i = 0; i < n; i++){
            scanf("%d", &b[i]);
            mx = max(mx, b[i]);
        }
        memset(a, 0, sizeof(int)*(mx+1));
        for (int i = 0; i < n; i++) a[b[i]] ++;
        Square(a, mx+1, c, mm);
        for (int i = 1; i <= mx; i++) c[i+i] -= a[i];
        for (int i = mx+mx; i ; i--) c[i] >>= 1;
        long long cnt = 0, sum = 0, num = 0;
        for (int i = mx+mx; i > mx; i--) sum += c[i];
        for (int i = mx; i > 0; i--){
            if (!a[i]){ sum += c[i]; continue;}
            cnt += (sum-(num+a[i])*n+num+a[i]+(num+a[i])*(num+a[i]-1)/2)*a[i] + (-num+n-a[i])*(a[i]-1)*a[i]/2 + (long long)a[i]*(a[i]-1)*(a[i]-2)/6;
            num += a[i]; sum += c[i];
        }
        printf("%.7f\n", 6.0*cnt/n/(n-1)/(n-2));
    }
    return 0;
}