HDU 4609 3-idiots (FFT)

题目传送门哦

---

写在前面

做题感想

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

---

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

蒟蒻我的理解：
FFT就是通过点值、插值等步骤将$O(n^2)$的多项式乘法(或者说卷积)优化到$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 100010


using 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;
}

---