FFT多项式计算

傅里叶变换是一个很大的概念，里面包含太多太多东西，整篇在傅里叶变换中只是一个小应用而已。就我的理解，描述快速傅里叶变换优化多项式乘法的原理和具体实践。

傅里叶变换

傅里叶变换是一种线性积分变换，用于信号在时域（或空域）和频域之间的变换。

一般情况下，若傅里叶变换一词不加任何限定语，则指的是连续傅里叶变换。（连续）傅里叶变换定义如下：傅里叶变换将可积函数f:R→Cf: \mathbb{R} \rightarrow \mathbb{C}f:R→C表示成复指数函数的积分或级数性质。
f^(ξ)=∫−∞∞f(x)e−2πixξdx\hat{f} \left( \xi \right) = \int_{-\infty}^{\infty} f \left( x \right) e^{-2\pi i x \xi } dx f^(ξ)=∫−∞∞f(x)e−2πixξdx
其中ξ\xiξ为任意实数。自变量xxx表示时间，变换变量ξ\xiξ表示频率。在多数情况下，傅里叶变换是可逆的，即可以通过f^\hat{f}f^得到其原函数fff，定义如下：
f(ξ)=∫−∞∞f^(ξ)e2πiξxdξf \left( \xi \right) = \int_{-\infty}^{\infty} \hat{f} \left( \xi \right) e^{2\pi i \xi x} d \xi f(ξ)=∫−∞∞f^(ξ)e2πiξxdξ

离散傅里叶变换

为了在科学计算和数字信号处理等领域使用计算机进行傅里叶变换，必须将函数xnx_{n}xn定义在离散点而非连续域内，且须满足有限性和周期性等提交。这种情况下，使用离散傅里叶变换（Discrete Fourier Transform，DFT），DFT是傅里叶变换在时域和频域上都呈离散的形式，函数定义如下：
x^n=∑k=0N−1xke−i2πNknn=0,⋯,N−1\hat{x}_{n} = \sum_{k=0}^{N-1}x_{k} e^{-i \frac{2\pi}{N} k n} \qquad n = 0, \cdots, N-1 x^n=k=0∑N−1xke−iN2πknn=0,⋯,N−1
其逆变换为：
xk=1N∑n=0N−1x^ne−i2πNnkk=0,⋯,N−1x_{k} =\frac{1}{N} \sum_{n=0}^{N-1} \hat{x}_{n} e^{-i \frac{2\pi}{N} n k} \qquad k = 0, \cdots, N-1 xk=N1n=0∑N−1x^ne−iN2πnkk=0,⋯,N−1

单位根

数学上，nnn次单位根是nnn次幂为1的负数。定义为ωn=1(n=1,2,3,⋯)\omega^{n}=1 \quad \left( n=1,2,3,\cdots\right)ωn=1(n=1,2,3,⋯)，ω\omegaω为nnn次单位根。单位的nnn次根有nnn个：e2πki/n(k=0,1,2,⋯,n−1)e^{2\pi k i / n} \quad \left(k=0,1,2,\cdots, n-1 \right)e2πki/n(k=0,1,2,⋯,n−1)。我们用ωn0,ωn1,⋯,ωnn−1\omega_{n}^{0}, \omega_{n}^{1}, \cdots, \omega_{n}^{n-1}ωn0,ωn1,⋯,ωnn−1表示这个nnn个nnn次根，ωni=e2πki/n\omega_{n}^{i} = e^{2\pi k i / n}ωni=e2πki/n。

结合欧拉公式，可以很容易得看出，ωnk=e2πkni=e2πk−1ni⋅e2πni=ωnk−1⋅ωn\omega_{n}^{k} = e ^{2\pi \frac{k}{n} i } = e ^{2\pi \frac{k-1}{n} i } \cdot e ^{ \frac{2\pi}{n} i } = \omega_{n}^{k-1} \cdot \omega_{n}ωnk=e2πnki=e2πnk−1i⋅en2πi=ωnk−1⋅ωn

对于n≥0,k≥0,d≥0n \geq 0, k \geq 0, d \geq 0n≥0,k≥0,d≥0，有ωdndk=e2πdki/dn=e2πki/n=ωnk\omega_{dn}^{dk}= e^{2\pi dk i/dn} = e^{2\pi k i /n}= \omega_{n}^{k}ωdndk=e2πdki/dn=e2πki/n=ωnk。

