快速傅里叶算法

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

傅里叶变换长什么样

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

F (ω) = F [f (t)] f (t) = F^{- 1} [F (ω)] = \int_{- \infty}^{+ \infty} f (t) e^{- iω t} d t = \frac{1}{2 π} \int_{- \infty}^{+ \infty} F (ω) e^{iω t} d ω

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

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

X (k) = DFT [x (n)] x (n) = IDFT [X (k)] = n = 0 \sum N - 1 x (n) e^{- i \frac{2 πn}{N} k} = \frac{1}{N} n = 0 \sum N - 1 X (k) e^{i \frac{2 πn}{N} k}

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

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

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

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

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

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

f (x) a_{0} a_{n} b_{n} = \frac{a _{0}}{2} + n = 1 \sum \infty [a_{n} cos (n x) + b_{n} sin (n x)] = \frac{1}{π} \int_{- π}^{π} f (x) d x = \frac{1}{π} \int_{- π}^{π} f (x) cos (n x) d x = \frac{1}{π} \int_{- π}^{π} f (x) sin (n x) d x

由欧拉公式：

cos x sin x = \frac{e ^{i x} + e ^{- i x}}{2} = \frac{e ^{i x} - e ^{- i x}}{2 i}

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

f (x) c_{n} = \frac{a _{n} - i b _{n}}{2} = n = - \infty \sum \infty c_{n} e^{in x} = \frac{1}{2 π} \int_{- π}^{π} f (x) e^{- in x} d x

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

f (x) c_{n} = n = - \infty \sum \infty c_{n} e^{i \frac{nπ x}{l}} = \frac{1}{2 l} \int_{- l}^{l} f (x) e^{- i \frac{nπ x}{l}} d x

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

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

f = \frac{ω}{2 π} = \frac{n}{2 π}

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

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

f (x) F (t) = \int_{- \infty}^{\infty} F (t) e^{i t x} d t = \frac{1}{2 π} \int_{- π}^{π} f (x) e^{- i t x} d x

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

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

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

F (t) f (x) = \int_{- \infty}^{\infty} f (x) e^{- i t x} d x = \frac{1}{2 π} \int_{- \infty}^{\infty} F (t) e^{i t x} d t

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

离散傅里叶变换

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

X (k) = DFT [x (n)] x (n) = IDFT [X (k)] = n = 0 \sum N - 1 x (n) e^{- i \frac{2 πn}{N} k} = \frac{1}{N} n = 0 \sum N - 1 X (k) e^{i \frac{2 πn}{N} k}

其中 $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 πk}{n}}, k = (0, 1, 2, ..., n - 1)

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

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

性质 1

ε_{k} = ε_{k % n}

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

性质 2

ε_{k} = (ε_{1})^{k}

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

性质 3

ε_{k + \frac{n}{2}} = - ε_{k}

推导：

ε_{k + \frac{n}{2}} = e^{i \frac{2 π ( k + \frac{n}{2} )}{n}} = e^{i \frac{2 πk + πn}{n}} = e^{i \frac{2 πk}{n}} \cdot e^{iπ} = - e^{i \frac{2 πk}{n}} = - ε_{k}

n 次本原单位根

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

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

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

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

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

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

ω_{N}^{d k} = ω_{\frac{N}{d}}^{k}

推导：

ω_{N}^{d k} = e^{i \frac{2 π}{N} d k} = e^{i \frac{2 π}{\frac{N}{d}} k} = ω_{\frac{N}{d}}^{k}

由 DFT 到 FFT

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

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

X (k) = n = 0 \sum N - 1 x (n) ω_{N}^{nk}

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

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

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

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

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

P (ω_{N}^{k}) Q (ω_{N}^{k}) = x (0) ω_{N}^{0 k} + x (2) ω_{N}^{1 k} + x (4) ω_{N}^{2 k} + ... + x (N - 2) ω_{N}^{\frac{N - 2}{2} k} = x (1) ω_{N}^{0 k} + x (3) ω_{N}^{1 k} + x (5) ω_{N}^{2 k} + ... + x (N - 1) ω_{N}^{\frac{N - 2}{2} k}

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

P (ω_{N}^{2 k}) ω_{N}^{k} Q (ω_{N}^{2 k}) = x (0) ω_{N}^{0 k} + x (2) ω_{N}^{2 k} + x (4) ω_{N}^{4 k} + ... + x (N - 2) ω_{N}^{(N - 2) k} = x (1) ω_{N}^{1 k} + x (3) ω_{N}^{3 k} + x (5) ω_{N}^{5 k} + ... + x (N - 1) ω_{N}^{(N - 1) k}

不难发现：

A (ω_{N}^{k}) = P (ω_{N}^{2 k}) + ω_{N}^{k} Q (ω_{N}^{2 k})

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

A (ω_{N}^{k}) = P (ω_{\frac{N}{2}}^{k}) + ω_{N}^{k} Q (ω_{\frac{N}{2}}^{k})

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

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

A (ω_{N}^{k + \frac{N}{2}}) = P (ω_{\frac{N}{2}}^{k}) + ω_{N}^{k + \frac{N}{2}} Q (ω_{\frac{N}{2}}^{k})

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

A (ω_{N}^{k + \frac{N}{2}}) = P (ω_{\frac{N}{2}}^{k}) - ω_{N}^{k} Q (ω_{\frac{N}{2}}^{k})

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

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

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

至此，我们成功把时间复杂度降到了 $O (n lo g_{2} n)$ ，把 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) = IDFT [X (k)] = \frac{1}{N} n = 0 \sum N - 1 X (k) e^{i \frac{2 πn}{N} k}

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

把 FFT 的输出当成输入给 IFFT 计算
IFFT 的输出需要乘一个系数 $\frac{1}{N}$
IFFT 的 n 次本原单位根取的是 $ε_{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 时，是我们递归操作的最小单元，此时可知：

X [0] X [1] = x [0] + ω_{2}^{0} \cdot x [1] = x [0] - ω_{2}^{0} \cdot x [1]

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

寒夏摸鱼站

Live your dream, and share your passion.