快速傅里叶变换

多项式乘法

记多项式$F(x), G(x)$进行乘法，得到一个新的多项式$H(x)$，记$C_{Poly}(i)$表示$Poly$这个多项式$x_i$这一项的系数是多少。

那么$C_H(i) = \sum_{j=0}^i C_F(j)C_G( i-j )$
我们也称这种形式的乘法为卷积。

朴素地计算$H(x)$的时间复杂度显然是$O(n^2)$的。
但是我们可以通过点值和插值运算，得到新多项式的值的信息，再得到新的多项式。

---

点值与插值

点值运算

一个$n$次多项式$f(x)$以由$n+1$个点值对$(x_0, f(x_0)), (x_1, f(x_1)), \cdots,( x_n, f( x_n ))$来表示，其中$x_i$各不相同。
因为$n+1$个点值对确定的$n$次多项式是唯一的。

证明
记系数列向量$A^T=( c_0, c_1, \cdots, c_n )$， 行向量$B=( f( x_0 ), f( x_1 ), \cdots, f( x_n ) )$
那么$B=V_nA$，其中$V_n$是$n$阶范德蒙德矩阵。
而范德蒙德矩阵可逆的充要条件是$x_i$互不相同。那么这里就可以求出唯一的逆矩阵$V^{-1}_n$，以及唯一的系数列向量$A=BV^{-1}_n$

插值运算

插值运算是点值运算的逆运算。
用矩阵形式来表达则是$A=BV^{-1}_n$

应用

假如我们可以通过某种手段获得两个多项式的点值表达，则他们的积的多项式的点值表达则可以$O(n)$的求出来。假如我们可以通过某种手段实现插值运算，则多项式问题就迎刃而解了。

注意这里假如两个多项式的次数为$n$，那么我们需要求$2(n+1)$个点值对。因为乘出来的多项式次数是$2n$的。

---

N次复数根

定义

主$n$次单位复数根是满足$w^n=1$的复数$w$，记为$w_n$
$n$次复数根恰好有$n$个，从图像上看它们平分了复平面上的单位圆。
数值上这些根为
$w_n^k=e^{2\pi \cdot \frac{k}{n}}$
其中$k=0, 1, \cdots, n-1$
特别的，$w_n^0=1$

欧拉幅角公式

$e^{i\varphi}=\cos \varphi + i\sin \varphi$

证明
$e^x=1+x+\frac{x^2}{2!}+\frac{x^3}{3!}\cdots$
$\sin x=x-\frac{x^3}{3!}+\frac{x^5}{5!}\cdots$
$\cos x=1-\frac{x^2}{2!}+\frac{x^4}{4!}\cdots$

$e^{i\varphi}=1+i\varphi-\frac{{\varphi}^2}{2!}-\frac{i\varphi^3}{3!}+\frac{\varphi^4}{4!}+\cdots$
$\cos \varphi=1-\frac{\varphi^2}{2!}+\frac{\varphi^4}{4!}\cdots$
$i\sin \varphi=i\varphi-\frac{i\varphi^3}{3!}+\frac{i\varphi^5}{5!}\cdots$
$e^{i\varphi}=\cos \varphi + i\sin \varphi$

性质

以下的性质中可以将$n$视作$2$的幂，因为实际操作时我们也是这么做的。

复数的乘法对应旋转，又由欧拉幅角公式可以看出这些复数的模长都为$1$，于是就有

$w_n^k=w_n^{ k-1 }w_n=(w_n)^k$

因为$w_n^n=w_n^0=1$，故

• $w_n^p=w_n^{p \mod n}$
• $w_n^aw_n^b=w_n^{(a+b)\mod n}$

消去引理：$w_{dn}^{dk}=w_n^k$

证明
$w_{dn}^{dk}=e^{2\pi \frac{dk}{dn}}=e^{2\pi \frac{k}{n}}=w_{n}{k}$

推论：

• $w_n^{\frac{n}{2}}=w_{2}=-1$
• $w_n^{\frac{n}{2}+k}=-w_n^k$
• $(w_n^k)^2=w_n^{2k}=w_{\frac{n}{2}}^k$

