算法分析-分治

1. 渐近符号

• $O$： "Big Oh"，$f = O(g)$ 表示，$\exists \space c > 0$，使得 $\forall \space n$，$f(n) \le c\cdot g(n)$。相当于 $f\le g$。
• $o$："Little Oh"，$f=o(g)$ 表示，$\lim_{n\rightarrow\infty}\frac{f(n)}{g(n)} = 0$，相当于 $f < g$ 。
• $\Omega$："Big Omega"，$f=\Omega (g)$ 表示 $g = O(f)$，相当于 $f\ge g$。
• $\omega$："Little Omega"，$f=\omega(g)$ 表示 $f = o(f)$，相当于 $f>g$。
• $\Theta$："Theta"，$f = \Theta(g)$，表示 $f = O(g)$ 且 $f=\Omega(g)$，相当于 $f = g$。

$O$ 还可以使用 $\frac{f(n)}{g(n)}$ 的 lim sup 定义。lim sup 直观的感受是，当趋近于正无穷的时候，函数值的上界会越来越小（如果存在 lim sup 的话），最后这个上界收敛到的值就是该函数的 lim sup。$\frac{f(n)}{g(n)}$ 存在 lim sup 的话，就有 $f = O(g)$ 的关系。

不过排除掉震荡函数和其他什么邪门的函数，如果我们直接算 $\lim_{n\rightarrow\infty}\frac{f(n)}{g(n)}$ 发现是存在极限的，那么必然可以得到 $f = O(g)$（注意我们把极限等于 $0$ 的情况也包含进去了，因为对于 $f$，满足 $o$ 关系的函数 $g$ 必然也能满足 $O$。

2. 主定理

对于 $T(n) = a\cdot T(\frac{n}{b}) + f(n)$，其中 $a$ 和 $b$ 均为常数，我们有如下定理：

• 情况一，如果存在 $\epsilon > 0$，满足 $f(n) = O(n^{(\log_ba) - \epsilon})$，则有 $T(n) = \Theta(n^{\log_b a})$。
• 情况二，如果存在 $k \ge 0$，满足 $f(n) = \Theta(n^{\log_b a}\log^kn)$，则有 $T(n)=\Theta(n^{\log_ba}\log^{k+1}n)$。
• 情况三，如果存在 $\epsilon > 0$，满足 $f(n) = \Omega(n^{(\log_ba)+\epsilon})$，同时满足正则条件，即存在 $c < 1$，满足当 $n$ 足够大时， $a\cdot f(\frac{n}{b})\le c \cdot f(n)$，则有 $T(n)=\Theta(f(n))$。

分治算法的时间复杂度一般使用递归树来分析。如果某个式子不满足主定理，那么可以考虑从递归树上进行推导。

同样地，主定理也可以使用递归树来理解。递归树的时间复杂度贡献来源有两种，一种是叶子节点的贡献，一种是非叶子节点的贡献。

对于叶子节点的贡献，递归树每层往下都分 $a$ 个叉，一共是 $\log_bn$ 层。叶子节点的时间复杂度贡献只有 $\Theta(1)$，所以全部叶子节点的贡献为 $\Theta(a^{\log_bn}) = \Theta(n^{\log_ba})$。

对于非叶子节点的贡献，稍微准确一点地来算应该是：

$$\sum_{j=0}^{\log_bn - 1}a^jf(\frac{n}{b^j})$$

这大概意义上来说，有：（这里不保对，仅供用于后续理解）

$$\frac{a^{j+1}\cdot f(\frac{n}{b^{j+1}})}{a^j\cdot f(\frac{n}{b^j})} = a\cdot \frac{f(\frac{n}{b^j\cdot b})}{f(\frac{n}{b^j})} \approx a\cdot \frac{f(\frac{n}{b})}{f(n)} = r$$

然后非叶子节点的贡献就是 $p\cdot f(n)$，$p = \sum_{i=0}^{\log_bn-1} r^i$。

