内容提要 多项式的表示方式 多项式与FFT 清华大学 多项式的基本运算 宋斌恒 运算的效率分析 FFT思想与实现 多项式 多项式加法 B(x) 在代数域F上的变量为x的多项式 n是幂次( degree)的上界,我们称k是上述多项 式的幂次,如果幂次最大的非零系数是a C(x) C=A+B 例1 若A(x)=6x3+7x2-10X+9 B(x)=2X3+4X-5, 则C(x)=4x3+7x2-6x+4 a,+b, for j
1 多项式与FFT 清华大学 宋斌恒 内容提要 • 多项式的表示方式 – 系数 – 值 • 多项式的基本运算 – 加减乘除 • 运算的效率分析 • FFT思想与实现 多项式 • 在代数域F上的变量为x的多项式 • n是幂次（degree）的上界，我们称k是上述多项 式的幂次，如果幂次最大的非零系数是ak， ∑ − = = 1 0 ( ) n j j j A x a x 多项式加法 ∑ − = = 1 0 ( ) n j j j B x b x ∑ − = = 1 0 ( ) n j j j C x c x C=A+B c = a + b for j = 0,...,n −1 j j j 示例1 • 若A(x) = 6x3 + 7x2 - 10x + 9 • B(x) = -2x3+ 4x - 5, • 则 C(x) = 4x3 + 7x2 - 6x + 4
示例2:乘法 乘法公式 6x3+7x2-10x+9 C(x)=A(x)BO 30x3-35×2+50X-45 24×4+28x3-40x2+36x 12x6-14x5+20x4-18x3 ab-k,j=0.…2n-2 12x6-14x5+44x4-20x3-75x2+86X-45 多项式表示法 系数表示法宜于求值和加法 算A()在点X0的值A(x Horners rule方法计算时间复杂度(mn =ao+x0(a,+X0(a2+…+X0(an2+X0an)…) 加法就是向量相加 a=(a0,a1,an.可以用一个n维向量来表示 而对于乘法 点值表示法 对于由系数表示的多项式A(x)和B(x),其 点值表示( point-value representation) 乘法复杂度是0(n2 幂次为n-1的多项式A(x)可以由点值对 其积的系数向量c,被称为是a和b的卷积 point-value pairs convolution并记做c=at f(xoyb),(x1,yb,,(xn3yn)来表示 由于多项式乘法和卷积运算都是基本运 其中xk互不相等,yk=A(X) 算,我们重点来计算卷积的高效率算法
2 示例2：乘法 • 6x3 + 7x2 - 10x + 9 • - 2x3 + 4x - 5 • ------------------------------------------------------------ • - 30x3 - 35x2 + 50x - 45 • 24x4 + 28x3 - 40x2 + 36x • - 12x6 - 14x5 + 20x4 - 18x3 • ------------------------------------------------------------ • - 12x6 - 14x5 + 44x4 - 20x3 - 75x2 + 86x - 45 乘法公式 C(x) = A(x)B(x) , 0,...,2 2 1 0 =∑ = − − = c a b − j n n k j k j k 多项式表示法 • a = (a0, a1, . . ., an-1). 可以用一个n维向量来表示。 ∑ − = = 1 0 ( ) n j j j A x a x 系数表示法宜于求值和加法 • 计算 A(x) 在点 x0 的值 A(x0). • Horner‘s rule方法计算时间复杂度 θ (n) : • A(x0) = a0 + x0(a1+x0(a2+ ... +x0(an-2+x0(an-1))...)). • 加法就是向量相加 而对于乘法 • 对于由系数表示的多项式 A(x) 和 B(x) ，其 乘法复杂度是 θ (n2), • 其积的系数向量 c, 被称为是a和b的卷积 convolution 并记做 c = a*b. • 由于多项式乘法和卷积运算都是基本运 算，我们重点来计算卷积的高效率算法 点值表示法 • 点值表示（ point-value representation） 幂次为n-1的多项式 A(x) 可以由点值对 （ point-value pairs ） • {(x0, y0), (x1, y1), . . ., (xn-1, yn-1)} 来表示 • 其中 xk 互不相等，yk = A(xk)
进一步思路 点表示定理 个问题可能有多种等价表示方法。而不 同的表示方法对思考、解决问题有极大的 影 =Ax对于k=0 ·我们可以从各种角度去观察、度量一个物 体,只要我们的观察数据可以完全唯一地 反映这个物体的本身,我们就可以认为这 些度量就是这个物体的一个表示。 证明 Vandermonde矩阵 Vandermonde矩阵 ∏x-x) Vo xi x8-1 eterminant【行列 a1 式】不为零 r2 和系数表示是等价 式 Lagrange公式分析 由此引发问题 复杂度(n2) 科学计算需要注意的问题 稳定性:如果点选择不好,结果会非常 一正确性 精确度 截断误差 计算溢出 差放大
3 进一步思路 • 一个问题可能有 多种等价表示方法。而不 同的表示方法对思考、解决问题有极大的 影响。 • 我们可以从各种角度去观察、度量一个物 体，只要我们的观察数据可以完全唯一地 反映这个物体的本身，我们就可以认为这 些度量就是这个物体的一个表示。 点表示定理 • 对于 {(x0, y0), (x1, y1), . . ., (xn-1, yn-1)} 存在 唯一 幂次的上界为n的多项式A(x) 使得 • yk = A(xk) 对于 k = 0, 1, . . ., n – 1成立. 证明： Vandermonde Vandermonde 矩阵 • 上面提到的矩阵是 Vandermonde 矩阵， 它可逆，其 determinant【行列 式】不为零。 • 以上事实表明，点表 示和系数表示是等价 的。 • Lagrange公式： Lagrange Lagrange 公式分析 • 复杂度 θ(n2) • 稳定性：如果点选择不好，结果会非常 差！！！ 由此引发问题 • 科学计算需要注意的问题 – 正确性 – 稳定性 – 精确度 – 截断误差 – 舍入误差 – 计算溢出 – 误差放大 – ……
点表示法的操作 点表示的乘法 加法:只要是 degree- bound【幂次的上 类似地点表示法也可以方便地进行乘法 界】一样的具有相同点表示的多项式的加 若C(x)=A(x)B(x),则C(x)=A(XB()对 法,可以直接把其相关点的值相加即可。 于所有Xk都成 其复杂度为e(n) ·但我们必须 个问题:即C的 degree bound是A和B的和。 我们需要扩充A和B的点值对。 乘法(续) 系数表示多项式的快速乘法 ·扩充A的点对 戈们可以利用线形复杂度的点对表示多项 {(xXoy6),(X1,y1)…,(x2n1y2n1 式的乘法来实现系数表示多项式的快速乘 扩充B的点对 {(xX0y0),(x1y1)…(x2n1,y2n) 计算给定点集的值【如果选用合适点,其复杂 度engn)】 则C的点对表示 计算点值的乘法【e(n)】 ((o,yoyo), (X,,y1),.(X2n-1,y2,-1y2n-1)1 点值表示插值成系数表示法【如果选用合适算 其复杂度e(nlgn)】 示意图 DFT和FFT离散和快速傅立叶变换 ldpc 对于多项式,我们取x 的值为 ,其中a吨 Evalation Tene Bia g n) Tmn用 的第一个根 Pume-valce the principal nth root of unity: n阶主根 -i i ag
4 点表示法的操作 • 加法：只要是degree－bound【幂次的上 界】一样的具有相同点表示的多项式的加 法，可以直接把其相关点的值相加即可。 • 其复杂度为θ(n) 点表示的乘法 • 类似地点表示法也可以方便地进行乘法， 若 C(x) = A(x)B(x), 则 C(xk) = A(xk)B(xk) 对 于所有 xk都成立 • 但我们必须面临一个问题：即 C的degree －bound是 A 和 B 的和。 • 我们需要扩充A和B的点值对。 乘法（续） • 扩充 A的点对, • {(x0,y0),(x1, y1),...,(x2n-1,y2n-1)}, • 扩充 B的点对 • {(x0,y' 0),(x1,y' 1),...,(x2n-1,y' 2n-1)}, • 则C 的点对表示 • {(x0,y0y' 0),(x1,y1y' 1),...,(x2n-1,y2n-1y' 2n-1)}. 系数表示多项式的快速乘法 • 我们可以利用线形复杂度的点对表示多项 式的乘法来实现系数表示多项式的快速乘 法。 – 计算给定点集的值【如果选用合适点，其复杂 度θ(nlgn）】 – 计算点值的乘法【θ(n）】 – 点值表示插值成系数表示法【如果选用合适算 法，其复杂度θ(nlgn）】 示意图 DFT和FFT离散和快速傅立叶变换 • 对于多项式，我们取x 的值为xi =ωn i ，其中ωn 是方程 • xn＝1 • 的第一个根 ωn ＝e2πι/n • the principal nth root of unity :n阶主根
主根的一些性质 DET 多项式 A4(x) (o+m2)2=(o)2 计算其在点:0,o0……,on,上的值 ∑(o)=0如果k不是n的倍数 a Discrete Fourier Transform(DFT) Fast Fourier Transform(FFT) ·我们称y=(6,y1,…,ynd)是系数向量a 利用EF可以在nm)时间内计算DFTn(a)而 aoa1…,an)的离散傅立叶变换,记 FFT使用了 divide-and-conquer策略,把A(x)的 y= DFT,(a) 系数分为奇偶标号两类,然后定义两个 degree- bound为n2的多项式A0(x)和A(x) A0(x) t ax t A1(x) FFT续 RECURSIVE-FFT(a b n is a power of 2. then return a 故有 A(x)=A0(x2)+xA1(x2 6a←(a0,a2,…,an-2) 因此计算A(x)化简为只需计算A0(x)和 A1(x)在点 RECURSIVE-FFT(all (o0)2,(on1)2 10fork←0ton/2-1 上的值 yfo+oyk 14 return y b y is assumed to be column vector
5 主根的一些性质 ( ) 0如果k不是n的倍数 ( ) ( ) 1 1 0 / 2 2 2 2 / 2 = = = = − = ∑ − = + n j k j n k n k n n n n k n dk dn ω ω ω ω ω ω ω DFT ∑ ∑ − = − − = = = 1 0 1 n 1 n 0 n 1 0 y , ,......, , ( ) n j kj k j n n n j j j a A x a x ω 计算其在点：ω ω ω 上的值 多项式 Discrete Fourier Transform (DFT) Discrete Fourier Transform (DFT) • 我们称 y = (y0, y1, . . . , yn-1) 是 系数向量a = (a0, a1, . . . , an-1) 的离散傅立叶变换，记 为： • y = DFTn(a). Fast Fourier Transform (FFT) Fast Fourier Transform (FFT) • 利用FFT可以在 θ(n lg n)时间内计算 DFTn(a) 而 不是直接法的 θ (n2) • FFT 使用了 divide-and-conquer 策略, 把 A(x) 的 系数分为奇偶标号两类，然后定义两个 degreebound 为 n/2 的多项式 A[0](x) 和 A[1](x): FFT续 • 故有 • A(x) = A[0](x2) + xA[1](x2), • 因此计算A(x)化简为只需计算 A[0](x) 和 A[1](x)在点 • （ωn 0）2,（ωn 1）2,……., （ωn n－1）2， • 上的值
复杂度 求解方程,计算a T(n)=27(m2)+n) →7(m)=6nlgn 7-1 od o o a-1 说明 定理 y=Vna,其中Vn是 Vandermonde矩阵x 对于任意k=0 的取值为an的幂次方 于是 第K)位置的 数是n a=DFT(y)=V1(y) 证明:如右图计 当时为1否则 结论 计算方法 利用FFT思想,我们同样可以得到复杂度为 e(nlgn)的快速计算方法 ∑yn 6
6 复杂度 • T(n)=2T(n/2)+θ(n) • ÎT(n)= θ(n lgn) 求解方程，计算a 说明 • y = Vna, 其中 Vn 是 Vandermonde 矩阵其x 的取值为 ωn 的幂次方 • 于是 • a＝DFTn -1(y)＝Vn -1(y) 定理 • 对于任意 j, k = 0, 1, . . . , n - 1, Vn -1 第 (j, k) 位置 的 数是 ωn -jk/n. • 证明：如右图计 算Vn -1Vn=I可得： • 当j=j’时为1否则 为0。 结论 计算方法 • 利用FFT思想，我们同样可以得到复杂度为 θ(n lg n)的快速计算方法
定理 FFT高效实现 任意给定阶为n的向量a和b,其中n是2的密,有 由于FFT在信号处理等方面有相当大的需 a8b= DFT(DFTn(a)-DFT2n(b) 求,故其高效实现有非常重要的意义。 过精心设计可以让复杂度系数变得更小 其中我们把a和b前面用0补足到2n长。 著名的蝴蝶操作:左边两个输 分治的递归调用图【n=8】 入,y1乘以onk 一 循环实现 FFT-BASE(A) 1n∈ length(a]∥ n is a power of2. 递归实现中的10-13 2fors∈1 to Ig n do 行的优化:相同乘法 do Ioy 5tork∈oton-1 by m do orj∈otom2-1co
7 定理 • 任意给定阶为n的向量a和b，其中n是2的密，有 • 其中我们把a和b前面用0补足到2n长。 FFT高效实现 • 由于FFT在信号处理等方面有相当大的需 求，故其高效实现有非常重要的意义。通 过精心设计可以让复杂度系数变得更小。 著名的蝴蝶操作：左边两个输 入 ，yk [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 Å 2s • 4 ωm Å e2π i/m • 5 for k Å 0 to n - 1 by m do • 6 ω Å 1 • 7 for j Å 0 to m/2 - 1 do • 8 t Å ω 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
ERATIVE-FFT(a Bit-reverse位反转] IT-REVERSE-COPY(a, A) n length(a] /n is a power of 2 ·001→100 011110 100→001 k的位反转记为rev(k) Aku+ 把lgn位长度的数字按照二进制数进行位反 转,递归图的顺序排列就是一个位反转排 BIT-REVERSE-COPY(a, A) Parallel ffT to n-1 do 3A[rev(∈a 总结 多项式的表示可以有多种 多项式的乘法可以利用FFT加快计算速度 FFT是信号处理中非常普遍的运算,其速度 直接影响许多仪器的性能 ·换一种思路,开阔你的眼界 ·FFT的实现细节
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 Å 2s • 5 ωm Å e2 π I /m • 6 ω Å 1 • 7 for j Å 0 to m/2 - 1 do • 8 for k j to n - 1 by m do • 9 t Å ω 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[ 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)] Å ak Parallel FFT 总结 • 多项式的表示可以有多种 • 多项式的乘法可以利用FFT加快计算速度 • FFT是信号处理中非常普遍的运算，其速度 直接影响许多仪器的性能 • 换一种思路，开阔你的眼界 • FFT的实现细节