多项式与快速傅里叶变换(FFT)

计算两个n次多项式相乘的最直接方法所需的时间为$\Theta (n^2)$，在本文中，我们将讨论快速傅里叶变换(FFT)，使得多项式相乘的时间复杂度降为$\Theta (n\text{lg}n)$。

一、多项式的表示

系数表示和点值表示

对于一个次数界为n的多项式$A(x)={\sum }_{j=0}^{n-1}{a}_j{x}^j$而言，其系数表示是一个由系数组成的向量$a=(a_0,a_1,...,a_{n-1})$。而$A(x)$的点值表示就是一个由n个点组成的集合$\{(x_0,y_0),(x_1,y_1),...(x_{n-1},y_{n-1})\}$，使得对$k=0,1,...,n-1$，所有$x_k$各不相同。

对于一个用系数表示的多项式来说，只要选取n个不同的点$x_0,x_1,...,x_{n-1}$，然后在$\Theta (n^2)$的时间内分别计算出各点对应$A(x)$的值，就可以根据这些值解出多项式的系数，这一过程（从一个多项式的点值确定表达式的系数）称为插值。值得注意的是，如果我们可以巧妙地选取$n$个点，那么计算$n$个点对应值的过程所需时间可以降为$\Theta (n\lg n)$。下面定理说明，当多项式的次数界等于已知的点值对的数目，插值才是正确的。

（插值多项式的唯一性） 对于任意n个点值对组成的集合$\{(x_0,y_0),(x_1,y_1),...(x_{n-1},y_{n-1})\}$，其中所有的$x_k$互不相同，那么存在唯一的一个次数界为n的多项式$A(x)$，满足$y_k=A(x_k)$，k=0,1,…,n-1。

