[HDU 4609] 3 Idiots FFT

好像有一段时间没写博客了。。
明天就期中考试。。希望能rp大爆发。。

传送门·Teleport!

题目大意

给N条长度不一的木棍,求随机选取三根木棍能组成三角形的概率。
精度到1e-7。

解题报告

附赠样例:

Input:
5
5
1 2 2 5 7
9
1 7 2 6 9 2 5 1 7
10
1 1 5 5 1 3 6 2 6 9
8
1 2 3 3 2 4 5 1
11
1 5 2 7 4 3 7 7 8 1 9
Output:
0.1000000
0.2619048
0.2416667
0.3214286
0.3818182

这题从排列和组合两种方向进行思考都可以,Hineven这里从组合出发。

合法方案 = 总方案 – 不合法方案
运用快速傅里叶变换来统计不合法方案。
把木棍长度看成指数,把该长度的木棍数量看成系数,构造多项式后复制两份相乘。
发现指数在多项式乘法中以加法运算的,但是多项式乘法中系数是乘积之和。
那么乘出来结果的指数为i的项的系数刚好是两个木棍的长度和为i的选取方案数。
然后发现会有两个人同时选了一根相同木棒的方案被统计了,明显当i为偶数时才会出现这种情况,减去即可。

然后就是结果了,注意开64位整数进行中间运算,FFT的话用只用double就可以。
最后说一句。。这题我数组开小。。写到其他数组里了,然后显示Wa,弄了好久整个人都不好了。。
注意FFT数组要开4倍而不是2倍。

接下来是代码, 迭代实现++, 速度快6倍!

代码

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <climits>
#include <vector>
#include <cmath>

typedef long long LL;
const int MAXN = 401000;
const double pi = acos(-1);

struct Complex {
    double r, img;
    Complex () {}
    Complex (double a, double b): r(a), img(b) {}
    inline Complex operator + (const Complex & b) {return Complex(r+b.r, img+b.img);}
    inline Complex operator - (const Complex & b) {return Complex(r-b.r, img-b.img);}
    inline Complex operator * (const Complex & b) {return Complex(r*b.r-img*b.img, r*b.img+img*b.r);}
};

inline int bitReverse (int num, int lvl) {
    int ret = 0, bit = 1;
    while(bit < lvl) {
        ret <<= 1;
        ret += (bool)(num & bit);
        bit <<= 1;
    }
    return ret;
}

Complex temp[MAXN<<1];
inline void fft (Complex * A, int N, double t) {
    for(int i = 0; i<N; i++) temp[bitReverse(i, N)] = A[i];
    for(int i = 2; i<=N; i<<=1) {
        Complex wn(cos(2.0*pi/(double)i), t*sin(2.0*pi/(double)i)), t1, t2;
        for(int j = 0; j<N; j+=i) {
            Complex w(1.0, 0.0);
            for(int k = 0; k<(i>>1); k++, w=w*wn) {
                t1 = temp[j+k], t2 = w*temp[j+(i>>1)+k];
                temp[j+k] = t1+t2, temp[j+(i>>1)+k] = t1-t2;
            }
        }
    }
    for(int i = 0; i<N; i++) A[i] = temp[i];
}

int stick[MAXN];
Complex X[MAXN<<1], Y[MAXN<<1];
int main () {
    int cas;
    scanf("%d", &cas);
    for(int tc = 1; tc <= cas; tc ++) {
        int N, M = 0, t;
        scanf("%d", &N);
        memset(stick, 0, sizeof stick);
        for(int i = 1; i<=N; i++)
            scanf("%d", &t), M = std :: max(M, t), stick[t] ++;
        memset(X, 0, sizeof X), memset(Y, 0, sizeof Y);
        for(int i = 0; i<=M; i++)
            X[i].r = Y[i].r = stick[i];
        int lvl = 1;
        while(lvl <= ((M+1)<<1)) lvl <<= 1;
        fft(X, lvl, 1.0), fft(Y, lvl, 1.0);
        for(int i = 0; i<lvl; i++) X[i] = X[i]*Y[i];
        fft(X, lvl, -1.0);
        LL deno = (LL)N*(LL)(N-1)*(LL)(N-2);
        LL ans = 0, cur = 0;
        for(int i = 0; i<lvl; i++) {
            cur += floor(X[i].r/(double)lvl+0.4);
            if(!(i&1)) cur -= stick[i>>1];
            ans += cur*stick[i]*3LL;
        }
        printf("%.7lf\n", (double)(deno-ans)/(double)deno);
    }
}