然后对于情况一，$f(n)=O(n^{\log_ba})$，此时 $f(n)$ 都比较小了，$p$ 往往不能把 $f(n)$ 放大太多，所以非叶节点的贡献小于叶子节点的贡献，最终的时间复杂度就由叶子节点决定，为 $\Theta(n^{\log_ba})$。

对于情况二，如果 $f(n) = \Theta(n^{\log_b a}\log^kn)$，说明 $f(n)$ 和叶子节点的贡献差不多。同时我们有如下观察：

$$\begin{align} a\cdot f(\frac{n}{b}) &= \Theta(a\cdot (\frac{n}{b})^{\log_ba}\log^k(\frac{n}{b})) \\ &= \Theta(a\cdot \frac{n^{\log_ba}}{a}\log^k(\frac{n}{b})) \\ &= \Theta(n^{\log_ba}\log^kn) \\ &= f(n) \end{align}$$

所以每一层的贡献都是一样的。总共有 $\Theta(\log n)$ 层，因此总时间复杂度为 $T(n)=\Theta(n^{\log_ba}\log^{k+1}n)$。

对于情况三，如果 $f(n) = \Omega(n^{\log_ba})$，说明单个 $f(n)$ 就比叶子节点的贡献大了。理论上总时间复杂度应该是 $\Theta(p\cdot f(n))$ 的。主定理又加了一个限制，要求 $a\cdot f(\frac{n}{b})\le c \cdot f(n)$，也就是 $a\cdot \frac{f(\frac{n}{b})}{f(n)} = r \le c < 1$，这就使得 $p$ 是对一个公比小于 $1$ 的等比数列求和，求和结果是收敛的，所以总时间复杂度是 $\Theta(f(n))$ 。不过我认为，实际上这个限制只是用来卡掉一些特殊情况，比如不太连续的那种邪门函数。对于一般的函数这个条件是自然满足的，毕竟如果不收敛的话，一致发散到最后一层叶子节点，很容易使得叶子节点的贡献大于 $f(n)$，从而使得情况三的前提就不成立了。

然后谈一下情况一和情况三中的 $\epsilon$。我的理解是，这是为了让 $O$ 和 $\Omega$ 与 $\Theta$ 不重叠。但是主定理没选择用 $o$ 和 $\omega$ 来限制，而是使用 $\epsilon$，这是因为 $o$ 和 $\omega$ 只是说和 $\Theta$ 有差距即可，有什么差距不做限制，而主定理的 $\epsilon$ 要求必须是差多项式因子的差距。

比如，如果 $f(n) = \frac{n^{\log_ba}}{\log n}$，那么 $\lim_{n\rightarrow\infty}\frac{\frac{n^{\log_ba}}{\log n}}{n^{(\log_ba)-\epsilon}}=\lim_{n\rightarrow\infty}\frac{n^\epsilon}{\log n}$，我们发现根本找不到一个 $\epsilon > 0$ 使得这个极限存在。

3. 算法

3.1 Karatsuba 算法

主要用于高精度乘法。

现在有两个长度为 $n$ 的大数 $X$ 和 $Y$，将其分别拆分成两个部分，得到：

$$X=x_{high}\cdot 10^\frac{n}{2} + x_{low} \\ Y = y_{high}\cdot 10^\frac{n}{2} + y_{low}$$

于是：

$$\begin{align} XY &= (x_{high}\cdot 10^\frac{n}{2} + x_{low})(y_{high}\cdot 10^\frac{n}{2} + y_{low}) \\ &= x_{high}y_{high}\cdot 10^n+(x_{high}y_{low}+x_{low}y_{high})\cdot 10^\frac{n}{2}+x_{low}y_{low} \end{align}$$

然后注意到：

$$\begin{align} &(x_{high}+x_{low})(y_{high}+y_{low})\\ &= x_{high}y_{high}+(x_{high}y_{low}+x_{low}y_{high})+x_{low}y_{low} \end{align}$$

$$\begin{align} & x_{high}y_{low}+x_{low}y_{high} \\ & = (x_{high}+x_{low})(y_{high}+y_{low}) - x_{high}y_{high} - x_{low}y_{low} \end{align}$$