证明 \quad设$A(x)$的系数向量为$(a_0,a_1,...,a_{n-1})$，将方程组(1)化为矩阵形式(2)：
$$\begin{array}{c}\{\begin{array}{ll}A(x_0)=a_0+a_1{x}_0+...+a_{n-1}{x}_0=y_0& \\ A(x_1)=a_0+a_1{x}_1+...+a_{n-1}{x}_1=y_1& \\ \cdots & \\ A(x_{n-1})=a_0+a_1{x}_{n-1}+...+a_{n-1}{x}_{n-1}=y_{n-1}& \end{array}\end{array}$$

$$\begin{array}{c}\left[\begin{array}{ccccc}1& x_0& x_0^2& \cdots & x_0^{n-1}\\ 1& x_1& x_1^2& \cdots & x_1^{n-1}\\ \cdots & & & & \\ 1& x_{n-1}& x_{n-1}^2& \cdots & x_{n-1}^{n-1}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ \cdots \\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ \cdots \\ y_{n-1}\end{array}\right]\end{array}$$

左边的矩阵记为$V(x_0,x_1,...,x_{n-1})$，根据范德蒙德矩阵的性质，矩阵$V$行列式的值$|V|=\prod _{0\le j<k\le n-1}(x_k-x_j)\ne 0$，因此，$V$是可逆的，故存在唯一的系数矩阵$a=V^{-1}y$。

对于多项式乘法，点值表达是方便的，设$C(x)=A(x)B(x)$，则对于任意点$x_k$，$C(x_k)=A(x_k)B(x_k)$，每次操作的时间为$\Theta (1)$。由于degree($C$)=degree($A$)+degree($B$)，当A和B的次数界都为n时，C的次数界为2n，因此需要对A和B的点值表达式进行扩展，使每个多项式都含2n个点值对。给定A的扩展点值表达$\{(x_0,y_0),(x_1,y_1),...(x_{n-1},y_{n-1})\}$和B的扩展点值表达$\{(x_0,y_0^{\prime }),(x_1,y_1^{\prime }),...(x_{n-1},y_{n-1}^{\prime })\}$，则C的点值表达为$\{(x_0,y_0{y}_0^{\prime }),(x_1,y_1{y}_1^{\prime }),...(x_{n-1},y_{n-1}{y}_{n-1}^{\prime })\}$。给定两个点值扩展形式的输入多项式，得到其相乘结果的点值形式所需时间为$\Theta (n)$。

---

系数形式表示多项式的快速乘法

我们需要利用基于点值形式表达的多项式线性时间乘法算法，来加速基于系数形式表达的多项式乘法运算。关键在于快速将一个多项式从系数形式转换为点值形式(求值)以及从点值形式转换为系数形式(插值)。

通过巧妙地选择求值点，我们可以把两种表示方法之间的转换所需的时间复杂度变为$\Theta (nlgn)$，在第二节中选择单位复数根作为求值点，对系数向量做离散傅里叶变换(DFT)，得到对应的点值表达，同样我们可以通过对这些点值执行逆DFT变换，获得其相应的系数向量，这样就实现了插值运算。

如下图所示，两个次数界为n的多项式的乘积是一个次数界为2n的多项式，因此，在输入多项式A和B进行求值前，先把这两个多项式在最高次系数前添加n个0，使其次数界增大为2n。

给定FFT，我们就有下面时间复杂度为$\Theta (n\lg n)$的算法，把两个次数界为n的多项式$A(x)$和$B(x)$进行乘法运算，其中输入和输出均以系数表达，假设n是2的幂(实际可以通过向高阶系数补0的方式来满足这个要求)：

1. 加倍次数界：通过加入n个系数为0的高阶系数，把多项式$A(x)$和$B(x)$变为次数界为2n的多项式，并构造其系数表达。
2. 求值：通过应用2n阶的FFT计算出$A(x)$和$B(x)$的长度为2n的点值表达，这些点值表达包含了两个多项式在2n次单位根处的取值。
3. 逐点相乘：把$A(x)$的值和$B(x)$的值逐点相乘，可以计算多项式$C(x)=A(x)B(x)$的点值表达，这个点值表达包含了$C(x)$在每个2n次单位根处的取值。
4. 插值：通过对2n个点值对应用FFT，计算其逆DFT，就可以构造出多项式$C(x)$的系数表达。

其中，执行第一步和第三步的时间为$\Theta (n)$，执行第二步和第四步的时间为$\Theta (n\lg n)$。

二、DFT与FFT

单位复数根

n次单位复数根是满足$w^n=1$的复数$w$。n次单位复数根恰好有n个：对于$k=0,1,\cdots ,n-1$，这些根是$e^{2\pi ik/n}$。这里我们通过欧拉定理用指数形式表示复数：$e^{iu}=\cos (u)+i\sin (u)$，下图说明$n$个单位复数根均匀地分布在以复平面的原点为圆心的单位圆周上，$w_n=e^{2\pi i/n}$称为主$n$次单位根，其他所有n次单位复数根都是$w_n$的幂次。

n个n次单位复数根$w_n^0,w_n^1,\cdots ,w_n^{n-1}$在乘法意义下形成一个群，该群与加法群(整数模n)具有同样的结构，即$w_n^n=w_n^0$意味着$w_n^j{w}_n^k=w_n^{(j+k)\%n}$。对于n次单位复数根有以下性质：

（消去引理） 1. 对于任何整数$n\ge 0,k\ge 0$以及$d>0$，有$w_{dn}^{dk}=w_n^k$。

（折半引理） 2. 如果$n>0$为偶数，那么n个n次单位复数根的平方的集合就是$n/2$个$n/2$次单位复数根的集合。
（求和引理） 3. 对任意整数$n\ge 1$和不能被$n$整除的非负整数$k$，有$$\sum _{j=0}^{n-1}(w_n^k{)}^j=0$$

---

DFT

我们希望计算次数界为n的多项式$A(x)={\sum }_{j=0}^{n-1}{a}_j{x}^j$在$n$个单位复数根$w_n^0,w_n^1,\cdots ,w_n^{n-1}$处的值。假设$A$以系数形式给出：$a=(a_0,a_1,\cdots ,a_{n-1})$，接下来对$k=0,1,\cdots ,n-1$，定义结果$y_k$:$$\begin{array}{c}{y}_k=A(w_n^k)=\sum _{j=0}^{n-1}{a}_j{w}_n^{kj}\end{array}$$向量$y=(y_0,y_1,\cdots ,y_{n-1})$就是系数向量$a=(a_0,a_1,\cdots ,a_{n-1})$的离散傅里叶变换(DFT)，记为$y=DF{T}_n(a)$。

---

FFT

通过使用一种称为快速傅里叶变换(FFT)的方法，利用单位复数根的性质，我们就可以在$\Theta (nlgn)$的时间内计算出$DF{T}_n(a)$，这里假设$n$恰好是2的整数幂。

FFT采用了分治策略，采用$A(x)$中偶数下标的系数和奇数下标的系数，分别定义两个新的次数界为$n/2$的多项式$A^{[0]}(x)$和$A^{[1]}(x)$：
$$\begin{gathered} A^{[0]}(x)=a_0+a_{2}x+a_4{x}^2+\cdots +a_{n-2}{x}^{n/2-1} \\ {A}^{[1]}(x)=a_1+a_{3}x+a_5{x}^2+\cdots +a_{n-1}{x}^{n/2-1} \end{gathered}$$
注意到，$A^{[0]}(x)$包含$A$中所有偶数下标的系数，$A^{[1]}(x)$包含$A$中所有奇数下标的系数，于是有
$$\begin{array}{c}A(x)=A^{[0]}(x^2)+x{A}^{[1]}(x^2)\end{array}$$
于是原问题转化为：

1. 求次数界为$n/2$的多项式$A^{[0]}(x)$和$A^{[1]}(x)$在点$(w_n^0{)}^2,(w_n^1{)}^2,\cdots ,(w_n^{n-1}{)}^2$处的取值
2. 根据式(4)综合上述结果得到$A(x)$的值，即$DF{T}_n$

根据折半引理，$(w_n^0{)}^2,(w_n^1{)}^2,\cdots ,(w_n^{n-1}{)}^2$可以看作仅由n/2个n/2次单位复数根所组成，每个根正好出现两次，因此，可以递归地对次数界为n/2的多项式$A^{[0]}(x)$和$A^{[1]}(x)$在$n/2$个$n/2$次单位复数根处进行求值。这些子问题与原问题形式相同，但规模减小为原问题的一半，即将一个n个元素的$DF{T}_n$划分为两个规模为$n/2$个元素的$DF{T}_{n/2}$。这一分解是下面递归FFT算法的基础，此算法计算出一个由$n$个元素组成的向量$a=(a_0,a_1,\cdots ,a_{n-1})$的$DF{T}_n$，其中$n$是2的幂。

RECURSIVE-FFT的执行过程如下：

第2-3行是递归的基础，一个元素的DFT就是它本身，因为在这种情形下，$y_0=a_0{w}_1^0=a_0$。

第6-7行定义多项式$A^{[0]}(x)$和$A^{[1]}(x)$的系数向量，第4、5、13行保证$w$可以正确更新，在执行第11-12行时，始终有$w=w_n^k$。第8-9行递归执行$DF{T}_{n/2}$，对于$k=0,1,\cdots ,n-1$，
$$\begin{gathered} y_k^{[0]}=A^{[0]}(w_{n/2}^k)=A^{[0]}(w_n^{2k}) \\ {y}_k^{[1]}=A^{[1]}(w_{n/2}^k)=A^{[1]}(w_n^{2k}) \end{gathered}$$
第11-12行综合了递归$DF{T}_{n/2}$的结果，对$y_0,y_1,\cdots ,y_{n/2-1}$，第11行推出：
$$\begin{array}{rl}{y}_k& =y_k^{[0]}+w_n^k{y}_k^{[1]}\\ & =A^{[0]}(w_n^{2k})+w_n^k{A}^{[1]}(w_n^{2k})\\ & =A(w_n^k)式(4)\end{array}$$
对$y_{n/2},y_{n/2+1},\cdots ,y_{n-1}$，设$k=0,1,\cdots ,n/2-1$，第12行推出：
$$\begin{array}{rl}{y}_{k+n/2}& =y_k^{[0]}-w_n^k{y}_k^{[1]}\\ & =A^{[0]}(w_n^{2k})+w_n^{k+n/2}{A}^{[1]}(w_n^{2k})折半引理\\ & =A^{[0]}(w_n^{2k+n})+w_n^{k+n/2}{A}^{[1]}(w_n^{2k+n})\\ & =A(w_n^{k+n/2})式(4)\end{array}$$
至此，RECURSIVE-FFT的正确性得以证明，它可以返回输入向量$a$的DFT。

关于RECURISVE-FFT的执行时间，我们可以通过该递归式得到：$T(n)=2T(n/2)+\Theta (n)=\Theta (nlgn)$。这里的$n$是输入向量$a$的长度，每次除递归调用外，所需时间为算法第10-13行的for循环所消耗$\Theta (n)$，因此上述递归式成立。

---

在单位复数根处插值

接下来我们需要在单位复数根处插值将多项式由点值表达转化为系数表达，首先将DFT写成矩阵方程，然后观察它逆矩阵的形式。

根据式(2)，我们可以把DFT写成矩阵乘积$y=V_{n}a$，其中$V_n$是一个由$w_n$适当幂次填充成的范德蒙德矩阵：
$$\begin{array}{c}\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ \cdots \\ y_{n-1}\end{array}\right]=\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& w_n^1& w_n^2& \cdots & w_n^{n-1}\\ 1& w_n^2& w_n^4& \cdots & w_n^{2(n-1)}\\ \cdots & & & & \\ 1& w_n^{n-1}& w_n^{2(n-1)}& \cdots & w_n^{(n-1)(n-1)}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ \cdots \\ a_{n-1}\end{array}\right]\end{array}$$
对$j,k=0,1,\cdots ,n-1$，$V_n$的$(k,j)$处元素为$w_n^{kj}$。对于逆运算$DF{T}^{-1}$，我们把$y$乘以$V_n$的逆矩阵$V_n^{-1}$来处理得到系数向量$a$。