对于∀n\forall n∀n满足n>0n > 0n>0且2∣n2 | n2∣n，有ωnk+n2=e2πin(k+n2)=e2πik/neπi=ωnk(cosπ+isinπ)=−ωnk\omega_{n}^{k + \frac{n}{2}} = e ^ {\frac{2 \pi i}{n} \left( k + \frac{n}{2}\right)} = e^{2 \pi i k / n} e^{\pi i} = \omega_{n}^{k} \left( cos \pi + i sin \pi \right) = -\omega_{n}^{k}ωnk+2n=en2πi(k+2n)=e2πik/neπi=ωnk(cosπ+isinπ)=−ωnk

多项式

形如a0+a1x+⋯+anxna_{0} + a_{1} x + \cdots + a_{n} x^{n}a0+a1x+⋯+anxn的多项式被称为多项式，记为A(x)=∑i=0n−1aixiA\left( x \right) = \sum_{i=0}^{n-1} a_{i} x^{i}A(x)=∑i=0n−1aixi，其中ai,i=0,1,⋯,na_{i}, i=0,1,\cdots,nai,i=0,1,⋯,n为系数，用系数表示多项式为a⃗=(a0,a1,⋯,an)\vec{a}=\left(a_{0}, a_{1}, \cdots, a_{n} \right)a=(a0,a1,⋯,an)，a⃗\vec{a}a就为多项式的系数表示法。

我们知道通过n+1n+1n+1个不同的点{(x0,y0),⋯,(xn,yn)}\left\{\left(x_{0},y_{0}\right),\cdots, \left(x_{n},y_{n}\right) \right\}{(x0,y0),⋯,(xn,yn)}就可以唯一确定一个nnn次多项式，{(x0,y0),⋯,(xn,yn)}\left\{\left(x_{0},y_{0}\right),\cdots, \left(x_{n},y_{n}\right) \right\}{(x0,y0),⋯,(xn,yn)}为多项式的点值表示法。

快速傅里叶变换

快速傅里叶变换（Fast Fourier Transform，FFT），是计算离散傅里叶变换及其逆变换的快算算法。按照DFT的定义，计算一个长为nnn的序列需要的时间复杂度为O(n2)O\left( n^{2} \right)O(n2)，而FFT通过把DFT矩阵分解为稀疏因子之积来快速计算变换，从而将时间复杂度降到O(nlogn)O\left( nlog n \right)O(nlogn)。

库里-图基快速傅里叶变换算法是最常见的FFT算法，基于分治策略，可将时间复杂度降为O(nlogn)O\left( nlog n \right)O(nlogn)。

假设一个多项式表示为A(x)=∑i=0n−1aixiA\left( x \right) = \sum_{i=0}^{n-1} a_{i} x^{i}A(x)=∑i=0n−1aixi，用系数表示法表示就是向量a⃗=(a0,a1,⋯,an)\vec{a}=\left(a_{0}, a_{1}, \cdots, a_{n} \right)a=(a0,a1,⋯,an)。用点值表示，将nnn个nnn次单位根ωn0,ωn1,⋯,ωnn−1\omega_{n}^{0}, \omega_{n}^{1}, \cdots, \omega_{n}^{n-1}ωn0,ωn1,⋯,ωnn−1带入到多项式为A(ωnk)=∑i=0n−1aiωnki,k=0,1,⋯,n−1A\left( \omega_{n}^{k} \right) = \sum_{i=0}^{n-1} a_{i} \omega_{n}^{ki}, k = 0,1,\cdots, n-1A(ωnk)=∑i=0n−1aiωnki,k=0,1,⋯,n−1。现在我们和DFT定义的公式比对，A(ωni)A\left(\omega_{n}^{i}\right)A(ωni)就是是aia_{i}ai的DFT，那么aia_{i}ai就是A(ωni)A\left(\omega_{n}^{i}\right)A(ωni)就是的IDFT。记为A(ωni)=DFT(ai)A\left(\omega_{n}^{i}\right) = DFT\left( a_{i} \right)A(ωni)=DFT(ai)。

