HDU4609 3-idiots

多校那时候我还年轻，没能理解FFT，只知道它能够进行信号在时域和频域上的转换。

现在学习了FFT在多项式乘法上的应用，感觉对FFT的理解更深了一层。

要搞懂FFT，首先要理解卷积。

先抄一下wikipedia

函数 f 与 g 的卷积记作，它是其中一个函数翻转并平移后与另一个函数的乘积的积分，是一个对平移量的函数。

积分区间取决于 f 的定义域与 g 的定义域。

对于定义在离散域的函数，卷积定义为

理解了卷积之后，我们可以把FFT先简单的理解为求卷积的O(nlogn)算法。

然后进入正题，本题中我们已知每一种长度的棒子各有多少个，换个思路想想，我们就知道了一个多项式的各项系数，把这个多项式与其自身做快速傅里叶变换，就可以得到新多项式的各项系数，对于新系数ki'，有k[i]' = sigma( k[j] + k[i-j] )，而这里的系数正是由两根棒子加起来的长度的种数，假设用cnt[i]表示和的长度为i的棒子组合有多少种。

求出数组cnt[]之后，这个题目已经做了一半了，下面是计数问题。

(1) 首先去掉cnt里重复的部分，一根棒子不能与自己组合，所以有cnt[a[i] * 2] --

(2) 任意两个棒子组合的顺序我们不需要考虑，而实际上他们算了两次，所以有cnt[i] /=2

(3) 枚举所有棒子，假设这根棒子a[i]是组成三角形其中最长的那根，我们首先求出另外两根的长度和比它的长的总种数，即sigma(cnt[j])   (a[i]<j<=crest*2，crest是最长的棒子长度)

(4) 另外两根中，这根被枚举的棒子不能再出现了，所以要减去n-1

(5)  另外两根中，可能有一根长度大于a[i]，另一根小于a[i]，所以要减去(n-i) * (i-1)

(6) 另外两根中，可能两根长度都大于a[i]，再减去(n-i) * (n-i-1) / 2

计算下所有可能的种数 tot = n * (n-1) * (n-2) / 6

把统计出的答案除以tot就行了。

注意：cnt 和 tot 都有可能超过int的范围，所以应该用long long 

贴代码：

#include <cstdio>#include <cstring>#include <cctype>#include <cstdlib>#include <ctime>#include <climits>#include <cmath>#include <iostream>#include <string>#include <vector>#include <set>#include <map>#include <list>#include <queue>#include <stack>#include <deque>#include <algorithm>using namespace std;typedef long long ll;const double PI = acos(-1.0);const int maxn = 100010;struct Complex{    double r, i;    Complex(double r=0, double i=0) : r(r), i(i) {}    Complex operator + (const Complex &o) {return Complex(r+o.r, i+o.i);}    Complex operator - (const Complex &o) {return Complex(r-o.r, i-o.i);}    Complex operator * (const Complex &o) {return Complex(r*o.r-i*o.i, r*o.i+o.r*i);}};void brc(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;    }}void fft(Complex *y, int len, int on){    brc(y, len);    for (int i = 2; i <= len; i <<= 1)    {        Complex wn(cos(on*2.0*PI/i), sin(on*2.0*PI/i));        for (int j = 0; j < len; j += i)        {            Complex w(1.0, 0.0);            for (int k = j; k < j+i/2; k++)            {                Complex u = y[k];                Complex t = w*y[k+i/2];                y[k] = u + t;                y[k+i/2] = u - t;                w = w * wn;            }        }    }    if (on == -1)        for (int i = 0; i < len; i++)            y[i].r /= len;}int T,n;int a[maxn];ll cnt[maxn*4], sum[maxn*4], ans, tot;Complex x1[maxn*4];int main(){    scanf("%d", &T);    while (T--)    {        scanf("%d", &n);        memset(cnt, 0, sizeof(cnt));        int len = 1, crest = 0;        for (int i = 1; i <= n; i++)        {            scanf("%d", &a[i]);            if (a[i] > crest) crest = a[i];            cnt[a[i]]++;        }        crest++;        while (len < crest*2) len <<= 1;        for (int i = 0; i < len; i++) x1[i] = Complex(cnt[i], 0.0);        fft(x1, len, 1);        for (int i = 0; i < len; i++) x1[i] = x1[i] * x1[i];        fft(x1, len, -1);        for (int i = 0; i < len; i++) cnt[i] = (ll)(x1[i].r + 0.5);        len = crest * 2;        for (int i = 1; i <= n; i++) cnt[a[i]*2]--;        for (int i = 1; i <= len; i++) cnt[i] /= 2;        sum[0] = 0;        for (int i = 1; i <= len; i++) sum[i] = sum[i-1] + cnt[i];        ans = 0;        for (int i = 1; i <= n; i++)        {            ans += sum[len] - sum[a[i]];            ans -= n - 1;            ans -= (ll)(i - 1) * (n - i);            ans -= (ll)(n - i) * (n - i - 1) / 2LL;        }        tot = (ll)n * (n - 1) * (n - 2) / 6LL;        printf("%.7lf\n", (double)ans / tot);    }return 0;}