所以我们只需要计算 $(x_{high}+x_{low})(y_{high}+y_{low})$ 、$x_{high}y_{high}$、$x_{low}y_{low}$ 这三部分即可。于是有 $T(n) = 3T(\frac{n}{2}) + \Theta (n)$，得到时间复杂度为 $\Theta (n^{\log_23}) \approx \Theta(n^{1.585})$。

在 Python 中，乘法使用 $\Theta(n^2)$ 的小学乘法算法。但是，如果运算符两边的操作数的位数均大于 KARATSUBA_CUTOFF，那么就使用 Karatsuba 算法。KARATSUBA_CUTOFF 被定义成 $70$。

3.2 Fibonacci 数列

定义 Fibonacci 数列的递归公式为：

$$\begin{cases} F_0 = 0 & \\ F_1 = 1 & (n\in N) \\ F_n = F_{n-1} + F_{n-2} & \end{cases}$$

直接递归

有 $T(n) = T(n-1) + T(n-2) + \Theta(1)$。

可以先通过放缩来得到上界和下界。令 $c = \Theta(1)$。

$$\begin{align} T(n) &= T(n-1) + T(n-2) + c \\ &\le 2T(n-1) + c \\ &\le 4T(n-2) + 3c \\ &\le 2^3\cdot T(n-3) + 7c \\ &... \\ &\le 2^n\cdot T(0) + (2^n - 1)c \\ T(n) & = O(2^n) \end{align}$$

$$\begin{align} T(n) &= T(n-1) + T(n-2) + c \\ &\ge 2T(n-2) + c \\ &\ge 4T(n-4) + 3c \\ &\ge 2^3\cdot T(n-6) + 7c \\ &... \\ &\ge 2^{\frac{n}{2}} \cdot T(0) + (2^{\frac{n}{2}} - 1)c \\ T(n) &= \Omega(2^{\frac{n}{2}}) = \Omega((\sqrt{2})^n) \end{align}$$

使用特征方程可以得到紧界。上式的齐次部分为：

$$T(n) - T(n-1) - T(n-2) = 0$$

，特称方程为：

$$r^2 - r - 1 = 0$$

得到两根：$\varphi_1 = \frac{1+\sqrt{5}}{2}$，$\varphi_2=\frac{1-\sqrt{5}}{2}$，齐次通解为：

$$A\cdot \varphi_1^n + B\cdot \varphi_2^n$$

由于非齐次项为常数，所以非齐次特解也为常数。于是非齐次通解为：

$$T(n) = A\cdot \varphi_1^n + B\cdot \varphi_2^n + \Theta(1)$$

注意到，$\left | \varphi_2 \right | < 1$，因此，当 $n \rightarrow \infty$ 时，$A\cdot \varphi_1^n$ 主导，$B\cdot \varphi_2^n \rightarrow 0$，即 $T(n) = \Theta(\varphi_1^n)$。

时间复杂度的式子和递归的式子高度一致。因此 $F_n = \Theta(\varphi_1^n)$ ，$F_n$ 的位数为 $\log_{10} F_n = \Theta(n)$。

这里我们假设了加法操作是常数的时间。但实际上加法操作是 $\Theta(n)$ 的。不过即使考虑了加法操作的时间，由于指数的增长速度很快，最后的时间复杂度还是指数的。

递推

时间复杂度显然是 $\Theta(n)$ 的。

考虑加法操作的耗时，该算法的总耗时是 $O(n^2)$ 的。

矩阵加速

时间复杂度为 $\Theta(\log n)$ 。

但是这中间需要进行矩阵乘法，进行元素相乘。两个 $n$ 位数相乘，只使用 $\Theta(n^2)$ 的小学乘法算法的话，最后得到的总耗时是 $O(n^2\log n)$ 的，看起来还不如递推。

不过我们还可以再精确地计算。矩阵加速时，会计算 $A^1$、$A^2$、$A^4$、$A^8$......，那么矩阵中数字的位数大概是 $1$、$2$、$4$、$8$......这样，然后进行 $O(n^2)$ 的乘法运算，所需时间大概为 $1^2+2^2+4^2+8^2+\dots+n^2$，中间有 $\Theta(\log n)$ 项。最后的结果大概是 $n^2$ 这一项主导，得到总耗时为 $O(n^2)$。

