利用快速傅里叶变换(FFT)计算多项式乘法作者:宋振华摘要本文将讨论快速傅里叶变换(FFT),利用FFT 设计一种算法,使多项式相乘的时间复杂度降低为()n n 2log Θ,以便在计算机上高效计算多项式乘法.关键词:快速傅里叶变换、多项式乘法目录一、引言 二、算法概述 三、引理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 设计一种算法,使多项式相乘的时间复杂度降低为()n n 2log Θ,以便在计算机上高效执行.二、算法概述通过FFT 进行离散傅里叶变换,将两个多项式由系数表达形式转化为点值表达形式.将这2个点值表达形式的多项式相乘,得到结果的点值表达形式.最后利用FFT 做DFT 的逆,得到结果的系数表达形式.三、引理1、多项式的表示(1)系数表达对于次数为n 的多项式()∑==ni i i x a x A 0,其系数表达是一个由系数组成的向量()Tn a a a a ,...,,10=.考虑用系数形式表示、次数为n 的多项式)(x A 、)(x B 的乘法运算. 记∑===ni ii x c x B x A x C 20)()()(. 有∑=-=ij j i j i b a c 0.可以看出,采取逐个计算i c ),20(N i n i ∈≤≤的方式进行求解,其时间复杂度为)(2n Θ. (2)点值表达 a.定义:一个次数为n 的多项式)(x A 的点值表达就是由一个有n 个点值对组成的集合)}(,,0|),{(i i i i x A y N i n i y x =∈≤≤,使得对于n k ,...,2,1,0=,所有的k x 各不相同.b.点值表达下多项式乘法对于n 次多项式)(x A ,)(x B ,设其点值表达式分别为:)(x A :)},(,),,(),,(),,{(221100n n y x y x y x y x ,)(x B :)},(,),,(),,(),,{(,,2211,00n n y x y x y x y x ，设)()()(x B x A x C =,由于)(x C 的次数为n 2,因此必须对)(x A 和)(x B 的点值表达式进行扩展,使每个多项式都包含12+n 个点值对.给定A ,B 的扩展点值表达:)(x A :)},(,),,(),,(),,{(22221100n n y x y x y x y x ,)(x B :)},(,),,(),,(),,{(,22,2211,00n n y x y x y x y x ，.则)(x C 的点值表达形式为:)},(,),,(),,(),,{(,222,222111,000n n n y y x y y x y y x y y x ，.因此,计算)(x C 点值表达式的时间复杂度为)(n Θ.` (3)插值求值运算的逆(从一个多项式的点值表达形式确定其系数表达形式)称为插值.下面将证明,n 个点求值运算与插值运算是定义完备的互逆运算.a.多项式的点值表达可以唯一确定多项式的系数表达形式.多项式的点值表达等价于矩阵方程:⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛n n n n n nnx x x x x x x x x x x x 22222121102001111⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛n a a a a 210=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛n y y y y 210. (3-1-3-a)记系数矩阵为V ,由Vandermonde 行列式可知,()∏≤<≤-=nj i jix x V 0.而j i n j i N j i ≠≤≤∈∀,,0,,,有j i x x ≠,因此0≠V ,即V 可逆.因此对于给定点值表达,我们能够唯一确定系数表达式,且y V a 1-=.b.对于次数为n 的多项式函数()x A ,插值算法基于如下Lagrange 公式:()()∑∏∏=≠≠--=nk kj jkk j jk x x x x y x A 0)(.(3-1-3-b)容易验证,Lagrange 公式的计算复杂度为)(2n Θ.因此,n 个点求值运算与插值运算是定义完备的互逆运算,它们将多项式的系数表达与点值表达进行相互转化.2、利用离散傅里叶变换(DFT)与FFT 导出结果的点值表达形式(1)单位复数根n 次单位复数根是满足1=n ω的复数ω.n 次单位复数根恰好有n 个:对于1,...,1,0-=n k ,这些根是nk i eπ2.由Euler 公式()()θθθsin cos i e i +=可知,这n 个单位复数根均匀地分布在以复平面原点为圆心的单位圆圆周上.ni n eπω2=称为主n 次单位根,所有其他n 次单位根都是n ω的幂次.n 个n 次单位复数根0nω,1n ω,2n ω,...,1-n n ω在乘法意义下构成一个群,该群与群>+<n n Z ,同构.由此,可以得到如下推论:a.0,0,0>≥≥∀d k n ,有knnk idkdneedndk iωωππ===22. (3-2-1-a) b.*,2N k k n ∈=∀,122-==ωωnn.(3-2-1-b)c.若*,2N k k n ∈=,},0|){(2N i n i i n ∈<≤ω=},20|){(22/N i n i i n ∈<≤ω. (3-2-1-c)对推论c 的证明:由a 可知,N k ∈∀,()2/2/2k n k n ωω=.所以对于20,n k N k <≤∈,有222222/)()(k n k n n n k n n k n n k n ωωωωωω====++.得证.d.N k N n ∈∈∀,*且k 不能被n 整除,有()∑-==100n i nk nω. (3-2-1-d)对推论d 的证明:()∑-=1n i n k nω=11)(--k n n k n ωω=11)(--kn k n n ωω=11)1(--k n k ω=0. (2)离散傅里叶变换(DFT)对于n 次多项式∑==ni i i x a x A 0)(,其系数形式为T n a a a a ),,,(10 =.(3-2-2-1) 设∑=+==ni kin i knk a A y 01)(ωω,N k n k ∈≤≤,0,(3-2-2-2)则向量T n y y y y ),,,(10 =(3-2-2-3)就是系数向量T n a a a a ),,,(10 =的离散傅里叶变换. (3)通过快速傅里叶变换(FFT)计算向量y .对于n 次多项式∑==ni i i x a x A 0)(,不妨假设N p n p ∈=+,21.对于1+n 不为2的整数次幂的情况类似,此处不再讨论.FFT 采取了分治策略,采用)(x A 中偶数下标与奇数下标的系数,分别定义两个新的2n次多项式: []21-124200)(n n x a x a x a a x A -++++= , (3-2-3-1)[]21-25311)(n n xa x a x a a x A ++++= . (3-2-3-2)注意到[])(0x A 包含A 所有偶数下标的系数,[])(1x A 包含A 中所有奇数下标的系数,于是有:[])()()(2]1[20x xA x A x A +=.(3-2-3-3)所以,求)(x A 在01+n ω,11+n ω,...,nn 1+ω处的值的问题转化为:a.求次数为2n 的多项式)(]0[x A ,)(]1[x A 在点201)(+n ω,211)(+n ω,..., 21)(n n +ω处的取值.b.根据(2-2-3-3)综合结果.根据(2-2-1-c),式(2-2-3-3)不是由1+n 个不同值组成,而是仅由21+n 个21+n 次单位复数根所组成,每个根正好出现2次.因此,我们递归地对次数为21+n 的多项式)(]0[x A ,)(]1[x A 在21+n 个21+n 次单位复数根处进行求值.这些子问题与原始问题形式相同,但规模变为一半.下面确定DFT 过程的时间复杂度.注意到除了递归调用外,每次调用需要枚举T n a a a a ),,,(10 =中所有元素,将a 划分为[]},,{200n a a a a =、[]},,,,{15311-=n a a a a a ,分别与多项式)(]0[x A ,)(]1[x A 相对应.其时间复杂度为)(n Θ.因此,对运行时间有下列递推式:)()2(2)(n nT n T Θ+=.(3-2-3-4)求解该递推式,有)log ()(2n n n T Θ=.(3-2-3-5)采取快速傅里叶变换,我们可以在)log (2n n Θ时间内,求出次数为n的多项式在1+n 次单位复数根处的值.3、利用FFT 计算逆DFT 将结果的点值表达形式化为系数表达形式(1)在单位复数根处插值根据(2-1-3-a),我们可以把DFT 写成矩阵乘积a V y n 1+=,其中n V 是一个由1+n ω适当幂次填充成的Vandermonde 矩阵.⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛n y y y y 210=⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛+++++++++2121121412112111111111n n nn n n n n n n nn n n ωωωωωωωωω⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛n a a a a 210 (3-3-1-1)易知01≠+n V .即1+n V 可逆. 下面证明11-+n V),(j i 处元素11+=-+n v ij n ij ω. (3-3-1-2)考虑111+-+n n V V ),(j i 处元素ij p :()∑∑=-+=+-++=⋅+=nk i j k n n k kjn ikn ij n n p 011111ωωω. (3-3-1-3)当i j =时,1=ij p ;当i j ≠时,根据(2-2-1-d)可知,0=ij p .满足n n n I V V =+-+111.(3-3-1-4)给定逆矩阵11-+n V ,可以求出T n a a a a ),,,(10 =.其中∑=-++=n k ki n k i y n a 0111ω.(3-3-1-5)(2)利用FFT 计算逆DFT比较(3-3-1-5)和(3-2-2-2)可知,对FFT 算法进行如下修改,即可计算出逆DFT:把a 与y 互换,用1-ω替换ω,并且将计算结果每个元素除以n .因此,我们也可以在)log (2n n Θ时间内计算出逆DFT.四、算法具体流程1、加倍多项式次数通过加入n 个系数为0的高阶项,把多项式)(x A 和)(x B 变为次数为n 2的多项式,并构造其系数表达. 2、求值通过应用)1(2+n 阶的FFT 计算出)(x A 和)(x B 长度为)1(2+n 的点值表达.这些点值表达中包含了两个多项式在)1(2+n 次单位根处的取值. 3、逐点相乘把)(x A 的值与)(x B 的值逐点相乘,可以计算出)()()(x B x A x C =的点值表达,这个表示中包含了)(x C 在每个)1(2+n 次单位根处的值. 4、插值通过对)1(2+n 个点值应用FFT,计算其逆DFT,就可以构造出多项式)(x C 的系数表达.由于1、3的时间复杂度为)(n Θ,2、4的时间复杂度为)log (2n n Θ, 因此整个算法的时间复杂度为)log (2n n Θ. 五、算法的实际应用:计算大整数乘法在密码学等领域中,经常需要进行大整数乘法运算.如果对整数q p ,逐位相乘然后相加,其时间复杂度为)log (log 1010q p ⋅Θ.当q p ,规模巨大时,这种算法将会十分低效.因此,我们采取快速傅里叶变换进行优化.设01111101010a a a a p n n n n +⋅++⋅+⋅=-- ,其中N i n i N a a i i ∈≤≤∈≤≤,0,,90(5-1)01111101010b b b b q m m m m +⋅++⋅+⋅=-- ,其中N i m i N b b i i ∈≤≤∈≤≤,0,,90(5-2)令多项式01111)(a x a x a x a x A n n n n ++++=-- ,(5-3)01111)(b x b x b x b x B m m m m ++++=-- .(5-4))()()(x B x A x C =. (5-5)注意到)10(A p =,)10(B q =. 因此)10()10()10(B A C pq ==.(5-6)将大整数相乘转化为多项式的乘法,应用快速傅里叶变换,即可得出结果.六、参考文献1、《大学数学学习指南——线性代数》(第二版),山东大学出版社,刘建亚,吴臻,秦静,史敬涛,许闻天,张天德,金辉,胡发胜,宿洁,崔玉泉,蒋晓芸编著2、《大学数学线性代数》(第二版),高等教育出版社,上海交通大学数学系线性代数课程组编著3、Introduction to Algorithms(Third Edition),ThomasH.Cormen,Charles E.Leiserson,Ronald L.Rivest,Clifford Stein4、《离散数学》,山东大学软件学院,徐秋亮,栾峻峰,卢雷,王慧,赵合计编著5、《数学分析 下册》(第二版),高等教育出版社,陈纪修,於崇华,金路编著6、《高等代数》,科学出版社,丘维声编著7、《复变函数论》(第三版),高等教育出版社,钟玉泉编著8、《快速傅里叶变换乘法的性能研究》,陕西师范大学,毛庆,李顺东。
![基于FFT的大整数乘法SSA算法多核并行化实现方法[发明专利]](https://img.taocdn.com/s1/m/e4e4438acaaedd3382c4d34e.png)
