3-idiots hdu4609 母函数+FFT 组合数学题

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

题意：1e5个数，求取三个数能形成三角形的概率。

题解（这怎么会是fft入门题QAQ）：

　　概率的算法就是三角形取法/总取法。总取法就是C（n，3）.

　　三角形取法如何计算？

part1：构造母函数F（日常套路），每一项的次数为长度，系数为该长度的木棍数量，用FFT算F^2 ，

　　　　得到的多项式就包含了任意取两跟棍子得到的所有长度的方案数：其中次数为两根棍长之和，系数为该长度的方案数，

part2：去重，考虑part1中得到的系数，并非方案数（举个例子，（a+b）^2==a^2+b^2+2ab）

　　　　首先，对于每个平方项，因为棍子不能重复使用，所以对于k^m要减去平方前k^(m/2）的系数，

　　　　其次，对于剩下的每一项，由于多项式乘法每个乘法都做了两遍，所以得到的系数为方案数的两倍。要/2.

　　　　我们将处理过的系数存入数组A。

part3：算方案数：为了避免重复，对于每一根的木棍，计算以它为最长边的三角形方案数：为该长度到最大长度的A的系数和，（前缀和O（1）算出）

　　　　

　　　　减去其中包含它的方案数，为1*（n-1）.

　　　　减去其中它不是最长的方案数，为（n-1-1）*比它长的木棍数-两根都比它长的方案数/2（乘法原理多算了一遍，容斥减去），也可以理解为 一根比它长一根比他短+两根都比他长/2.（其中比它长的木棍数并不需要树状数组，只要排个序，他就是n-i）

end

总之，fft只是一个开始，后面巨麻烦orz

坑：套miskcoo dalao的fft，结果发现别人的是魔改版的，需要深刻地理解fft才会用orz

ac代码：

#include<iostream>
#include<sstream>
#include<fstream>
#include<vector>
#include<list>
#include<deque>
#include<queue>
#include<stack>
#include<map>
#include<set>
#include<bitset>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cctype>
#include<cmath>
#include<ctime>
#include<iomanip>
#include<complex>
using namespace std;
//const double eps(1e-8);
typedef long long lint;

const double PI = acos(-1.0);

const int MaxL = 18, MaxN = 1 << MaxL;
typedef complex<double> complex_t;
complex_t f[MaxN], g[MaxN];
complex_t eps[MaxN], inv_eps[MaxN];
void init_eps(int p)
{
    double pi = acos(-1);
    //double angle = 2.0 * pi / p;
    for (int i = 0; i != p; ++i)
        eps[i] = complex_t(cos(2.0 * pi * i / p), sin(2.0 * pi * i / p));
    for (int i = 0; i != p; ++i)
        inv_eps[i] = conj(eps[i]);
}



void transform(int n, complex_t *x, complex_t *w)
{
    for (int i = 0, j = 0; i != n; ++i)
    {
        if (i > j) std::swap(x[i], x[j]);
        for (int l = n >> 1; (j ^= l) < l; l >>= 1);
    }

    for (int i = 2; i <= n; i <<= 1)
    {
        int m = i >> 1;
        for (int j = 0; j < n; j += i)
        {
            for (int k = 0; k != m; ++k)
            {
                complex_t z = x[j + m + k] * w[n / i * k];
                x[j + m + k] = x[j + k] - z;
                x[j + k] += z;
            }
        }
    }
}


int branch[100010];
int num[200010];
complex_t a[264000];//(1 << 17 = 131072, 1 << 18 = 262144)
lint A[200010];
lint sumA[200010];//表示A[i]的前缀和

int main()
{
    int T, n;
    scanf("%d", &T);
    while (T--)
    {
        scanf("%d", &n);
        int maxBranch = 0;
        for (int i = 0; i < n; i++)
        {
            scanf("%d", branch + i);
            maxBranch = max(maxBranch, branch[i]);
        }
        memset(num, 0, sizeof(num));
        for (int i = 0; i < n; i++)
            num[branch[i]]++;
        for (int i = 0; i <= maxBranch; i++)
            a[i] = num[i];
        int len = 1;
        while (len <= maxBranch) len <<= 1;
        len <<= 1;
        for (int i = maxBranch + 1; i < len; i++)
            a[i] = 0.0;
        init_eps(len);
        transform( len,a ,eps);
        for (int i = 0; i < len; i++)
            a[i] = a[i] * a[i];
        transform( len,a, inv_eps);
        for (int i = 0; i <= 2 * maxBranch; i++)
            A[i] = (lint)(a[i].real() + 0.5)/len;
        for (int i = 0; i <= 2 * maxBranch; i += 2)
            A[i] -= num[i >> 1];
        for (int i = 0; i <= 2 * maxBranch; i++)
            A[i] /= 2;
        //到现在为止A[i]表示的是取两根不同的branch的长度和为i的组合种数
        sumA[0] = 0;
        for (int i = 1; i <= 2 * maxBranch; i++)
            sumA[i] = sumA[i - 1] + A[i];
        lint ans = 0;
        sort(branch, branch + n);
        for (int i = 1; i <= n; i++)//以第i根作为边最长的
        {
            lint tmp = sumA[2 * maxBranch] - sumA[branch[i]];//另外两条边长度和要大于branch[i]
            tmp -= (lint)(n - i)*(n - 2);//比它长
            tmp += (lint)(n - i)*(n - i - 1) / 2;//两条都比他长
            tmp -= n - 1;//另外两条的组合中包括它自己的组合
            ans += tmp;
        }
        double p = ans * 6. / n / (n - 1) / (n - 2);
        printf("%.7f
", p);
    }
    cin >> n;
    return 0;
}
/*
*/
/*
2
4
1 3 3 4
4
2 3 3 4

*/