实际上，如果两数相乘的时间复杂度如果为 $M(n)$ 的话，那么总耗时就是 $M(n)$。

3.3 BFPRT

BFPRT 用于在长度为 $n$ 的无序数组中找出前 $k$ 大的数，即 Top k 问题。

朴素的做法是直接排序，时间复杂度是 $\Omega(n\log n)$ 的。

使用堆可以把时间复杂度做到 $O(n+k\log n)$，其中 $O(n)$ 的建堆时间，之后每次从堆中取出一个数需要 $O(\log n)$ 的时间，总共进行 $k$ 次。但是，如果 $k$ 比较靠近中间，比如 $k=\frac{n}{2}$，时间复杂度就又变成了 $O(n\log n)$。

快速选择算法基于快速排序的思想。首先选取一个主元，然后把主元两边的元素交换，将序列随机分成三个部分：小于主元的部分、等于主元的部分（即主元本身）和大于主元的部分，此时主元本身是排好序的了，那么我们可以根据主元的位置来判断接下来递归左右哪个部分。

当主元选取的比较好（恰好是中位数时），每次递归都会排除掉一半的序列，时间复杂度为 $O(n)$。但是如果选取的不好（比如选到最值或者接近最值），那么每次只会排除掉几个，时间复杂度为 $O(n^2)$。随机选取可以做到期望的时间复杂度为 $O(n)$。

BFPRT 算法，也被称为 Median of medians 算法，能够最坏情况下的时间复杂度仍为 $O(n)$。

首先这个算法将序列从左到右，每 $5$ 个元素分成一组，共得到 $\frac{n}{5}$ 组。每组分别找到该组内 $5$ 个元素之间的中位数，这样共产生 $\frac{n}{5}$ 个中位数。然后，我们对这 $\frac{n}{5}$ 个中位数跑一遍当前算法（当前算法就是求 top k，令 $k=\frac{n}{2}$ 就是求中位数），即可得到中位数之间的中位数。然后选择这个中位数作为主元，后续的过程就跟快速选择算法一样了。

取中位数的中位数作为主元，可以确保这个主元的选择不会太差。如上图，如果 $a$ 指向 $b$，说明一定存在 $a\le b$ 的关系。没连边的节点之间的大小关系不确定。每个竖行表示一组。白色节点是该组内的中位数。节点 $x$ 表示是中位数的中位数，也是我们最后要作为主元的节点。

然后我们看到，所有能到达 $x$ 的点都是必然小于等于 $x$ 的，这些点位于图中的阴影区域。这个区域内的点的数量大概为 $\frac{1}{2}\cdot \frac{n}{5}\cdot 3=\frac{3n}{10}$，也就是至少有 30% 的节点是小于等于 $x$ 的。

所有 $x$ 能到达的点也都是必然大于等于 $x$ 的，同理，大概至少有 30% 的节点是大于等于 $x$ 的。

那么我们选到的这个中位数的中位数，每次递归怎么也能排掉 30% 的序列。

于是我们有：

$$T(n) \le T(\frac{n}{5}) + T(\frac{7n}{10}) + \Theta(n)$$

其中 $T(\frac{n}{5})$ 是对中位数找到中位数的时间，$T(\frac{7n}{10})$ 是继续往下递归的时间。

这个形式无法使用主定理，我们只能去猜测这个时间复杂度并带入验证。

假设 $T(n)\le cn$，带入上式，有：

$$\begin{align} T(n) & \le T(\frac{n}{5}) + T(\frac{7n}{10}) + \Theta(n) \\ & \le \frac{cn}{5}+\frac{7cn}{10}+an \\ &= \frac{9cn}{10}+an \\ &= cn + (-\frac{cn}{10}+an) \end{align}$$