折半引理

$n$次复数根的平方组成的集合恰好是$\frac{n}{2}$次复数根组成的集合

证明
$w_{\frac{n}{2}}^k=( w_n^k )^2=( w_n^{ k+\frac{n}{2}})^2$

求和引理

$\sum_{j=0}^{n-1}(w_n^k)^j=0$

证明
$\sum_{j=0}^{n-1}(w_n^k)^j=0=\frac{1-( w_n^k )^n}{1-w_n^k}=\frac{1-w_n^{kn}}{1-w_n^k}=\frac{1-1^k}{1-w_n^k}=0$

注意这里$k$不能是$n$的倍数，不然分母就为$0$了。

---

计算点值表达

如何快速地计算点值表达呢？

朴素地计算点值显然是$O(n^2)$的。
但是我们可以巧妙地选取复数根作自变量的值，加速运算。

核心思想

我们选取$n$次复数根作为自变量的值来求点值表达。
那么计算$A(w_n^k)=c_0+c_1w_n^k+\cdots+c_{n-1}(w_n^k)^{n-1}$
令$A_0(x)=c_0+c_2x^2+\cdots+c_{2k}x^k$
$A_1(x)=c_1+c_3x+\cdots+c_{2k+1}x^k$
那么$A(w_n^k)=A_0(( w_n^k )^2)+w_n^kA_1(( w_n^k )^2)$
又根据折半引理，要求的自变量数减少了一半，但是总的多项式项数不变。

直到要求的自变量数只有一个，也就是$w_1^0=1$的时候，每个多项式的项数也只有一项常数项了，于是我们可以直接得到。根据上面的式子来递归计算就可以了。

具体实现

由于在实际实现程序的时候不能对整个数组进行复制，我们要考虑一下如何使用一些技巧来避免进行一整个数组的复制。

首先假如是递归地实现，大概是以下的流程。

• 如果多项式只剩下一项了，那么返回常数项。否则
• 将多项式分成奇数项组成的多项式和偶数项组成的多项式，以当前所有自变量的平方作为自变量递归计算两个多项式的值
• 对所有的自变量合并两个分多项式所得的结果

注意到在每一层的自变量都是$w_n^k$的形式，其中$k$取遍$0$到$n-1$。那么我们就没必要再开多余的空间来记录有哪些自变量。我们只需要简单地知道多项式的项数$n$就可以得到所有的自变量了。

于是递归版本大概长这个样子

• 如果多项式只剩下一项了，那么返回常数项。否则
• 将多项式分成奇数项组成的多项式和偶数项组成的多项式，以当前单位复数根的平方传递递归计算两个多项式的值
• 对复数根的幂对应合并两个分多项式所得的结果

但是这里还是需要多开空间来存储递归计算的结果，这一部分的处理也是这个算法比较巧妙的地方。

假如我们可以按照某种顺序来储存，并用某种顺序计算，满足
- 寻址快速$( O(1) )$
- 线性空间
- 无后效性

那么一个可行的算法就出来了。

一个简单而又有效的方式，是按照它编号的二进制翻转后的顺序进行储存系数和答案。（有趣的是，翻转再翻转就是本身。也就是说我们进行两次这样的运算以后我们得到的是原顺序，之后会提到这个性质）
当$n=8$的时候，大概是这样的一个顺序：

$0, 4, 2, 6, 1, 5, 3, 7$

而此时每一层存储的答案数组的意义是不同的。假如当前正在处理第$k$层$[ l, r=l+\frac{n}{2^{k-1}} )$这个区间。那么左半边的答案数组$f_{i}(l\leq i < mid)$代表的实际上是$[ l, mid=l+\frac{n}{2^k} )$区间这个多项式，以$w_{\frac{n}{2^k}}^{ i\mod{2^k} }$为自变量的值。

