HDU 4609 3-idiots(FFT)

来源:互联网 发布:手机淘宝运费怎么设置 编辑:程序博客网 时间:2024/06/06 10:39

转载请注明出处,谢谢http://blog.csdn.net/ACM_cxlove?viewmode=contents    by---cxlove

题意 :给出n条边,问选出三条边能组成三角形的概率 

http://acm.hdu.edu.cn/showproblem.php?pid=4609

第一次搞FFT,理论还不是非常清楚,首先要了解卷积。

我只是来存代码的,具体的可以看kuangbin巨巨的解释

http://www.cnblogs.com/kuangbin/archive/2013/07/24/3210565.html 

num[i]表示长度为i的边有几条,求一次num与num的卷积之后,num[i]表示两条边和为i的有多少对。

然后 需要去重一下,最后就可以 O(n)统计了,去重的地方需要注意,blog里有讲很详细


#include <iostream>#include <cstdio>#include <cstring>#include <cmath>#include <algorithm>using namespace std;//FFT copy from kuangbinconst double pi = acos (-1.0);// Complex  z = a + b * i  struct Complex {    double a, b;    Complex(double _a=0.0,double _b=0.0):a(_a),b(_b){}    Complex operator + (const Complex &c) const {        return Complex(a + c.a , b + c.b);    }    Complex operator - (const Complex &c) const {        return Complex(a - c.a , b - c.b);    }    Complex operator * (const Complex &c) const {        return Complex(a * c.a - b * c.b , a * c.b + b * c.a);    }};//len = 2 ^ kvoid change (Complex y[] , int len) {    for (int i = 1 , j = len / 2 ; i < len -1 ; i ++) {        if (i < j) swap(y[i] , y[j]);        int k = len / 2;        while (j >= k) {            j -= k;            k /= 2;        }        if(j < k) j += k;    } }// FFT // len = 2 ^ k// on = 1  DFT    on = -1 IDFTvoid FFT (Complex y[], int len , int on) {    change (y , len);    for (int h = 2 ; h <= len ; h <<= 1) {        Complex wn(cos (-on * 2 * pi / h), sin (-on * 2 * pi / h));        for (int j = 0 ; j < len ; j += h) {            Complex w(1 , 0);            for (int k = j ; k < j + h / 2 ; k ++) {                Complex u = y[k];                Complex t = w * y [k + h / 2];                y[k] = u + t;                y[k + h / 2] = u - t;                w = w * wn;            }        }    }    if (on == -1) {        for (int i = 0 ; i < len ; i ++) {            y[i].a /= len;        }    }}const int N = 100005;typedef long long LL;int n , a[N];LL sum[N << 2] , num[N << 2];Complex x1[N << 2];int main () {    #ifndef ONLINE_JUDGE        freopen("input.txt" , "r" , stdin);    #endif    int t;    scanf ("%d", &t);    while (t --) {        memset (num , 0 , sizeof(num));        scanf ("%d", &n);        for (int i = 0 ; i < n ; i ++) {            scanf ("%d", &a[i]);            num[a[i]] ++;        }        sort (a , a + n);        int len = a[n - 1] + 1;        int l = 1;        while (l < len * 2) l <<= 1;        for (int i = 0 ; i < len ; i ++) {            x1[i] = Complex (num[i] , 0);        }        for (int i = len ; i < l ; i ++) {            x1[i] = Complex (0 , 0);        }        FFT(x1 , l , 1);        for (int i = 0 ; i < l ; i ++) {            x1[i] = x1[i] * x1[i];        }        FFT(x1 , l , -1);        for (int i = 0 ; i < l ; i ++) {            num[i] = (LL)(x1[i].a + 0.5);        }        l = 2 * a[n - 1];        for (int i = 0 ; i < n ; i ++) {            num[a[i] << 1] --;        }        for (int i = 1 ; i <= l ; i ++) {            num[i] /= 2;        }        sum[0] = 0;        for (int i = 1 ; i <= l ; i ++) {            sum[i] = sum[i - 1] + num[i];        }        double ans = 0;        for (int i = 0 ; i < n ; i ++) {            ans += sum[l] - sum[a[i]];            ans -= n - 1;            ans -= (double)i * (n - 1 - i);            ans -= (double)(n - i - 1) * (n - i - 2) / 2;        }        printf ("%.7f\n", ans * 6.0 / n / (n - 1.0) / (n - 2.0));    }    return 0;}