下面我们就依据上面陈列的基础基础，对A(ωnk)A\left(\omega_{n}^{k}\right)A(ωnk)进行公式推导：
A(ωnk)=∑i=0n−1aiωnki=∑i=0n2−1a2iωn2ki+ωnk∑i=0n2−1a2i+1ωn2ki=∑i=0n2−1a2iωn2ki+ωnk∑i=0n2−1a2i+1ωn2ki\begin{matrix} A\left(\omega_{n}^{k}\right) &= \sum_{i=0}^{n-1} a_{i}\omega_{n}^{ki} \qquad \qquad \qquad \qquad \qquad \\ \\ & = \sum_{i=0}^{\frac{n}{2}-1} a_{2i}\omega_{n}^{2ki} + \omega_{n}^{k} \sum_{i=0}^{\frac{n}{2}-1} a_{2i+1}\omega_{n}^{2ki} \\ \\ & = \sum_{i=0}^{\frac{n}{2}-1} a_{2i}\omega_{\frac{n}{2}}^{ki} + \omega_{n}^{k} \sum_{i=0}^{\frac{n}{2}-1} a_{2i+1}\omega_{\frac{n}{2}}^{ki} \; \; \; \end{matrix} A(ωnk)=∑i=0n−1aiωnki=∑i=02n−1a2iωn2ki+ωnk∑i=02n−1a2i+1ωn2ki=∑i=02n−1a2iω2nki+ωnk∑i=02n−1a2i+1ω2nki
又由于
A(ωnk+n2)=∑i=0n−1aiωn(k+n2)i=∑i=0n2−1a2iωn2ki+ωnk+n2∑i=0n2−1a2i+1ωn2ki=∑i=0n2−1a2iωn2ki−ωnk∑i=0n2−1a2i+1ωn2ki\begin{matrix} A\left(\omega_{n}^{k + \frac{n}{2}}\right) &= \sum_{i=0}^{n-1} a_{i}\omega_{n}^{\left(k+ \frac{n}{2}\right)i}\qquad \qquad \qquad \qquad \quad \\ \\ & = \sum_{i=0}^{\frac{n}{2}-1} a_{2i}\omega_{n}^{2ki} + \omega_{n}^{k + \frac{n}{2}} \sum_{i=0}^{\frac{n}{2}-1} a_{2i+1}\omega_{n}^{2ki} \\ \\ & = \sum_{i=0}^{\frac{n}{2}-1} a_{2i}\omega_{\frac{n}{2}}^{ki} - \omega_{n}^{k} \sum_{i=0}^{\frac{n}{2}-1} a_{2i+1}\omega_{\frac{n}{2}}^{ki} \qquad \end{matrix} A(ωnk+2n)=∑i=0n−1aiωn(k+2n)i=∑i=02n−1a2iωn2ki+ωnk+2n∑i=02n−1a2i+1ωn2ki=∑i=02n−1a2iω2nki−ωnk∑i=02n−1a2i+1ω2nki
可以看出，对于k<n2k < \frac{n}{2}k<2n时，只要代入ωn2,ωn22,⋯,ωn2n2−1\omega_{\frac{n}{2}}, \omega_{\frac{n}{2}}^{2}, \cdots, \omega_{\frac{n}{2}}^{\frac{n}{2} - 1}ω2n,ω2n2,⋯,ω2n2n−1，就可以求出A(ωnk)A\left(\omega_{n}^{k}\right)A(ωnk)和A(ωnk+n2)A\left(\omega_{n}^{k + \frac{n}{2}}\right)A(ωnk+2n)。
A(ωnk)=A1(ωn2k)+ωnk⋅A2(ωn2k)A(ωnk+n2)=A1(ωn2k)−ωnk⋅A2(ωn2k)(公式1)\begin{matrix} A\left(\omega_{n}^{k}\right) \quad =& A_{1}\left(\omega_{\frac{n}{2}}^{k}\right) + \omega_{n}^{k} \cdot A_{2}\left(\omega_{\frac{n}{2}}^{k}\right) \\ \\ A\left(\omega_{n}^{k + \frac{n}{2}}\right) = & A_{1}\left(\omega_{\frac{n}{2}}^{k}\right) - \omega_{n}^{k} \cdot A_{2}\left(\omega_{\frac{n}{2}}^{k}\right) \end{matrix} (公式1) A(ωnk)=A(ωnk+2n)=A1(ω2nk)+ωnk⋅A2(ω2nk)A1(ω2nk)−ωnk⋅A2(ω2nk)(公式1)
不递归的伪代码如下所示：