只要 $-\frac{cn}{10}+an\le 0$，即可使得 $T(n)\le cn$ 的假设成立。那么也就是说 $c \ge 10a$。对于任意常数 $a$，我们都能满足，于是假设成立，因此 BFPRT 的时间复杂度为 $O(n)$，且是最坏情况下的时间复杂度。

为什么选择 5 个一组

选择其他的数字其实也没多大的问题。不过首先最好是要选择奇数，因为偶数的中位数计算比较麻烦，容易造成偏差。假如我们选择 $3$ 个，同样的推导过程，有：

$$\begin{align} T(n) &= T(\frac{n}{3}) + T(\frac{2n}{3}) + O(n) \\ &\le \frac{cn}{3} + \frac{2cn}{3} + an \\ &= cn + an \end{align}$$

我们发现 $an$ 是非负的，因此 $T(n) \le cn$ 的假设不成立。三个一组不太可行。

实际上只要是大于等于五个一组的都可以，但是五个一组的常数会小一点，所以一般我们选择五个一组。

为什么不一直分组找中位数

既然我们通过了五个一组划分，然后找到了组内的中位数，那么看起来还可以再对这些中位数再进行五个一组这样的操作，这么一直进行下去，毕竟再递归 BFPRT 算法本身去找中位数，其内部还是会进行这么个过程。

实际上分组划分找中位数得到的只是近似中位数，递归调用 BFPRT 找到的是精确的中位数。如果不断地分组划分的话，最后找到的中位数会是反复近似多次的结果，很可能无法保证后续时间复杂度的正确性了。

3.4 平面最近点对

我们充分发扬人类智慧：

将所有点全部绕原点旋转同一个角度，然后按 $x$ 坐标排序

根据数学直觉，在随机旋转后，答案中的两个点在数组中肯定不会离得太远

所以我们只取每个点向后的 $5$ 个点来计算答案

这样速度快得飞起，在 $n=1000000$ 时都可以在 1s 内卡过

当然，这是错的。

我们充分发扬人类智慧：

将所有点全部绕原点旋转同一个角度，然后按 $x\times y$ 排序

根据数学直觉，在随机旋转后，答案中的两个点在数组中肯定不会离得太远

所以我们只取每个点向前的 $50$ 个点来计算答案

这样速度快得飞起，在 $n=400000$ 时都可以在 124ms 内卡过

先把所有点按照 $x_i$ 为第一关键字，$y_i$ 为第二关键字排序，并以点 $p_m$（$m=\left \lfloor \frac{n}{2} \right \rfloor$）把点集分成两个部分，然后对两个部分分别递归下去，求出点集各自内部的最近点对，设距离为 $h_1$ 和 $h_2$，取其中的最小值为 $h$。

然后再考虑跨过两个点集之间的点的最近距离。首先距离大于 $h$ 的肯定不用考虑了，于是我们把所有横坐标与 $x_m$ 的差小于 $h$ 的点放入集合 $B$ 中，即图中绿色的点。

然后为了避免两个点之间相互考虑，对于点 $p_i$，只考虑纵坐标大于等于 $y_i$ 的点，也就是对于图中的点 $P$ ，只考虑红色区域内的 $B$ 的点，即下图中红色的点，这相当于一个剪枝。

（下图中，中间的黑色粗线是左右两部分的分割线。红色部分区域由两个 $h\times h$ 的正方形组成）

然后直接在 $B$ 集合里面暴力枚举就好了。只要在 $y$ 坐标上进行如上这种剪枝，把遍历范围限制住（而不是直接在整个 $B$ 里枚举），时间复杂度就是对的。

为什么对？上图中左半部分的 $h\times h$ 区域里面的点，已经在分治时处理了，所以点之间两两距离必然是 $\ge h$ 的。（是同一正方形部分内的点之间，不是横跨两部分正方形的点之间）

我们把 $h\times h$ 的正方形再分成 $4$ 个 $\frac{h}{2} \times \frac{h}{2}$ 的小正方形，一个小正方形内最多只有一个点，因为小正方形内的最大距离为对角线，即 $\frac{\sqrt{2}h}{2}<h$。