对$j,k=0,1,\cdots ,n-1$，$V_n^{-1}$的$(j,k)$处元素为$w_n^{-kj}/n$。

证明 \quad 已知$V_n^{-1}{V}_n=I_n$。考虑矩阵乘法结果在$(j,j^{\prime })$处的元素：
$$[V_n^{-1}{V}_n{]}_{j{j}^{\prime }}=\sum _{k=0}^{n-1}(w_n^{-kj}/n)(w_n^{k{j}^{\prime }})=\sum _{k=0}^{n-1}(w_n^{k(j^{\prime }-j)}/n)$$
如果$j^{\prime }=j$，则结果为${\sum }_{k=0}^{n-1}\frac{1}{n}=1$；否则，根据求和引理($j^{\prime }-j$不能被$n$整除)，结果为0。

因此，给定逆矩阵$V_n^{-1}$，可以推导出$DF{T}_n^{-1}(y)$：
$$\begin{array}{c}{a}_j=\frac{1}{n}\sum _{k=0}^{n-1}{y}_k{w}_n^{-kj}\end{array}$$
将上式(6)与式(3)：$y_k=A(w_n^k)={\sum }_{j=0}^{n-1}{a}_j{w}_n^{kj}$对比，可以得到只需要对FFT算法做出如下修改就可以在$\Theta (nlgn)$的时间内计算出$DF{T}_n^{-1}$：将$a$与$y$互换，用$w_n^{-1}$替换$w_n$，并将结果的每个元素除以$n$。

