快速傅里叶算法

所以对应到傅里叶级数中，我们就可以知道每个 $a_n$ 、 $b_n$ 实质上就是将 $f(x)$ 用正弦 / 余弦函数分解后，频率为 $\frac{n}{2\pi}$ 的正弦 / 余弦函数的系数。至于为什么我们要把这些系数用 $c_n$ 来表示？~~其实就是纯属为了省纸。~~ 两者在算数意义上是一致的，也仅仅是算数意义上一致，写成复数形式好看而已，不必深究。

为了得到连续的频率，我们只需要把求和换成积分，弃用 $c_n$ （因为已经不是离散了），把它换成一个连续的函数 $F(t)$ 即可：

\begin{align*} f(x)&=\int_{-\infty}^\infty F(t)e^{itx}\mathrm{d}t\\ F(t)&=\frac{1}{2\pi}\int_{-\pi}^\pi f(x)e^{-itx}\mathrm{d}x \end{align*}

走到这一步，我们得到的式子看起来已经和傅里叶变换差不多了，但还是有一点点不一样，比如系数 $\frac{1}{2\pi}$ 的位置，但是这个很好理解，因为 $f(x)$ 、 $F(t)$ 都出现在了对方的积分式里，由积分的性质可知这个系数是可以转移给对方的。

那么只剩下把 $F(t)$ 的积分上下限换成无穷了。这个过程实际上就是对整个函数进行趋于无穷的缩放变换，把函数作用域提升到整个实数域。

所以我们最终得到了如下式子：

\begin{align*} F(t)&=\int_{-\infty}^\infty f(x)e^{-itx}\mathrm{d}x\\ f(x)&=\frac{1}{2\pi}\int_{-\infty}^\infty F(t)e^{itx}\mathrm{d}t \end{align*}

这下成功我们从傅里叶级数开始，得到了傅里叶变换。我们把傅里叶变换抽取有限点计算，就得到了离散傅里叶变换。

离散傅里叶变换

我们复习一下离散傅里叶变换的公式：

\begin{align*} X(k)=\mathrm{DFT}[x(n)]&=\sum_{n=0}^{N-1}x(n)e^{-i\frac{2\pi n}{N}k}\\ x(n)=\mathrm{IDFT}[X(k)]&=\frac{1}{N}\sum_{n=0}^{N-1}X(k)e^{i\frac{2\pi n}{N}k} \end{align*}

其中 $X(k)$ 为 $x(n)$ 的象函数， $x(n)$ 为 $X(k)$ 的象原函数

只要你认识公式，我相信你可以直接实现出一份 $O(n^2)$ 的 DFT 算法，只需要对于每个 $X(k)$ 跑一次循环即可。

但是在竞赛中，这还是太慢了，所以我们需要用到一些小技巧来升级 DFT。

n次单位根

在数学上有这么一类特殊的方程： $x^n=1$ ，由于 n 次方程有 n 个复数根，易得：

x=e^{i\frac{2\pi k}{n}},k=(0,1,2,...,n-1)

形如 $e^{i\frac{2\pi k}{n}}$ 的东西就叫 n次单位根，我们用 $\varepsilon_k$ 来表示，这些根都位于复平面的单位圆上，把单位圆均分成 n 份。

n 次单位根有一些特殊的性质，而这些性质可以用来加速 DFT：

性质 1

\varepsilon_k=\varepsilon_{k\%n}

这个性质可以把式子转换成三角函数，由周期性证得。

性质 2

\varepsilon_k=(\varepsilon_1)^k

这个性质可以把式子转换成自然底数幂证得。

性质 3

\varepsilon_{k+\frac{n}{2}}=-\varepsilon_k

推导：

\begin{align*} \varepsilon_{k+\frac{n}{2}}&=e^{i\frac{2\pi(k+\frac{n}{2})}{n}}\\ &=e^{i\frac{2\pi k+\pi n}{n}}\\ &=e^{i\frac{2\pi k}{n}}\cdot e^{i\pi}\\ &=-e^{i\frac{2\pi k}{n}}\\ &=-\varepsilon_k \end{align*}

n 次本原单位根

这个东西和 n 次单位根比起来好像就是多了两个字，但是两者的意思却完全不一样。

如果 n 个 n 次单位根可以由一个 n 次单位根的 0 到 (n-1) 次幂求得，那么我们称这个 n 次单位根为 n 次本原单位根，我们记为 $\omega_n$ 。