于是，对于点 $P$ 最多只会枚举到 $2\times 4$ 个点。

在具体实现上，为了做上面的剪枝，需要再把集合 $B$ 中的点按 $y$ 坐标排序，于是 $T(n) = 2\cdot T(\frac{n}{2}) + O(n\log n)$，时间复杂度为 $O(n\log^2n)$。

可以使用归并排序。在分治找最小距离的同时再把坐标按 $y$ 排序（两部分合并之后，最开始 $x$ 的排序关系就没用了，再变成 $y$ 的排序是没问题的），这样就是 $T(n)=2\cdot T(\frac{n}{2})+O(n)$，时间复杂度为 $O(n)$。

3.5 棋盘覆盖问题

有一个 $2^k\times 2^k$ 的棋盘，其中有一个方格已经被覆盖了，然后要求你用下图三种东西去把剩余的没被覆盖的方格覆盖掉。

做法是把棋盘横着和竖着分成四部分，然后在没有方格覆盖的三个部分上应用上图的图形，这样搞完之后分成的四个部分都是只有一个方格被覆盖了，化成了规模更小的子问题。

3.6 循环赛日程表

有如下的设计规则：

• 每位选手必须与其他 $n-1$ 位选手各比赛 $1$ 次。
• 每位选手每天只能比赛 $1$ 次。
• 循环赛一共进行 $n-1$ 天。

这个问题可以用分治解决，对于 $2n$ 个选手，可以先把所有他们分成两部分，每部分 $n$ 个人。在前 $n-1$ 天，两个部分同时先自己内部解决，然后后 $n$ 天，两部分之间再相互比赛。

根据这个思想，可以构造出一个日程表：

然后我们发现这个表是有点规律的。实际写代码的时候根据这个规律去生成方案。

3.7 FFT

对于多项式 $f(x)=\sum_{i=0}^{n}a_ix^i,g(x)=\sum_{i=0}^{n}b_ix^i$，定义其乘积 $fg$ 为：

$$(fg)(x)=\left(\sum_{i=0}^{n}a_ix^i\right)\left(\sum_{i=0}^{n}b_ix^i\right)$$

朴素做法是 $O(n^2)$ 的时间复杂度。使用 FFT 可以做到 $O(n\log n)$ 的时间复杂度。

对于 $n$ 次多项式，可以用 $n+1$ 个系数来唯一表示，即系数表示。还可以用 $n+1$ 个多项式上不同的点来表示（插值），即点值表示，正确性可以用范德蒙行列式证明。

点值表示下的多项式可以快速地计算多项式乘积，只需要把相同 $x$ 下的 $y$ 值相乘即可得到多项式相乘后，新的多项式在 $x$ 下的 $y$ 值。

FFT 求多项式乘积的过程如下：

• 将 $f$ 通过快速傅里叶变换，从系数表示转为点值表示。通过取特殊的点，该过程可以做到 $O(n\log n)$。

• 将 $g$ 也按上述的方法从系数表示转为点值表示。

• 然后可以 $O(n)$ 的时间计算出 $fg$ 的点值表示。

• 再通过对 $fg$ 进行快速傅里叶逆变换，从点值表示转为系数表示，该过程耗时 $O(n\log n)$。

单位根

在复平面上画一单位圆，从 $x$ 轴正方向开始逆时针在单位圆上取点，将单位圆平均分成 $n$ 部分，这 $n$ 个点即为 $n$ 个复数。设幅角为正且最小的复数为 $\omega_n$，称为 $n$ 次单位根，即：

$$\omega_{n}=\cos \frac{2\pi}{n}+{\rm i}\sin \frac{2\pi}{n}$$

单位根有如下性质：

• 由欧拉公式可得：

  $$\omega_{n}^{k}=\cos \frac{2k\pi}{n}+{\rm i}\sin \frac{2k\pi}{n}\\$$

• $\omega_n^0=\omega_n^n=1$

• $\omega_{rn}^{rk} = \omega_n^k$

• $\omega _{n}^{k+\frac{n}{2}} = -\omega _{n}^{k}$