---

快速傅里叶变换的代码实现

根据上述分析，求值和插值的过程可以共用一份代码，唯一不同的是$w_n$和结果向量是否需要除以长度。在求值过程，$w_n=e^{2\pi i/n}=\cos (\frac{2\pi }{n})+\sin (\frac{2\pi }{n})i$；在求值过程中，$w_n=e^{-2\pi i/n}=\cos (-\frac{2\pi }{n})+\sin (-\frac{2\pi }{n})i=\cos (\frac{2\pi }{n})-\sin (\frac{2\pi }{n})i$。可以发现，二者仅在虚部符号相反，因此可以在算法中加入一个参数$op$乘以$w_n$虚部的值：当$op=1$，表示进行FFT；当$op=-1$，表示进行IFFT。

下面是快速傅里叶变换算法的代码：

void FFT(plural A[], int n, int op)
{
	if (n == 1)  return;
	plural A1[n] = { 0 };
	plural A2[n] = { 0 };
	for (int i = 0; i < n / 2; i++) {
		A1[i] = A[i * 2];
		A2[i] = A[i * 2 + 1];
	}
	FFT(A1, n / 2, op);
	FFT(A2, n / 2, op);
	plural w, w_n;
	w_n.real = cos(2.0 * Pi / n);
	w_n.img = sin(2.0 * Pi / n) * op;
	w.real = 1;    w.img = 0;
	for (int i = 0; i < n / 2; i++)
	{
		A[i] = add(A1[i], mul(A2[i], w));
		A[i + n / 2] = sub(A1[i], mul(A2[i], w));
		w = mul(w, w_n); 
	}
}

