非递归的快速傅立叶变换

快速傅里叶变换没什么好说的, 只是需要注意避免递归.

避免递归

上面那个递归版本的FFT实际上常数很大，因为要复制数组。我们观察一下递归树，假设这是一个长度为8的多项式。

+===============================================+
| 000   001   010   011   100   101   110   111 |
+-----------------------------------------------+
| 000   010   100   110 | 001   011   101   111 |
+-----------------------------------------------+
| 000   100 | 010   110 | 001   101 | 011   111 |
+-----------------------------------------------+
| 000 | 100 | 010 | 110 | 001 | 101 | 011 | 111 |
+===============================================+

实际上，从叶子往上做的话就是把每个数的二进制位颠倒过来。

 

struct Complex
{
  double real, imag;
  Complex(const double r = .0, const double i = .0): real(r), imag(i) { }
};
const Complex ZERO(.0, .0);
inline Complex operator+(const Complex& lhs, const Complex& rhs) { return Complex(lhs.real+rhs.real, lhs.imag+rhs.imag); }
inline Complex operator-(const Complex& lhs, const Complex& rhs) { return Complex(lhs.real-rhs.real, lhs.imag-rhs.imag); }
inline Complex operator*(const Complex& lhs, const Complex& rhs) { return Complex(lhs.real*rhs.real-lhs.imag*rhs.imag, lhs.real*rhs.imag+lhs.imag*rhs.real); }
inline int bit_reverse(const int x, const int n)
{
  int res = 0;
  for (int i = 0; i < n; ++i)
    res |= (x>>i&1)<<(n-i-1);
  return res;
}
void fft(Complex y[], const Complex a[], const int n, const int rev)
{
  const int len = 1<<n;
  for (int i = 0; i < len; ++i) y[i] = a[bit_reverse(i, n)];
  for (int d = 1; d <= n; ++d)
  {
    const int m = 1<<d;
    const Complex wn(std::cos(2*PI/m*rev), std::sin(2*PI/m*rev));
    for (int k = 0; k < len; k += m)
    {
      Complex w(1., .0);
      for (int j = 0; j < m/2; ++j)
      {
        const Complex u = y[k+j], t = w*y[k+j+m/2];
        y[k+j] = u+t, y[k+j+m/2] = u-t;
        w = w*wn;
      }
    }
  }
  if (rev == -1)
    for (int i = 0; i < len; ++i) y[i].real /= len, y[i].imag = .0;
}
void convolution(Complex c[], Complex a[], Complex b[], const int la, const int lb)
{
  static Complex tmp1[MAXN], tmp2[MAXN];
  int n = 0, len = 1<<n;
  for (; len < 2*la || len < 2*lb; ++n) len <<= 1;
  std::fill(a+la, a+len, ZERO);
  std::fill(b+lb, b+len, ZERO);
  fft(tmp1, a, n, 1);
  fft(tmp2, b, n, 1);
  for (int i = 0; i < len; ++i) tmp1[i] = tmp1[i]*tmp2[i];
  fft(c, tmp1, n, -1);
}