• $\overline{\omega_n^k} = \omega_n^{n-k}$

快速傅里叶变换

任何多项式都可以化成 $n$ 次多项式 $f(x)$，其中 $n$ 为 $2$ 的幂次（只需要在原来的多项式后面加系数为 $0$ 的高次项即可）。

然后按下标的奇偶性，把多项式中的项分成两部分：

$$\begin{align}f(x)= \space &(a_0+a_2{x^2}+a_4{x^4}+\dots+a_{n-2}x^{n-2})\\&+(a_1x+a_3{x^3}+a_5{x^5}+ \dots+a_{n-1}x^{n-1})\end{align}\\$$

令

$$\begin{align} f_1(x)&=a_0+a_2{x^1}+a_4{x^2}+\dots+a_{n-2}x^{\frac{n}{2}-1}\\ f_2(x)&=a_1+a_3{x}+a_5{x^2}+ \dots+a_{n-1}x^{\frac{n}{2}-1} \end{align}\\$$

则

$$f(x)=f_1(x^2)+xf_2(x^2)\\$$

现在我们想快速地求出 $f(w_n^i)$ （$i=0\dots n-1$）来得到多项式对应的点值表示。我们把要求的这个分成两部分：$\omega_n^k$ 和 $\omega_n^{k+\frac{n}{2}}$（$k<\frac{n}{2}$）。

带入 $x=\omega_n^k$（$k<\frac{n}{2}$）可得

$$\begin{align}f(\omega_n^k)&=f_1(\omega_n^{2k})+\omega_n^kf_2(\omega_n^{2k})\\ &=f_1(\omega_{\frac{n}{2}}^{k})+\omega_n^kf_2(\omega_{\frac{n}{2}}^{k}) \end{align}\tag1$$

带入 $x=\omega_n^{k+\frac{n}{2}}$（$k<\frac{n}{2}$）可得

$$\begin{align}f(\omega_n^{k+\frac{n}{2}})&=f_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}f_2(\omega_n^{2k+n})\\ &=f_1(\omega_n^{2k}\cdot\omega_n^n)-\omega_n^kf_2(\omega_n^{2k}\cdot\omega_n^n)\\&=f_1(\omega_n^{2k})-\omega_n^kf_2(\omega_n^{2k}) \\ &=f_1(\omega_{\frac{n}{2}}^{k})-\omega_n^kf_2(\omega_{\frac{n}{2}}^{k}) \end{align}\tag2$$

于是我们只需要求出 $f_1(\omega_{\frac{n}{2}}^{k})$ 和 $f_2(\omega_{\frac{n}{2}}^{k})$ 即可得到 $f(\omega_n^k)$ 和 $f(\omega_n^{k+\frac{n}{2}})$，相当于只需要求解两个规模砍半的子问题，$T(n)=2\cdot T(\frac{n}{2})+O(n)$，得到时间复杂度为 $O(n\log n)$。

快速傅里叶逆变换

傅里叶变换的过程可以用矩阵表示：

$$\begin{bmatrix}1 & 1 & 1 & \cdots & 1 \\1 & (\omega_n^1)^1 & (\omega_n^1)^2 & \cdots & (\omega_n^1) ^ {n-1}\\1 & (\omega_n^2)^1 & (\omega_n^2)^2 & \cdots & (\omega_n^2) ^ {n-1} \\ \vdots & \vdots &\vdots & \ddots & \vdots \\ 1 & (\omega_n^{n-1})^1 & (\omega_n^{n-1})^2 & \cdots & (\omega_n^{n-1}) ^ {n-1}\end{bmatrix} \begin{bmatrix}a_0 \\ a _ 1 \\ a _ 2 \\ \vdots \\ a _ {n-1}\\ \end{bmatrix} = \begin{bmatrix}y_0 \\ y _ 1 \\ y _ 2 \\ \vdots \\ y _ {n-1}\\ \end{bmatrix}\quad\\$$

然后我们看这个公式：

$$x^n-1=(x-1)(x^{n-1}+x^{n-2}+\cdots+1)\\$$