调用方式如下：

	memset(A, 0, sizeof(A));
	memset(B, 0, sizeof(B));
	// 输入 A、B的系数向量
	input(A, B);
    FFT(A, n, 1);   // 计算A在单位复数根处的值
	FFT(B, n, 1);	// 计算B在单位复数根处的值
	// 以点值形式表示的多项式相乘,结果存储在A中
	for (int i = 0; i < n; i++)
		A[i] = mul(A[i], B[i]);
	// 插值操作，由点值形式转化为系数形式
	FFT(A, n, -1);
	for(int i = 0; i < n; i++)
		 printf("%d ",(int)(A[i].real/n + 0.5));

三、高效FFT的实现

FFT的一种迭代实现

首先我们注意到，在RECURSIVE-FFT中，第10-13行的for循环中包含了$w_n^k{y}_k^{[1]}$的两次计算，我们称该值为公用子表达式，我们可以改变循环，使其仅计算一次，并将其存储到变量中。

$$\begin{gathered} \text{for}k=0\text{to}n/2-1 \\ t=w{y}_k^{[1]} \\ {y}_k=y_k^{[0]}+t \\ {y}_{k+(n/2)}=y_k^{[0]}-t \\ w=w{w}_n \end{gathered}$$

在这个循环中，我们将循环因子$w_n$乘以$y_k^{[1]}$，把所得结果存入变量$t$中，然后从$y_k^{[0]}$增加或减去$t$，这一系列操作被称为一个蝴蝶操作(butterfly operation)，下图说明了操作执行步骤：

接下来将算法的递归结构修改为迭代结构。下图是输入一个长度n=8的向量RECURSIVE-FFT执行过程中的递归结构(一棵向量树)，树中的每一个节点对应每次过程的递归调用，由相应的输入向量标记。每次RECURSIVE-FFT调用产生两个递归调用，除非该调用已收到了1个元素的向量，第一次调用作为左孩子，第二次调用作为右孩子。

只要把初始向量$a$中的元素按照其在叶中出现次序进行安排，就可以对过程RECURSIVE-FFT的执行进行追踪，不过是自底向上而不是自顶向下。首先，我们成对地提取出元素，利用一次蝴蝶操作计算出每对的DFT，然后对其DFT取代这对元素，这样向量中就包含了$n/2$个二元素的DFT。下一步，按对取出这$n/2$个DFT，通过两次DFT蝴蝶操作计算出具有$n/4$个四元素的DFT，并用一个具有四个元素的DFT取代对应的两个二元素的DFT。继续这样的操作，直到向量包含两个具有$n/2$个元素的DFT，这时，我们可以综合$n/2$次蝴蝶操作，合成具有$n$个元素的DFT了。

我们用一个数组A[0 … n-1]，初始时该数组包含输入向量$a$中的元素，其顺序为上图递归树中树叶顺序(得到这一序列的过程称为位逆序置换)。因为需要在树的每一层进行组合，于是引入变量$s$来计算树的层次，其取值范围为从$1$到$\lg n$，因此，该算法有以下结构：
$$\begin{gathered} 1\text{for}s=1\text{to}\lg n \\ 2\text{for}k=0\text{to}n-1\text{by}{2}^s \\ 3\text{combine the two}{2}^{s-1}\text{-element DFTs in}A[k\cdots k+2^{s-1}-1]\text{and}A[k+2^{s-1}\cdots k+2^s-1] \\ \text{ into one}{2}^s\text{-element DFT in }A[k\cdots k+2^s-1] \end{gathered}$$