如所有的 4 次单位根为： $1$ 、 $i$ 、 $-1$ 、 $-i$ 。其中 $i$ 和 $-i$ 通过求 0 到 3 次幂可以导出所有的 4 次单位根，而 $1$ 和 $-1$ 却不行，所以 4 次本原单位根为 $i$ 和 $-i$ 。

对于某个 n 来说，它的 n 次本原单位根可能不止一个，但是对于所有的 n 来说，有两个本原单位根都是一样的，那就是 $\varepsilon_1$ 和 $\varepsilon_{-1}$ ，即 $e^{i\frac{2\pi}{n}}$ 和 $e^{-i\frac{2\pi}{n}}$ 。

原因也很简单，由 n 次单位根性质 2，当然可以导出所有的 n 次单位根，只不过 $\varepsilon_1$ 和 $\varepsilon_{-1}$ 只是顺着数和倒着数的区别。

在这里我们需要用到一个 n 次本原根的性质：

\omega_N^{dk}=\omega_{\frac{N}{d}}^k

推导：

\begin{align*} \omega_N^{dk}&=e^{i\frac{2\pi}{N}dk}\\ &=e^{i\frac{2\pi}{\frac{N}{d}}k}\\ &=\omega_{\frac{N}{d}}^k \end{align*}

由 DFT 到 FFT

有了上面的前置知识，我们开始对 DFT 进行加速。

我们令 $\omega_N=\varepsilon_{-1}=e^{-i\frac{2\pi}{N}}$ ，所以 DFT 公式可以写作：

X(k)=\sum_{n=0}^{N-1}x(n)\omega_N^{nk}

将 $X(k)$ 展开，得：

X(k)=x(0)\omega_N^{0k}+x(1)\omega_N^{1k}+x(2)\omega_N^{2k}+...+x(N-1)\omega_N^{(N-1)k}

我们得到了一个 (N-1) 次多项式，当然这里的变量是 $\omega_N^k$ ，所以我们令：

A(\omega_N^k)=x(0)\omega_N^{0k}+x(1)\omega_N^{1k}+x(2)\omega_N^{2k}+...+x(N-1)\omega_N^{(N-1)k}

然后我们可以添加常数为 0 的项来 把这个多项式的项数变成偶数。我们再设两个式子：

\begin{align*} P(\omega_N^k)&=x(0)\omega_N^{0k}+x(2)\omega_N^{1k}+x(4)\omega_N^{2k}+...+x(N-2)\omega_N^{\frac{N-2}{2}k}\\ Q(\omega_N^k)&=x(1)\omega_N^{0k}+x(3)\omega_N^{1k}+x(5)\omega_N^{2k}+...+x(N-1)\omega_N^{\frac{N-2}{2}k} \end{align*}

这时如果我们拿 $(\omega_N^k)^2$ ，即 $\omega_N^{2k}$ 来替换 $P(\omega_N^k)$ 和 $Q(\omega_N^k)$ 的变量，再给 $Q(\omega_N^k)$ 乘上一个 $\omega_N^k$ ，我们又能得到：

\begin{align*} P(\omega_N^{2k})&=x(0)\omega_N^{0k}+x(2)\omega_N^{2k}+x(4)\omega_N^{4k}+...+x(N-2)\omega_N^{(N-2)k}\\ \omega_N^kQ(\omega_N^{2k})&=x(1)\omega_N^{1k}+x(3)\omega_N^{3k}+x(5)\omega_N^{5k}+...+x(N-1)\omega_N^{(N-1)k} \end{align*}

不难发现：

A(\omega_N^k)=P(\omega_N^{2k})+\omega_N^kQ(\omega_N^{2k})

又由 n 次本原单位根性质，得：

A(\omega_N^k)=P(\omega_{\frac{N}{2}}^k)+\omega_N^kQ(\omega_{\frac{N}{2}}^k)

看见这个 $\frac{N}{2}$ ，相信大家应该知道我们的加速思路是往哪里走了吧。

我们考虑二分算法，所以用 $\omega_N^{k+\frac{N}{2}}$ 代替 $\omega_N^k$ ，得：

A(\omega_N^{k+\frac{N}{2}})=P(\omega_{\frac{N}{2}}^k)+\omega_N^{k+\frac{N}{2}}Q(\omega_{\frac{N}{2}}^k)

再利用 n 次单位根性质 3，得：

A(\omega_N^{k+\frac{N}{2}})=P(\omega_{\frac{N}{2}}^k)-\omega_N^kQ(\omega_{\frac{N}{2}}^k)