• 当 $x=1$ 时，$x^{n-1}+x^{n-2}+\cdots+1 = n$
• 当 $x\neq1$ 时， $x^{n-1}+x^{n-2}+\cdots+1 = \frac{x^n-1}{x-1}$

然后我们有：

$$\begin{align}\begin{bmatrix}1&(\omega_n^{n-i})^1&({\omega_n^{n-i}})^2&\cdots&(\omega_{n}^{n-i})^{n-1}\end{bmatrix}\begin{bmatrix}1\\\ (\omega_n^{j})^1\\(\omega_n^{j})^2\\\vdots\\(\omega_n^{j})^{n-1}\end{bmatrix}&=\sum_{k=0}^{n-1}(\omega_n^{n-i})^k(\omega_n^j)^k\\ &=\sum_{k=0}^{n-1}(\omega_n^{n-i+j})^k \end{align}$$

套用上面的公式：

• 当 $\omega_n^{n-i+j} = 1$，即 $i=j$ 时，上式为 $n$
• 当 $\omega_n^{n-i+j} \neq 1$，即 $i\neq j$ 时，$(\omega_n^{n-i+j})^n = \omega_n^{n(n-i-j)}=\omega_1^{n-i-j}=1$，于是上式为 $0$。

这就表明

$$\mathbf C=\frac1n\begin{bmatrix}1&1&1&\cdots&1\\ 1&(\omega_n^{n-1})^1&({\omega_n^{n-1}})^2&\cdots&(\omega_{n}^{n-1})^{n-1}\\ 1&(\omega_n^{n-2})^1&({\omega_n^{n-2}})^2&\cdots&(\omega_{n}^{n-2})^{n-1}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&(\omega_n^1)^1&{(\omega_n^1)}^2&\cdots&(\omega_{n}^1)^{n-1}\end{bmatrix}\\$$

为

$$\begin{bmatrix}1 & 1 & 1 & \cdots & 1 \\1 & (\omega_n^1)^1 & (\omega_n^1)^2 & \cdots & (\omega_n^1) ^ {n-1}\\1 & (\omega_n^2)^1 & (\omega_n^2)^2 & \cdots & (\omega_n^2) ^ {n-1} \\ \vdots & \vdots &\vdots & \ddots & \vdots \\ 1 & (\omega_n^{n-1})^1 & (\omega_n^{n-1})^2 & \cdots & (\omega_n^{n-1}) ^ {n-1}\end{bmatrix}\\$$

的逆矩阵。

根据单位根的性质：

$$\begin{align} \mathbf C&=\frac1n\begin{bmatrix}1&1&1&\cdots&1\\ 1&(\omega_n^{n-1})^1&({\omega_n^{n-1}})^2&\cdots&(\omega_{n}^{n-1})^{n-1}\\ 1&(\omega_n^{n-2})^1&({\omega_n^{n-2}})^2&\cdots&(\omega_{n}^{n-2})^{n-1}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&(\omega_n^1)^1&{(\omega_n^1)}^2&\cdots&(\omega_{n}^1)^{n-1}\end{bmatrix}\\ &=\frac1n\begin{bmatrix}1 & 1 & 1 & \cdots & 1 \\1 & (\overline{\omega_n^1})^1 & (\overline{\omega_n^1})^2 & \cdots & (\overline{\omega_n^1}) ^ {n-1}\\1 & (\overline{\omega_n^2})^1 & (\overline{\omega_n^2})^2 & \cdots & (\overline{\omega_n^2}) ^ {n-1} \\ \vdots & \vdots &\vdots & \ddots & \vdots \\ 1 & \overline{(\omega_n^{n-1}})^1 & (\overline{\omega_n^{n-1}})^2 & \cdots & (\overline{\omega_n^{n-1}}) ^ {n-1}\end{bmatrix}\\\end{align}\\$$

也就是说，在做快速傅里叶逆变换的时候，只需要在做快速傅里叶变换的基础上，把本来要带入的单位根取共轭复数，然后再把变换的结果乘 $\frac{1}{n}$ 就可以了。