可以通过如下方式实现上述算法：从子程序RECURSIVE-FFT中复制for循环，让$y^{[0]}$与$A[k\cdots k+2^{s-1}-1]$一致，让$y^{[1]}$与$A[k+2^{s-1}\cdots k+2^s-1]$一致。在每次蝴蝶操作中，使用的旋转因子依赖于$s$的值，它是$w_{2^s}$的幂，我们引入临时变量$u$使得它能恰当地执行蝴蝶操作(数据不会被覆盖)。用循环主体代替上述算法第三行内容，就可以得到下面的伪代码，代码首先调用辅助过程BIT-REVERSE，把向量$a$按我们所需要的初始顺序复制到数组$A$中。

接下来考虑BIT-REVERSE是如何将输入向量$a$按希望顺序放入数组$A$中。在递归树中，叶出现顺序是一个位逆序置换，也就是说，如果让$rev(k)$为$k$的二进制表示各位逆序所形成的$\lg n$位数，那么我们希望把向量中的元素$a_k$放在数组的$A[rev(k)]$位置上。如下表所示，当$k=3$，$rev(k)=110$，也就是说$a_3$出现在$A[6]$处。

index	0	1	2	3	4	5	6	7
$a$	$a_0$	$a_1$	$a_2$	$a_3$	$a_4$	$a_5$	$a_6$	$a_7$
二进制	000	001	010	011	100	101	110	111
$A$	$a_0$	$a_4$	$a_2$	$a_6$	$a_1$	$a_5$	$a_3$	$a_7$
二进制rev(index)	000	100	010	110	001	101	011	111

注意到在树的最顶层，最低位为0的下标在左子树中，最低位为1的下标在右子树中。在每一层去掉最低位后，我们沿着树往下继续这一过程，直到在叶子处得到由位逆序置换给出的顺序。于是BIT-REVERSE过程如下：
$$\begin{gathered} \text{BIT-REVERSE}(a,A) \\ 1n=a.length \\ 2\text{for}k=0\text{to}n-1 \\ 3A[rev(k)]=a_k \end{gathered}$$

这种迭代方式实现的时间位$\Theta (n\lg n)$，这是由于迭代了$n$次，每次迭代可以在$O(\lg n)$的时间内把一个介于0和n-1的$\lg n$位数逆序。

在实际中通常可以通过预处理，在$\Theta (n)$的时间内完成位逆序置换。在$a$中偶数下标处的位逆序置换结果二进制首位为0，奇数位下标处的位逆序置换结果二进制首位为1；如果按照0-1，2-3的方式分组，每两个相邻的二进制数除首位外其余位都相同，并且除首位外是上一级子问题的解；$A$中前一半($A[0]$到$A[n/2-1]$)二进制数除最后一位外是上一级子问题的解。综上，$rev[k]$有如下关系：
$$rev[k]=(rev[k/2]>>1)|((k\&1)\times \frac{n}{2})$$

---

并行FFT电路

我们可以利用迭代FFT算法的性质来产生一个高效的并行FFT算法。如下图所示，输入长度为$n=8$的向量$a$，用一个并行FFT电路计算FFT。该电路开始时对输入进行位逆序置换，其后电路分$\lg n$级，每一级由$n/2$个并行执行的蝴蝶操作组成。电路的深度定义为任意的输入和任意的输出之间最大的可以达到的计算元素的数目，因此，并行FFT电路的深度为$\Theta (\lg n)$。

并行FFT电路的最外层执行位逆序操作，其余部分模拟迭代的ITERATIVE-FFT过程。因为最外层for循环的每次迭代执行$n/2$次独立的蝴蝶操作，于是电路并行地执行它们。在ITERATIVE-FFT内每次迭代的值$s$对应于上图中一个阶段的蝴蝶操作，对于$s=1,2,\cdots ,\lg n$，阶段$s$有$n/2^s$组蝴蝶操作，每组中有$2^{s-1}$个蝴蝶操作。蝴蝶操作用到的旋转因子对应于ITERATIVE-FFT中所用的旋转因子：在阶段$s$，我们使用$w_m^0,w_m^1,\cdots ,w_m^{m/2-1}$，其中$m=2^s$。