看！ $A(\omega_N^k)$ 和 $A(\omega_N^{k+\frac{N}{2}})$ 之间只差了一个负号！

所以我们只需要算出 $\omega_N^k$ 、 $P(\omega_{\frac{N}{2}}^k)$ 和 $Q(\omega_{\frac{N}{2}}^k)$ 就能得到 $A(\omega_N^k)$ 和 $A(\omega_N^{k+\frac{N}{2}})$ 。

又因为 P 多项式的系数是 A 多项式的偶次项，Q 多项式的系数是 A 多项式的奇次项，所以我们可以拆分 A 多项式的系数，就能得到 P、Q 多项式。而且由于 P、Q 多项式的构型又与 A 多项式完全一致，所以我们求解 P、Q 多项式的时候又可以把它们当成新的 A 多项式来求解。这就是一个标准的二分算法。

至此，我们成功把时间复杂度降到了 $O(n\log_2n)$ ，把 DFT 变成了 FFT。我们可以基于这种递归思路来实现 FFT 了，代码如下：

void FFT(std::complex<double>* src, int n) {
  if (n == 1) // 只有一个元素，无需计算
    return;
  int mid = n >> 1;
  static std::complex<double> dst[MAXN];
  for (int i = 0; i < mid; ++i) { // 前半为偶次项，后半为奇次项
    dst[i] = src[i << 1];
    dst[i + mid] = src[i << 1 | 1];
  }
  memcpy(src, dst, sizeof(dst));
  FFT(src, n >> 1); // 计算偶次项
  FFT(src + mid, n >> 1); // 计算奇次项
  for (int i = 0; i < mid; ++i) { // 计算原式
    // n次本原单位根的i次幂
    std::complex<double> w(cos(2 * PI * i / n), -sin(2 * PI * i / n));
    dst[i] = src[i] + w * src[i + mid];
    dst[i + mid] = src[i] - w * src[i + mid];
  }
  memcpy(src, dst, sizeof(dst));
}

快速傅里叶逆变换（IFFT）

公式为：

x(n)=\mathrm{IDFT}[X(k)]=\frac{1}{N}\sum_{n=0}^{N-1}X(k)e^{i\frac{2\pi n}{N}k}

我们发现 FFT 和 IFFT 的不同仅为：

把 FFT 的输出当成输入给 IFFT 计算
IFFT 的输出需要乘一个系数 $\frac{1}{N}$
IFFT 的 n 次本原单位根取的是 $\varepsilon_1$

仅此而已，所以我们只需要稍微改动一下 FFT 的代码，就可以让它同时满足计算 FFT 和 IFFT 的需要：

// op = 1: FFT, op = -1: IFFT
void FFT(std::complex<double>* src, int n, int op) {
  if (n == 1) // 只有一个元素，无需计算
    return;
  int mid = n >> 1;
  static std::complex<double> dst[MAXN];
  for (int i = 0; i < mid; ++i) { // 前半为偶次项，后半为奇次项
    dst[i] = src[i << 1];
    dst[i + mid] = src[i << 1 | 1];
  }
  memcpy(src, dst, sizeof(dst));
  FFT(src, n >> 1); // 计算偶次项
  FFT(src + mid, n >> 1); // 计算奇次项
  for (int i = 0; i < mid; ++i) { // 计算原式
    // n次本原单位根的i次幂，op控制共轭
    std::complex<double> w(cos(2 * PI * i / n), -op * sin(2 * PI * i / n));
    dst[i] = src[i] + w * src[i + mid];
    dst[i + mid] = src[i] - w * src[i + mid];
  }
  memcpy(src, dst, sizeof(dst));
}

当然在递归函数里我们没法把最后一步乘以 $\frac{1}{N}$ 解决，所以需要你在函数外部计算。

还可以更快？

当然！因为这份代码的常数有点大，我们需要用一些技巧来让它更快。

不要用 `std::complex`

有些情况下 STL 也是会被卡的，特别是在没开 O2 优化时。我们最好自己实现一个复数类，只需要实现乘法、加法和减法操作即可。

struct cpx {
  double r, i;
  cpx(double a = 0.0, double b = 0.0): r(a), i(b) {}
  cpx operator+(const cpx& tmp) {
    return cpx(r + tmp.r, i + tmp.i);
  }
  cpx operator-(const cpx& tmp) {
    return cpx(r - tmp.r, i - tmp.i);
  }
  cpx operator*(const cpx& tmp) {
    return cpx(r * tmp.r - i * tmp.i, r * tmp.i + i * tmp.r);
  }
};

