1 多项式与多项式与FFT 清华大学清华大学 宋斌恒宋斌恒 内容提要内容提要 ?多项式的表示方式 –系数 –值 ?多项式的基本运算 –加减乘除 ?运算的效率分析 ? FFT思想与实现 多项式多项式 ?在代数域F上的变量为x的多项式 ? n是幂次（degree）的上界，我们称k是上述多项 式的幂次，如果幂次最大的非零系数是a k ， ∑ ? = = 1 0 )( n j j j xaxA 多项式加法多项式加法 ∑ ? = = 1 0 )( n j j j xbxB ∑ ? = = 1 0 )( n j j j xcxC C=A+B 1,...,0for ?=+= njbac jjj 示例示例1 ?若A(x) = 6x 3 + 7x 2 -10x + 9 ? B(x) = -2x 3 + 4x -5, ?则C(x) = 4x 3 + 7x 2 -6x + 4. 2 示例示例2：乘法：乘法 ? 6x 3 + 7x 2 -10x + 9 ? -2x 3 + 4x -5 ? ------------------------------------------------------------ ? -30x 3 -35x 2 + 50x - 45 ? 24x 4 + 28x 3 -40x 2 + 36x ? -12x 6 -14x 5 + 20x 4 -18x 3 ? ------------------------------------------------------------ ? -12x 6 -14x 5 + 44x 4 -20x 3 -75x 2 + 86x - 45 乘法公式乘法公式 )()()( xBxAxC = 22,...,0 , 1 0 ?== ∑ ? = ? njbac n k kjkj 多项式表示法多项式表示法 ? a = (a 0 , a 1 , . . ., a n-1 ).可以用一个n维向量来表示。 ∑ ? = = 1 0 )( n j j j xaxA 系数表示法宜于求值和加法系数表示法宜于求值和加法 ?计算A(x) 在点x 0 的值A(x 0 ). ? Horner‘s rule方法计算时间复杂度θ (n) : ? A(x 0 ) = a 0 + x 0 (a 1 +x 0 (a 2 + ... +x 0 (a n-2 +x 0 (a n-1 ))...)). ?加法就是向量相加 而对于乘法而对于乘法 ?对于由系数表示的多项式A(x) 和B(x) ，其 乘法复杂度是θ (n 2 ), ?其积的系数向量c, 被称为是a和b的卷积 convolution并记做c = a*b. ?由于多项式乘法和卷积运算都是基本运 算，我们重点来计算卷积的高效率算法 点值表示法点值表示法 ?点值表示（point-value representation） 幂次为n-1的多项式A(x) 可以由点值对 （point-value pairs） ? {(x 0 , y 0 ), (x 1 , y 1 ), . . ., (x n-1 , y n-1 )} 来表示 ?其中x k 互不相等，y k = A(x k ) 3 进一步思路进一步思路 ?一个问题可能有多种等价表示方法。而不 同的表示方法对思考、解决问题有极大的 影响。 ?我们可以从各种角度去观察、度量一个物 体，只要我们的观察数据可以完全唯一地 反映这个物体的本身，我们就可以认为这 些度量就是这个物体的一个表示。 点表示定理点表示定理 ?对于{(x 0 , y 0 ), (x 1 , y 1 ), . . ., (x n-1 , y n-1 )} 存在 唯一幂次的上界为n的多项式A(x) 使得 ? y k = A(x k ) 对于k = 0, 1, . . ., n – 1成立. 证明：证明：Vandermonde矩阵 矩阵 ?上面提到的矩阵是 Vandermonde矩阵， 它可逆，其 determinant【行列 式】不为零。 ?以上事实表明，点表 示和系数表示是等价 的。 ? Lagrange公式： Lagrange 公式分析公式分析 ?复杂度θ(n 2 ) ?稳定性：如果点选择不好，结果会非常 差！！！ 由此引发问题由此引发问题 ?科学计算需要注意的问题 –正确性 –稳定性 –精确度 –截断误差 –舍入误差 –计算溢出 –误差放大 –…… 4 点表示法的操作点表示法的操作 ?加法：只要是degree－bound【幂次的上 界】一样的具有相同点表示的多项式的加 法，可以直接把其相关点的值相加即可。 ?其复杂度为θ(n) 点表示的乘法点表示的乘法 ?类似地点表示法也可以方便地进行乘法， 若C(x) = A(x)B(x), 则C(x k ) = A(x k )B(x k ) 对 于所有x k 都成立 ?但我们必须面临一个问题：即C的degree －bound是A 和B的和。 ?我们需要扩充A和B的点值对。 乘法（续）乘法（续） ?扩充A的点对, ? {(x 0 ,y 0 ),(x 1 , y 1 ),...,(x 2n-1 ,y 2n-1 )}, ?扩充B的点对 ? {(x 0 ,y' 0 ),(x 1 ,y' 1 ),...,(x 2n-1 ,y' 2n-1 )}, ?则C 的点对表示 ? {(x 0 ,y 0 y' 0 ),(x 1 ,y 1 y' 1 ),...,(x 2n-1 ,y 2n-1 y' 2n-1 )}. 系数表示多项式的快速乘法系数表示多项式的快速乘法 ?我们可以利用线形复杂度的点对表示多项 式的乘法来实现系数表示多项式的快速乘 法。 –计算给定点集的值【如果选用合适点，其复杂 度θ(nlgn）】 –计算点值的乘法【θ(n）】 –点值表示插值成系数表示法【如果选用合适算 法，其复杂度θ(nlgn）】 示意图示意图 DFT和和FFT离散和快速傅立叶变换离散和快速傅立叶变换 ?对于多项式，我们取x 的值为x i =ω n i ，其中ω n 是方程 ? x n ＝1 ?的第一个根 null ω n ＝e 2πι/n ? the principal n th root of unity :n阶主根 5 主根的一些性质主根的一些性质 的倍数不是如果nk 0)( )()( 1 1 0 222/ 2 2/ = = ?== = ∑ ? = + n j jk n k n nk n n n k n dk dn ω ωω ωω ωω DFT ∑ ∑ ? = ? ? = = = 1 0 1 n 1 n 0 n 1 0 y ,,......,, )( n j kj njk n n j j j a xaxA ω ωωω上的值计算其在点： 多项式 Discrete Fourier Transform (DFT) ?我们称y = (y 0 , y 1 , . . . , y n-1 ) 是系数向量a = (a 0 , a 1 , . . . , a n-1 ) 的离散傅立叶变换，记 为： ? y = DFT n (a). Fast Fourier Transform (FFT) ?利用FFT可以在θ(n lg n)时间内计算DFT n (a) 而 不是直接法的θ (n 2 ) ? FFT 使用了divide-and-conquer 策略, 把A(x) 的 系数分为奇偶标号两类，然后定义两个degree- bound 为n/2 的多项式A [0] (x) 和A [1] (x): FFT续续 ?故有 ? A(x) = A [0] (x 2 ) + xA [1] (x 2 ), ?因此计算A(x)化简为只需计算A [0] (x) 和 A [1] (x)在点 ?（ω n 0 ） 2 ,（ω n 1 ） 2 ,……., （ω n n－1 ） 2 ， ?上的值 6 复杂度复杂度 ? T(n)=2T(n/2)+θ(n) ? ?T(n)= θ(n lgn) 求解方程，计算求解方程，计算a 说明说明 ? y = V n a, 其中V n 是Vandermonde矩阵其x 的取值为ω n 的幂次方 ?于是 ? a＝DFT n -1 (y)＝V n -1 (y) 定理定理 ?对于任意j, k = 0, 1, . . . , n -1, V n -1 第(j, k) 位置的 数是ω n -jk /n. ?证明：如右图计 算V n -1 V n =I可得： ?当j=j’时为1否则 为0。 结论结论计算方法计算方法 ?利用FFT思想，我们同样可以得到复杂度为 θ(n lg n)的快速计算方法 7 定理定理 ?任意给定阶为n的向量a和b，其中n是2的密，有 ?其中我们把a和b前面用0补足到2n长。 FFT高效实现高效实现 ?由于FFT在信号处理等方面有相当大的需 求，故其高效实现有非常重要的意义。通 过精心设计可以让复杂度系数变得更小。 著名的蝴蝶操作：左边两个输著名的蝴蝶操作：左边两个输 入入，，y k [1] 乘以乘以ω n k 分治的递归调用图分治的递归调用图【【n=8】】 循环实现循环实现 ?递归实现中的10－13 行的优化：相同乘法 只做一次 FFT-BASE(A) ? 1 n ? length[a] // n is a power of 2. ? 2 for s ?1 to lg n do ? 3 m ?2 s ? 4 ω m ?e 2π i/m ? 5 for k ?0 to n - 1 by m do ? 6 ω ?1 ? 7 for j ? 0 to m/2 - 1 do ? 8t ?ω A[k + j + m/2] ? 9 u ? A[ k + j] ? 10 A[k + j] ? u + t ? 11 A[k + j + m /2] ?u - t ? 12 ω ? ωω m 8 ITERATIVE-FFT(a) ? 1 BIT-REVERSE-COPY (a, A) ? 2 n ? length[a] //n is a power of 2. ? 3 for s ?1 to lg n do ? 4 m ?2 s ? 5 ω m ?e 2 π I /m ? 6 ω ? 1 ? 7 for j ?0 to m/2 - 1 do ? 8 for k j to n - 1 by m do ? 9t ?ω A[k + m/2] ? 10 u ?A[k] ? 11 A[k] ?u + t ? 12 A[k + m/2] ?u - t ? 13 ω ?ωω m ? 14 return A Bit-reverse[位反转位反转] ? 001?100 ? 010?010 ? 011?110 ? 100?001 ? ……，k的位反转记为rev（k） ?把lgn位长度的数字按照二进制数进行位反 转，递归图的顺序排列就是一个位反转排 列 BIT-REVERSE-COPY(a, A) ? 1 n ? length[a] ? 2 for k ? 0 to n -1 do ? 3 A[rev(k)] ? a k Parallel FFT 总结总结 ?多项式的表示可以有多种 ?多项式的乘法可以利用FFT加快计算速度 ? FFT是信号处理中非常普遍的运算，其速度 直接影响许多仪器的性能 ?换一种思路，开阔你的眼界 ? FFT的实现细节