首先来看第一次递归。显然奇数项都在右边，偶数项都在左边。
第二次递归也是（编号整除二以后），每一次递归是将区间劈成两半，左边恰好是偶数项右边恰好是奇数项。这个结果是显然的，因为二进制翻转后最低位为$0$的编号总是排在左边，最低位为$1$的编号总是排在右边，相当于把当前这一层的奇数项和偶数项分开成为各自连续的一段。
于是分开多项式的问题解决了。

接下来就是转移的问题。因为$(w_{2n}^k)^2=(w_{2n}^{k+n})^2$，因此它们来自于上一层的转移是相同的。我们将所有的求值两两分成一组来考虑，它们是相互独立的。

$A( w_{2n}^k )=A0( w_n^k )+w_{2n}^kA1( w_n^k )$
$A( w_{2n}^{k+n} )=A0( w_n^k )+w_{2n}^{k+n}A1( w_n^k )=A0( w_n^k )-w_n^kA1(w_n^k )$

至此，层间的转移就基本上解决了。

大致的程序实现应该是这样的

• 将所有的编号翻转并建立起对应关系
• 递归求值，记当前层求值的区间为$[ l, r )$
• 如果项数为$1$，将该位的值记为常数项并返回
• 否则递归两侧，得到一半的多项式的对应复数根为自变量的值后，按照上面的式子得到这一层的值。具体来说
  
  1. 枚举当前转移的是哪一个分治区间。（ 若当前每个区间的长度为$2^k( k>0 )$，那么总共有$\frac{n}{2^k}$组，每一组内有$2^{ k-1 }$对转移 ）
  2. 枚举当前分治区间是哪一个转移，转移之。

高效率的实现方法

注意不要搞混枚举的顺序，假如我们先枚举了是第几个转移，然后再枚举区间，就会使得内存的读取十分的零散，使得效率大大降低。

至此，求点值的部分已经实现完毕了。

---

计算插值运算

理论部分

在实际操作的时候，点值运算和差值运算实际上相差不多。它相当于是通过将点值矩阵乘上范德蒙德矩阵的逆矩阵，然后得到系数矩阵。
而重点就在于得到范德蒙德矩阵的逆矩阵，不妨如下构造一个。

令$V$为$n\times n$的范德蒙德矩阵，我们令$V^{-1}_{i, j}=\frac{1}{ nV_{i, j} }$即可得到$V$的逆矩阵$V^{-1}$

证明
主要是要证明$V^{-1}V=E$

• 当$i\neq j$时，$( V^{-1}V )_{ i, j }=\sum_{ k=0 }^{ n-1 }V^{-1}_{i, k}V_{k, j}$
  又$V_{k, j}=( w_n^j )^k, V^{-1}_{ i, k }=\frac{1}{ nV_{i, k} }=\frac{1}{ n( w_n^k )^i }=\frac{1}{ n( w_n^i )^k }$
  因此原式$=\sum_{ k=0 }^{ n-1 } ( w_n^j )^k\cdot \frac{1}{ n( w_n^i )^k }=\frac{1}{n}\sum_{k=0}^{ n-1 } ( w_n^{ i-j } )^k$
  因为$w_n^{ i-j }\neq 0$，根据求和引理，原式$=0$
• 当$i=j$时，$原式=\frac{1}{n}\sum_{k=0}^{ n-1 } ( w_n^{i-j} )^k=\frac{1}{n}\sum_{k=0}^{ n-1 } 1^k=1$

故$V^{-1}V=E$，证毕。

实际操作

程序实现的时候，我们并不需要重新实现一段求插值的代码。因为求点值和求插值十分的相似，我们只需要修改一小部分的参数就可以了。具体来说

• 令$n$次单位复数根为$w_n^-1$
• 重新使用上述递归过程
• 将结果的每一位除$n$

就完成了。

---

快速傅里叶变换

上述的点值表达和插值运算综合在一起就称快速傅里叶变换(fast Fourier transfer)。
理论部分到此就结束了。
talk is cheap, show me your code.

---

延伸

• 模意义下的FFT
  
  • 模数为$k2^n+1$，且是质数
  • 模数是质数
  • 模数是一般数
  • 快速威尔士变换（fast Walsh transform），即FFT的布尔运算形式