把递归变成非递归

不像普通的二分算法，只需要移动计算区间即可，FFT 的计算是专门针对奇偶项的，每次二分会对数组进行重排，很明显不能直接展开。

所以我们还是动手实际分析一下。

当 n=2 时，是我们递归操作的最小单元，此时可知：

\begin{align*} X[0]&=x[0]+\omega_2^0\cdot x[1]\\ X[1]&=x[0]-\omega_2^0\cdot x[1] \end{align*}

把上式转换成一张图，如下：

我们能很直观地看出 $X[0]$ 和 $X[1]$ 都由 $x[0]$ 和 $x[1]$ 转化而来，且哪个乘上了 $\omega_2^0$ 也知道了。我们称这种结构为 蝶形结，这也是为什么这种变换叫蝴蝶变换 —— 它看起来就像一只蝴蝶。

如果你可以理解上面这张图，那么我们现在更上一层，看看 n=4 的情况：

左半部分还是我们熟悉的 n=2 的情况，但是似乎 下标顺序 好像有些不同？

我们再看右半部分：中间部分就是下标 0~3 分成了 偶数和奇数部分！

没错，这个中间层就是第一次奇偶分项得到的下标结果。同样，右半部分也是两项得一项的原式计算部分。

我们再来看看 n=8 的情况：

有点眼花？我们还是一层一层来，看右边的中间层，很明显下标也是奇偶分项的；再看左边的中间层：它被分成了两个小单元，每个单元依旧是上一层按 项序号顺序（不是下标）奇偶分项。

把整张图从右往左看，就是整个递归分项的过程；从左往右看，就是回溯计算的过程。

现在我们得到了最左边（最底层）的下标顺序，那么我们就可以用循环代替递归来逐层向上计算，得到最终的 $X[k]$ 了！

但是好像这个下标顺序好像也看不出什么端倪来……我们把它们转化成二进制：

\begin{align*} 0:(0)_{10}&\to(000)_2\\ 1:(4)_{10}&\to(100)_2\\ 2:(2)_{10}&\to(010)_2\\ 3:(6)_{10}&\to(110)_2\\ 4:(1)_{10}&\to(001)_2\\ 5:(5)_{10}&\to(101)_2\\ 6:(3)_{10}&\to(011)_2\\ 7:(7)_{10}&\to(111)_2 \end{align*}

看起来还是没啥端倪，那我们把序号也一起二进制化：

\begin{align*} (000)_2:(0)_{10}&\to(000)_2\\ (001)_2:(4)_{10}&\to(100)_2\\ (010)_2:(2)_{10}&\to(010)_2\\ (011)_2:(6)_{10}&\to(110)_2\\ (100)_2:(1)_{10}&\to(001)_2\\ (101)_2:(5)_{10}&\to(101)_2\\ (110)_2:(3)_{10}&\to(011)_2\\ (111)_2:(7)_{10}&\to(111)_2 \end{align*}

好像下标就是序号的二进制倒位！于是我们能得到新的算法代码：

// 要求n为2的倍数，len为n的位数-1
int n, len;
int rvs[n];
// 递推计算二进制倒位
for (int i = 1; i < n; ++i)
  rvs[i] = (rvs[i >> 1] >> 1) | ((i & 1) << (len - 1));
// FFT(op == true) or IFFT(op == false)
void FFT(cpx* A, int n, int len, bool op) {
  // Reverse
  for (int i = 0; i < n; ++i)
    if (i < rvs[i])
      std::swap(A[i], A[rvs[i]]);
  // Start calculating
  for (int i = 1; i < n; i <<= 1) {
    // FFT need sin(PI / i), but IFFT need -sin(PI / i)
    cpx wn(cos(PI / i), sin(PI / i) * (op ? 1.0 : -1.0));
    for (int j = 0; j < n; j += (i << 1)) {
      cpx w(1.0, 0.0);
      for (int k = 0; k < i; ++k, w = w * wn) {
        cpx x = A[j + k], y = w * A[i + j + k];
        A[j + k] = x + y;
        A[i + j + k] = x - y;
      }
    }
  }
  // If is IFFT
  if (!op)
    for (int i = 0; i < n; ++i) {
      A[i].r /= n;
      A[i].i /= n;
    }
}

至此，您已经学会了快速傅里叶算法！

傅里叶变换长什么样

从傅里叶级数到傅里叶变换