HDU 4609 3-idiots (FFT)

题目传送门哦

---

写在前面

做题感想

一道FFT果题写了一个中午，调了一个下午，题目还叫3-idiots，我感觉自己就是一idiot。。。

---

关于快速傅立叶变换(FFT)

蒟蒻我的理解：
FFT就是通过点值、插值等步骤将O(n2)的多项式乘法(或者说卷积)优化到O(nlogn)的算法。
FFT就是快速实现离散傅立叶变换(DFT)、逆离散傅立叶变换(IDFT)的过程。

该方法主要分为几个步骤：
①先将两个系数表达的多项式转换为点值表达(DFT)，这个过程是O(nlogn)的。
②将两个点值表达的多项式乘起来，O(n)。
③逆点值过程，求插值得到系数表达的卷积结果。O(nlogn)。

其实在①之前还有使次数界变为2的幂(当然先使次数界相同)，加倍次数界等预处理。

深入学习FFT，推荐看《算法导论》第30章(第三版)，或者看神犇罗指导的博客(Here)，里面还有讲到CZT(项数不变的DFT)、NTT(快速数论变换)等更高级的东西。

关于FFT的具体实现方法，这里就不多讲了。有直接递归实现分治的写法，亦有迭代的写法。这里推荐迭代写法，因为既可以避免爆栈(好像并不会。。)好像还跑得快一点。递归的写法更简明，迭代更难写一点。但总之理解后熟记模版，做几题慢慢就会熟练了的(自我安慰×)。。

---

题目解法

题目大意就是给你N根木棍，求任选三根成三角形的概率。(3<=N<=100000)

我一开始想的是FFT+线段树。FFT搞每两个的和，然后线段树查找一段区间累计答案。FFT求和就是将读入转成多项式，加法就是将乘法理解为指数的加法，出现次数根据乘法原理(算排列)就是系数相乘。于是长度就是指数，次数就是系数。然后看网上别人的做法中并没有线段树。原来最后直接做一前缀和就行了，要线段树干啥？真是数据结构毁我青春，越学越idiot。

弄完FFT后枚举再除以组合数，我一开始想去统计合法情况，但这样其实很麻烦。FFT大神tututu告诉我只用算非法的再用总的减就行了。因为非法只用违反三条不等式(任意两边…)中的一条，且不会算重(举几个例子就懂了)，然而算合法的可能算重。例如1+4>5，1+5>4就算重了。(我真是CAI)。还有，由于每根棍只出现一次，以及排列改成组合，最后要去掉一些答案。

然后我就赶紧码，码完编译通过输出nan。真强！发现FFT码错了，w跟wm搞反了，然后又调试了很久，修改了N个地方，终于1A了。

果题码半天，真被自己的智商给感动了。

---

程序

#include <cstdio>#include <cstdlib>#include <cmath>#include <cstring>#include <iostream>#include <algorithm>#define MAXN 100010using namespace std;const double PI = acos(-1.0);typedef long long LL;int cnt, n, nG, P[MAXN];int br[MAXN];LL sum[MAXN<<2];struct Complex{    double real, image;    Complex() {}    Complex(double _real, double _image){      real = _real;      image = _image;    }    friend Complex operator + (Complex A, Complex B){return Complex(A.real + B.real, A.image + B.image);}    friend Complex operator - (Complex A, Complex B){return Complex(A.real - B.real, A.image - B.image);}    friend Complex operator * (Complex A, Complex B){return Complex(A.real * B.real - A.image * B.image, A.real * B.image + A.image * B.real);}}a[MAXN<<2], b[MAXN<<2], c[MAXN<<2];void Get_P(){    for(int i = 1, t = 1; i < MAXN; i++){      while(t < i)  t <<= 1;      P[i] = t;    }}void Reverse(Complex *A){    for(int i = 0; i < n-1; i++){      int j = 0;      for(int k = 1, tmp = i; k < n; k <<= 1, tmp >>= 1)        j = ((j << 1) | (tmp & 1));      if(j > i)  swap(A[i], A[j]);    }}void FFT(Complex *A, int n, int DFT){    Reverse(A);    for(int s = 1; (1<<s) <= n; s++){      int m = (1 << s);      Complex wm = Complex(cos(DFT*PI*2/m), sin(DFT*PI*2/m));      for(int k = 0; k < n; k += m){        Complex w = Complex(1, 0);        for(int j = 0; j < (m>>1); j++){          Complex u = A[k+j], t = w * A[k+j+(m>>1)];          A[k+j] = u + t;          A[k+j+(m>>1)] = u - t;          w = w * wm;         }      }    }    if(DFT == -1)  for(int i = 0; i < n; i++)  A[i].real /= n, A[i].image /= n;}int main(){    freopen("hdu4609.in", "r", stdin);    freopen("hdu4609.out", "w", stdout);    Get_P();    scanf("%d", &nG);    while(nG --){      scanf("%d", &cnt);      memset(a, 0, sizeof(a));      n = 0;      for(int i = 1; i <= cnt; i++){        scanf("%d", &br[i]);        a[br[i]].real ++;        n = max(n, br[i]);      }      n = (P[n+1] << 1);      memcpy(b, a, sizeof(b));      FFT(a, n, 1);      FFT(b, n, 1);      for(int i = 0; i < n; i++)  c[i] = a[i] * b[i];      FFT(c, n, -1);      for(int i = 0; i < n; i++)  sum[i] = (LL)(c[i].real+0.5);      for(int i = 1; i <= cnt; i++)  sum[br[i]*2] --;      for(int i = 0; i < n; i++)  sum[i] /= 2;      for(int i = 1; i < n; i++)  sum[i] += sum[i-1];      LL ans, tmp = 1;      for(int i = cnt; i > cnt-3; i--)  tmp = tmp * i;  tmp /= 6ll;      ans = tmp;      for(int i = 1; i <= cnt; i++)  ans -= sum[br[i]];      printf("%.7lf\n", (double)ans / tmp);    }    return 0;}

---