快速傅立叶变换(FFT，Fast Fourier Transform )是信号处理的常见手段，其可将时域信号变换到频域信号，并且时域卷积运算对应于频域时变为简单乘法。 两个多项式的乘积，其系数的运算实际上也是卷积运算，所以可以使用FFT来计算多项式的乘法。 关于网上FFT算法的描述大多是在工程领域，这里从算法竞赛的角度分析如何简洁高效地实现FFT。
一，为什么要用FFT计算多项式的乘法是众所周知的。 如果使用普通的for循环直接计算系数，则算法的时间复杂度为o(n^2)级。 算法慢的理由是使用系数公式进行计算。 其实，另一种方法可以完全表达多项式。 那是点值公式。 通俗地说，只要给出n个不同的Xi并求出对应的Yi，这个n对(Xi，Yi )就可以准确地决定一个多项式的系数。 那么，如何根据点值公式计算多项式乘法的系数呢？ 首先，让我们计算n个a(Xi )和n个a(Xi )。 于是，a ) Xi ) b ) Xi )是多项式相乘后在Xi处的结果。 然后，有n个a ) Xi ) b ) Xi )，可以唯一地确定一个多项式。 但是，如果只是随意选择n个点，计算这n个点的值的时间复杂度仍然是o(n ) 2，没有任何改善。 但是，如果巧妙地选择这n个点，就可以加速这个计算过程，把计算这n个点的时间复杂度降低到o(nlogn )。 这n个点是与n次单位的多根对应的点。 第二，n次单位的多根性质FFT算法之所以能很快利用n次单位根的性质加速计算。 以下是使用的三个性质。 1 .消除引理的一个推理如下。 2 .折半引理的一个推论如下。 3 .引理以上三个引理相加。 这里不证明，读者很少能用n次单位根的定义来证明。 让我们主要谈谈这三个定理是如何加速计算n个点的值的。
三、DFT和FFT DFT的定义如下。 从DFT的定义可以看出，实际上要计算n个n次单位的多根中的y值。 进行DFT变换可以得到n个点的值。 另一方面，FFT是利用n次单位的多根性质来加速DFT计算的方法。
四、FFT的分而治之战略FFT采用分而治之战略，假设一个多项式分为两部分，下标为偶数部分的多项式a0(x )和下标为奇数部分的多项式a1 ) x )，则计算a ) x )具体如下： 根据折半引理，上述n个可能的值相差n/2个，所以可以递归地使次数为n/2
五、利用迭代法实现的FFT，应该假设我们已经解决了部分问题，用第三等式合并结果。 具体的合并方式用伪代码记述如下。 fork=0ton/2-1t=w * y0 [ k ] y [ k ]=y0 [ k ] ty [ k ] [ n/2 ]=y0
以上伪代码中for语句下的前三个步骤称为蝴蝶操作，其作用在此极为重要，可以用下图表示。 在上图中，每个交叉位置都是蝴蝶的操作，计算出的两个中间结果是下一个蝴蝶操作的基础。 从上图也可以看出，FFT的时间复杂度是o(nlogN )，因为只计算logn层，各层使用的w值都相同，总共只输入了n个。
仔细观察，第一次操作蝴蝶时使用的两个操作数不相邻，却输入x0，x1，xn-1。 怎样才能快速变成蝴蝶操作所需的排列呢？ 实际上，正如在FFT的分割统治战略中已经提到的那样，下标根据奇偶性分为两个子多项式。 那么，继续划分直到只有一个要素消失为止，把所有要素按顺序排列成一列，就是上图的顺序。 例如，序列0、1、2、3、4、5、6和7在第一次操作结束时为0、2、4、6、1、3、5和7，在第二次操作结束时为0、4、2、6、1、5和3 和上图一样。 这种置换称为位逆序置换(Bit Reversal )。 那么，如何快速实现这一点呢？ 看看这些数的二进制形式，看看用红色标记的1的变化，就知道以下规律。 如果当前j的最高一位是1，则给定该1，继续检查下一位是否为1，如果是，则继续为0，直到没有连续的1，使得正好没有1的位置为1； 否则，将最高的下一位设为1。 这样就可以得到生成以下j的算法。 代码如下。 for(intI=0，j=0； i n； I ) if(jI ) swap ) a[j]，a[j]； int k=n； while(j ) k=1) ) j=~k； j |=k； 最后，得到n个积分值后，如何恢复系数？ 答案是用逆矩阵求系数。 因为把这n个点的值写在矩阵中的话，计算如下。
可以证明，该矩阵的系数矩阵v是wwddx矩阵，其逆必然存在，各元素的逆等于下式。
这样，最终系数就可以顺利恢复了。
六、FFT算法具体实现FFT加速多项式乘法的具体过程如下。 (1)0) 2个以上
项式的最高次前面补0，得到两个2n次的多项式，设系数向量分别为v1,v2。 （2）求值：用FFT计算f1=DFT(v1),f2=DFT(v2)。这里得到的f1,f2就是上面所说的n个点值。 （3）乘法：把两个向量f1,f2的每一维对应相乘，得到向量f，它对应两个多项式乘积的点值表示。 （5）插值：用FFT计算v=IDFT(f)，其中v就是乘积的系数向量。 下面直接给出一份简洁而高效的FFT算法实现的代码和利用FFT加速多项式相乘的实现过程 const double PI = acos(-1.0);typedef complex<double> CD;// Cooley-Tukey的FFT算法，迭代实现。inverse = false时计算逆FFTinline void FFT(vector<CD> &a, bool inverse) {int n = a.size();for (int i = 0, j = 0; i < n; i++) // 原地快速bit reversal{if (j > i) swap(a[i], a[j]);int k = n;while (j & (k >>= 1)) j &= ~k;j |= k;}double pi = inverse ? -PI : PI;for (int step = 1; step < n; step <<= 1) { // 把每相邻两个“step点DFT”通过一系列蝴蝶操作合并为一个“2*step点DFT”double alpha = pi / step;// 为求高效，我们并不是依次执行各个完整的DFT合并，而是枚举下标k// 对于一个下标k，执行所有DFT合并中该下标对应的蝴蝶操作，即通过E[k]和O[k]计算X[k]// 蝴蝶操作参考：http://en.wikipedia.org/wiki/Butterfly_diagramfor (int k = 0; k < step; k++) { // 计算omega^k. 这个方法效率低，但如果用每次乘omega的方法递推会有精度问题。// 有更快更精确的递推方法，为了清晰起见这里略去CD omegak = exp(CD(0, alpha*k));for (int Ek = k; Ek < n; Ek += step << 1) { // Ek是某次DFT合并中E[k]在原始序列中的下标int Ok = Ek + step; // Ok是该DFT合并中O[k]在原始序列中的下标CD t = omegak * a[Ok]; // 蝴蝶操作：x1 * omega^ka[Ok] = a[Ek] - t; // 蝴蝶操作：y1 = x0 - ta[Ek] += t; // 蝴蝶操作：y0 = x0 + t}}}if (inverse)for (int i = 0; i < n; i++) a[i] /= n;}// 用FFT实现的快速多项式乘法inline vector<double> operator * (const vector<double>& v1, const vector<double>& v2) {int s1 = v1.size(), s2 = v2.size(), S = 2;while (S < s1 + s2) S <<= 1;vector<CD> a(S, 0), b(S, 0); // 把FFT的输入长度补成2的幂，不小于v1和v2的长度之和for (int i = 0; i < s1; i++) a[i] = v1[i];FFT(a, false);for (int i = 0; i < s2; i++) b[i] = v2[i];FFT(b, false);for (int i = 0; i < S; i++) a[i] *= b[i];FFT(a, true);vector<double> res(s1 + s2 - 1);for (int i = 0; i < s1 + s2 - 1; i++) res[i] = a[i].real(); // 虚部均为0return res;}