傅里叶变换可以完成时域与频域的相互转换，在竞赛中常用于大数乘法与多项式乘法。

傅里叶变换长什么样

（连续）傅里叶变换/逆变换公式如下：

其中 F(ω)F(\omega) 为 f(t)f(t) 的象函数，f(t)f(t) 为 F(ω)F(\omega) 的象原函数

这是一个无穷区间反常积分，众所周知计算机是没法处理连续区间数据的，所以计算机使用的是另一种傅里叶变换 —— 离散傅里叶变换 （DFT）：

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

离散傅里叶变换本质上就是从象原函数中 等距抽取有限数量 的采样点，将积分转换成求和处理。
而我们的快速傅里叶变换 （FFT） 就是基于离散傅里叶变换算法的优化加速。

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

傅里叶变换源于法国数学家及物理学家 —— 让·巴普蒂斯·约瑟夫·傅里叶 的发现：

能将满足一定条件的某个函数表示成正弦/余弦函数或者它们的积分的线性组合

这就是 傅里叶级数，也是高数书上所告诉我们的：

由欧拉公式：

我们可以把上面这个式子写成更加简单的复数形式：

以上是对于 f(x)f(x) 是一个周期为 2π2\pi 的函数而言的，对于周期为 2l2l 的函数，只需要进行一次放缩变换即可：

傅里叶变换是连续的，而傅里叶级数仍然是离散的，所以我们要先设法把傅里叶级数给变成连续的（当然变成连续的之后就不能叫傅里叶级数了）。

我们现在重新思考一下变化的 nn 代表了什么。因为 nn 会在正弦 / 余弦函数中以 xx 的系数出现，也就是作为正弦 / 余弦函数的 角速度，由于角速度与频率之间有关系：

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

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

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

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

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

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

离散傅里叶变换

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

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

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

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

n次单位根

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

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

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

性质 1

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

性质 2

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

性质 3

推导：

n 次本原单位根

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

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

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

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

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

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

推导：

由 DFT 到 FFT

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

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

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

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

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

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

不难发现：

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

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

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

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

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

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

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

至此，我们成功把时间复杂度降到了 O(nlog⁡2n)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）

公式为：

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

1. 把 FFT 的输出当成输入给 IFFT 计算
2. IFFT 的输出需要乘一个系数 1N\frac{1}{N}
3. IFFT 的 n 次本原单位根取的是 ε1\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));
}

当然在递归函数里我们没法把最后一步乘以 1N\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 时，是我们递归操作的最小单元，此时可知：

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

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

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

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

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

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

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

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

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

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

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

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

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

// 要求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;
    }
}

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