利用快速傅里叶变换(FFT)计算多项式乘法

利用快速傅里叶变换(FFT)计算多项式乘法作者:宋振华摘要本文将讨论快速傅里叶变换(FFT),利用FFT设计一种算法,使多项式相乘的时间复杂度降低为nn2log,以便在计算机上高效计算多项式乘法.关键词:快速傅里叶变换、多项式乘法目录一、引言二、算法概述三、引理1、多项式的表示(1)系数表达形式(2)点值表达形式(3)插值2、利用离散傅里叶变换(DFT)与FFT导出结果的点值表达形式(1)单位复数根(2)离散傅里叶变换(DFT)(3)通过快速傅里叶变换(FFT)计算向量y3、利用FFT计算逆DFT,将结果的点值表达形式化为系数表达形式(1)在单位复数根处插值(2)利用FFT计算逆DFT四、算法具体流程五、算法的实际应用:计算大整数乘法六、参考文献一、引言基于FFT的离散傅里叶变换(DFT)技术,是当今信息传输、频谱分析等领域中最重要的数学工具之一.在计算机编程中,我们经常需要计算两个多项式函数的乘积.对于两个n次多项式函数,计算其乘积最直接方法所需时间为2n.本文将讨论快速傅里叶变换(FFT),利用FFT设计一种算法,使多项式相乘的时间复杂度降低为nn2log,以便在计算机上高效执行.二、算法概述通过FFT进行离散傅里叶变换,将两个多项式由系数表达形式转化为点值表达形式.将这2个点值表达形式的多项式相乘,得到结果的点值表达形式.最后利用FFT做DFT的逆,得到结果的系数表达形式.三、引理1、多项式的表示(1)系数表达对于次数为n的多项式niiixaxA0,其系数表达是一个由系数组成的向量Tnaaaa,...,,10.考虑用系数形式表示、次数为n的多项式)(xA、)(xB的乘法运算.记niiixcxBxAxC20)()()(.有ijjijibac0.可以看出,采取逐个计算ic),20(Nini的方式进行求解,其时间复杂度为)(2n.(2)点值表达a.定义:一个次数为n的多项式)(xA的点值表达就是由一个有n个点值对组成的集合)}(,,0|),{(iiiixAyNiniyx,使得对于nk,...,2,1,0,所有的kx各不相同.b.点值表达下多项式乘法对于n次多项式)(xA,)(xB,设其点值表达式分别为:)(xA:)},(,),,(),,(),,{(221100nnyxyxyxyx,)(xB:)},(,),,(),,(),,{(,,2211,00nnyxyxyxyx，设)()()(xBxAxC,由于)(xC的次数为n2,因此必须对)(xA和)(xB的点值表达式进行扩展,使每个多项式都包含12n个点值对.给定A,B的扩展点值表达:)(xA:)},(,),,(),,(),,{(22221100nnyxyxyxyx,)(xB:)},(,),,(),,(),,{(,22,2211,00nnyxyxyxyx，.则)(xC的点值表达形式为:)},(,),,(),,(),,{(,222,222111,000nnnyyxyyxyyxyyx，.因此,计算)(xC点值表达式的时间复杂度为)(n.`(3)插值求值运算的逆(从一个多项式的点值表达形式确定其系数表达形式)称为插值.下面将证明,n个点求值运算与插值运算是定义完备的互逆运算.a.多项式的点值表达可以唯一确定多项式的系数表达形式.多项式的点值表达等价于矩阵方程:nnnnnnnxxxxxxxxxxxx22222121102001111naaaa210=nyyyy210.(3-1-3-a)记系数矩阵为V,由Vandermonde行列式可知,njijixxV0.而jinjiNji,,0,,,有jixx,因此0V,即V可逆.因此对于给定点值表达,我们能够唯一确定系数表达式,且yVa1.b.对于次数为n的多项式函数xA,插值算法基于如下Lagrange公式:nkkjjkkjjkxxxxyxA0)(.(3-1-3-b)容易验证,Lagrange公式的计算复杂度为)(2n.因此,n个点求值运算与插值运算是定义完备的互逆运算,它们将多项式的系数表达与点值表达进行相互转化.2、利用离散傅里叶变换(DFT)与FFT导出结果的点值表达形式(1)单位复数根n次单位复数根是满足1n的复数.n次单位复数根恰好有n个:对于1,...,1,0nk,这些根是nkie2.由Euler公式sincosiei可知,这n个单位复数根均匀地分布在以复平面原点为圆心的单位圆圆周上.nine2称为主n次单位根,所有其他n次单位根都是n的幂次.n个n次单位复数根0n,1n,2n,...,1nn在乘法意义下构成一个群,该群与群nnZ,同构.由此,可以得到如下推论:a.0,0,0dkn,有knnkidkdneedndki22.(3-2-1-a)b.*,2Nkkn,122nn.(3-2-1-b)c.若*,2Nkkn,},0|){(2Niniin=},20|){(22/Niniin.(3-2-1-c)对推论c的证明:由a可知,Nk,2/2/2knkn.所以对于20,nkNk,有222222/)()(knknnnknnknnkn.得证.d.NkNn,*且k不能被n整除,有100ninkn.(3-2-1-d)对推论d的证明:10ninkn11)(knnkn11)(knknn11)1(knk0.(2)离散傅里叶变换(DFT)对于n次多项式niiixaxA0)(,其系数形式为Tnaaaa),,,(10.(3-2-2-1)设nikiniknkaAy01)(,Nknk,0,(3-2-2-2)则向量Tnyyyy),,,(10(3-2-2-3)就是系数向量Tnaaaa),,,(10的离散傅里叶变换.(3)通过快速傅里叶变换(FFT)计算向量y.对于n次多项式niiixaxA0)(,不妨假设Npnp,21.对于1n不为2的整数次幂的情况类似,此处不再讨论.FFT采取了分治策略,采用)(xA中偶数下标与奇数下标的系数,分别定义两个新的2n次多项式:21-124200)(nnxaxaxaaxA,(3-2-3-1)21-25311)(nnxaxaxaaxA.(3-2-3-2)注意到)(0xA包含A所有偶数下标的系数,)(1xA包含A中所有奇数下标的系数,于是有:)()()(2]1[20xxAxAxA.(3-2-3-3)所以,求)(xA在01n,11n,...,nn1处的值的问题转化为:a.求次数为2n的多项式)(]0[xA,)(]1[xA在点201)(n,211)(n,...,21)(nn处的取值.b.根据(2-2-3-3)综合结果.根据(2-2-1-c),式(2-2-3-3)不是由1n个不同值组成,而是仅由21n个21n次单位复数根所组成,每个根正好出现2次.因此,我们递归地对次数为21n的多项式)(]0[xA,)(]1[xA在21n个21n次单位复数根处进行求值.这些子问题与原始问题形式相同,但规模变为一半.下面确定DFT过程的时间复杂度.注意到除了递归调用外,每次调用需要枚举Tnaaaa),,,(10中所有元素,将a划分为},,{200naaaa、},,,,{15311naaaaa,分别与多项式)(]0[xA,)(]1[xA相对应.其时间复杂度为)(n.因此,对运行时间有下列递推式:)()2(2)(nnTnT.(3-2-3-4)求解该递推式,有)log()(2nnnT.(3-2-3-5)采取快速傅里叶变换,我们可以在)log(2nn时间内,求出次数为n的多项式在1n次单位复数根处的值.3、利用FFT计算逆DFT将结果的点值表达形式化为系数表达形式(1)在单位复数根处插值根据(2-1-3-a),我们可以把DFT写成矩阵乘积aVyn1,其中nV是一个由1n适当幂次填充成的Vandermonde矩阵.nyyyy210=2121121412112111111111nnnnnnnnnnnnnnnaaaa210(3-3-1-1)易知01nV.即1nV可逆.下面证明11nV),(ji处元素11nvijnij.(3-3-1-2)考虑111nnVV),(ji处元素ijp:nkijknnkkjniknijnnp0101111.(3-3-1-3)当ij时,1ijp;当ij时,根据(2-2-1-d)可知,0ijp.满足nnnIVV111.(3-3-1-4)给定逆矩阵11nV,可以求出Tnaaaa),,,(10.其中nkkinkiyna0111.(3-3-1-5)(2)利用FFT计算逆DFT比较(3-3-1-5)和(3-2-2-2)可知,对FFT算法进行如下修改,即可计算出逆DFT:把a与y互换,用1替换,并且将计算结果每个元素除以n.因此,我们也可以在)log(2nn时间内计算出逆DFT.四、算法具体流程1、加倍多项式次数通过加入n个系数为0的高阶项,把多项式)(xA和)(xB变为次数为n2的多项式,并构造其系数表达.2、求值通过应用)1(2n阶的FFT计算出)(xA和)(xB长度为)1(2n的点值表达.这些点值表达中包含了两个多项式在)1(2n次单位根处的取值.3、逐点相乘把)(xA的值与)(xB的值逐点相乘,可以计算出)()()(xBxAxC的点值表达,这个表示中包含了)(xC在每个)1(2n次单位根处的值.4、插值通过对)1(2n个点值应用FFT,计算其逆DFT,就可以构造出多项式)(xC的系数表达.由于1、3的时间复杂度为)(n,2、4的时间复杂度为)log(2nn,因此整个算法的时间复杂度为)log(2nn.五、算法的实际应用:计算大整数乘法在密码学等领域中,经常需要进行大整数乘法运算.如果对整数qp,逐位相乘然后相加,其时间复杂度为)log(log1010qp.当qp,规模巨大时,这种算法将会十分低效.因此,我们采取快速傅里叶变换进行优化.设01111101010aaaapnnnn,其中NiniNaaii,0,,90(5-1)01111101010bbbbqmmmm,其中NimiNbbii,0,,90(5-2)令多项式01111)(axaxaxaxAnnnn,(5-3)01111)(bxbxbxbxBmmmm.(5-4))()()(xBxAxC.(5-5)注意到)10(Ap,)10(Bq.因此)10()10()10(BACpq.(5-6)将大整数相乘转化为多项式的乘法,应用快速傅里叶变换,即可得出结果.六、参考文献1、《大学数学学习指南——线性代数》(第二版),山东大学出版社,刘建亚,吴臻,秦静,史敬涛,许闻天,张天德,金辉,胡发胜,宿洁,崔玉泉,蒋晓芸编著2、《大学数学线性代数》(第二版),高等教育出版社,上海交通大学数学系线性代数课程组编著3、IntroductiontoAlgorithms(ThirdEdition),ThomasH.Cormen,CharlesE.Leiserson,RonaldL.Rives