void get_rev(int bit)
{for(int i=0; i<(1<<bit); i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}
/**
n为2的幂次方
op==1，DFT
op==-1，IDFT
*/
void fft(complex<double>  *a,int n,int op)
{for(int i=0; i<n; i++){if(i<rev[i]){complex<double> tmp = a[i];a[i] = a[rev[i]];a[rev[i]] = tmp;}}for(int step=1; step<n; step<<=1){complex<double>  wn=exp(complex<double> (0,op*PI/step));for(int j=0; j<n; j+=step<<1){complex<double> wnk(1,0);for(int k=j; k<j+step; k++){complex<double> x=a[k];complex<double> y=wnk*a[k+step];a[k]=x+y;a[k+step]=x-y;wnk*=wn;}}}if(op==-1){for(int i=0; i<n; i++)a[i]/=n;}
}

费了好大经历才把代码和原理结合在一起，在网上搜了好久都没有理解为何代码要这样写，后来在离散傅里叶变换和快速傅里叶变换课件中搜到一张图，瞬间明白，其中的子操作被称为蝴蝶操作。如下图所示，展示了一个8个点的DFT是怎样分解成4个2个点的DFT。

蝴蝶操作见下图，类似于蝴蝶形状，是不是和上述公式1的计算相符合。

FFT实现多项式乘法

当现在有多项式A(x)A\left(x\right)A(x)和多项式B(x)B\left(x\right)B(x)，令C(x)=A(x)⋅B(x)C\left(x\right) = A\left(x\right) \cdot B\left(x\right)C(x)=A(x)⋅B(x)，那么如何求CCC的系数cjc_{j}cj呢？

先引入一个定理，函数卷积的傅里叶变换是函数傅里叶变换的乘积。C(x)C\left(x\right)C(x)中次数为jjj的项由A(x)A \left(x\right)A(x)中次数为kkk的项和B(x)B\left(x\right)B(x)中次数j−kj-kj−k为的项相乘得到。那么我们先通过C(ωn)=A(ωn)B(ωn,)C\left(\omega_{n}\right) = A\left(\omega_{n}\right) B\left(\omega_{n},\right)C(ωn)=A(ωn)B(ωn,)，然后在通过IDFT计算cjc_{j}cj。

    fft(a,n,1);fft(b,n,1);for(int i=0; i<s; i++) a[i]*=b[i];fft(a,n,-1);

备注

欧拉公式eix=cosx+i⋅sinxe^{ix} = cos x + i \cdot sin xeix=cosx+i⋅sinx

参考

傅里叶变换
离散傅里叶变换
单位根
快速傅里叶变换
离散傅里叶变换和快速傅里叶变换

利用快速傅里叶计算多项式相乘相关推荐

多项式相乘快速算法原理及相应C代码实现
最近认真研究了一下算法导论里面的多项式乘法的快速计算问题,主要是用到了FFT,自己也实现了一下,总结如下. 1.多项式乘法两个多项式相乘即为多项式乘法,例如:3*x^7+4*x^5+1*x^2+5与 ...
算法学习三：使用霍纳规则计算多项式
霍纳规则中的算法思想在<算法导论>第二章的思考题中,描述了利用霍纳规则计算多项式的方法.以前自己在写程序的时候都是傻傻的简单粗暴地直接上了,看到这个算法的时候眼前一亮,就多看了一些,果然 ...
【PAT甲级多项式相乘】1009 Product of Polynomials (25 分) C++ 全部AC
题目思路维护三个数组: arrA[1001]存储第一行数据 arrB[1001]存储第二行数据 c[1000000]存储计算结果数组下标表示多项式的指数,数组存的内容表示多项式的系数将arrA ...
信息学奥赛一本通 1012：计算多项式的值 | OpenJudge NOI 1.3 07
[题目链接] ybt 1012:计算多项式的值 OpenJudge NOI 1.3 07:计算多项式的值 [题目考点] 1. 计算表达式书写了解*的运算优先级比+高. 了解()可以改变运算优先级 2 ...
利用FFT计算非平稳随机信号的WVD分布
up目录一.理论基础二.核心程序三.测试结果一.理论基础 fft: 快速傅里叶变换 (fast Fourier transform), 即利用计算机计算离散傅里叶变换(DFT)的高效.快速计算 ...
python利用近似公式计算π_python如何利用公式计算π
python利用公式计算π的方法:首先导入数学模块及时间模块:然后计算Pi精确到小数点后几位数,代码为[print('n{:=^70}'.format('计算开始'))]:最后完成计算,代码为[pri ...
BUAA(2021春)多项式相乘
BUAA数据结构第三次编程题--多项式相乘看前须知第三次上机题汇总题目内容问题描述输入形式输出形式样例样例说明题解易错点和难点参考代码补充测试的数据看前须知要点介绍和简要声 ...
霍纳法则——计算多项式的值
例:,计算x=7时p(x)的值. 平时计算我们都会想到直接将x=代入方程中直接求解.但是这样子的话计算量很大,效率不高.而使用霍纳法则,则可以提高计算的效率. 一.霍纳法则的步骤 1.首先建立一个二维 ...
python计算csv文件内的数据_Python利用pandas计算多个CSV文件数据值的实例
功能:扫描当前目录下所有CSV文件并对其中文件进行统计,输出统计值到CSV文件 pip install pandas import pandas as pd import glob,os,sys in ...
Matlab计算多项式的值（数值）
MATLAB 中,多项式用一个行向量表示,行向量的元素值为多项式系数按幂次的降序排列: 例如多项式, P(x) = 2*x^4 + 3*x^3 - 2*x^2 + 7*x + 11 可表示为, p = ...

利用快速傅里叶计算多项式相乘