多项式与快速傅里叶变换
记录重要的,有利于理解的部分
多项式的表示
系数表达
对一个次数界(degree bound)为n的多项式A ( x ) = ∑ j = 0 n − 1 a j x j A(x) = \sum_{j=0}^{n-1}a_{j}x^{j}A(x)=∑j=0n−1ajxj而言,其系数表达是一个由系数组成的向量
a = ( a 0 , a 1 , . . . , a n − 1 ) a = (a_{0},a_{1},...,a_{n-1})a=(a0,a1,...,an−1)
求值运算
- 霍纳法则:时间复杂度 O ( N ) O(N)O(N)
乘法运算
卷积(convolution):c = a ⊗ b c = a \otimes bc=a⊗b
时间复杂度:O ( N ) O(N)O(N)
点值表达
一个表达式可以有很多不同的点值表达,采取n个不同的点x 0 , x 1 , . . . , x n x_{0},x_{1},...,x_{n}x0,x1,...,xn构成的集合
{ ( x 0 , y 0 ) , ( x 1 , y 1 ) , . . . , ( x n − 1 , y n − 1 ) } \lbrace(x_{0},y_{0}),(x_{1},y_{1}),...,(x_{n-1},y_{n-1})\rbrace{(x0,y0),(x1,y1),...,(xn−1,yn−1)}
其中,
y k = A ( x k ) y_{k} = A(x_{k})yk=A(xk)
插值运算:求值运算的逆,从一个多项式的点值表达确定其系数表达形式。
插值多项式具有唯一性
解法一:列出范德蒙德矩阵,利用LU分解算法,求解时间复杂度O ( n 3 ) O(n^{3})O(n3)
解法二:基于拉格朗日公式,求解时间复杂度O ( n 2 ) O(n^{2})O(n2)
A ( x ) = ∑ k = 0 n − 1 y k ∏ j ≠ k ( x − x j ) ∏ j ≠ k ( x k − x j ) A(x) = \sum^{n-1}_{k = 0}y_{k}\frac{\prod_{j\not=k}(x-x_{j})}{\prod_{j\not=k}(x_{k}-x_{j})}A(x)=k=0∑n−1yk∏j=k(xk−xj)∏j=k(x−xj)
系数形式表示的多项式的快速乘法
我们可以采用任何点作为求值点,但通过精心地挑选求值点,可以把两种表示法之间转化所需的时间复杂度变为O ( n l g n ) O(nlgn)O(nlgn)。
- DFT:离散傅里叶变换
- 逆DFT
细节:两个次数界为n的多项式乘积是一个次数界为2n的多项式,因此,在对输入多项式A和B进行求值以前,首先把这两个多项式在高次系数前添加n个0,使其次数界加倍为2n。因为现在这些向量有2n个元素,我们可以采用 “2n次单位复数根” ,记为 ω 2 n \omega_{2n}ω2n。
过程:
- 加倍次数界
- 求值(通过应用2n阶的FFT计算出A(x)和B(x)的长度为2n的点值表达)
- 逐点相乘
- 插值(通过对2n个点值对应用FFT,计算其逆DFT)
DFT与FFT
上一节提到通过采用“2n次单位复数根”点可以将时间复杂度降低到O ( n l g n ) O(nlgn)O(nlgn)。
单位复数根
n次单位复数根是满足ω n = 1 \omega^{n} = 1ωn=1的复数ω \omegaω。n次单位复数根恰好有n个:对于k = 0 , 1 , . . . , n − 1 k = 0,1,...,n-1k=0,1,...,n−1,这些根分别是是ω 2 π i k / n \omega^{2\pi ik/n}ω2πik/n。
n个单位复数根均匀地分布在以复平面的原点为圆心的单位半径圆周上。值ω n = e 2 π i / n \omega_{n} = e^{2\pi i/n}ωn=e2πi/n称为主n次单位根,所有其他n次复数单位根都是ω n \omega_{n}ωn的幂次(0次,1次,…,n-1次,n次;其中 ω n 0 = ω n n \omega_{n}^{0} = \omega_{n}^{n}ωn0=ωnn)。
引理30.3(消去引理):对任何整数n ≥ 0 , k ≥ 0 , d > 0 , n\ge0,k\ge0,d>0,n≥0,k≥0,d>0,
ω d n d k = ω n k \omega^{dk}_{dn} = \omega^{k}_{n}ωdndk=ωnk
感觉作用就是可以成比例变化。
ω n n / 2 = ω 2 = − 1 \omega_{n}^{n/2} = \omega_{2} = -1ωnn/2=ω2=−1
引理30.5(折半定理):如果 n > 0 n>0n>0 为偶数,那么n个n次单位复数根的平方的集合就是n/2个n/2次单位复数根的集合。
( ω n k ) 2 = ω 2 ( n / 2 ) 2 k = ω n / 2 k , k = 0 , 1 , . . . , n 2 − 1 (\omega^{k}_{n})^{2} = \omega^{2k}_{2(n/2)} = \omega^{k}_{n/2},k = 0,1,...,\frac{n}{2}-1(ωnk)2=ω2(n/2)2k=ωn/2k,k=0,1,...,2n−1( ω n k + n / 2 ) 2 = ω n 2 k + n = ω n 2 k ω n n = ω n 2 k = ω n / 2 k , k = 0 , 1 , . . . , n 2 − 1 (\omega^{k+n/2}_{n})^{2} = \omega^{2k+n}_{n} = \omega^{2k}_{n}\omega^{n}_{n} = \omega^{2k}_{n} = \omega^{k}_{n/2} , k = 0,1,...,\frac{n}{2}-1(ωnk+n/2)2=ωn2k+n=ωn2kωnn=ωn2k=ωn/2k,k=0,1,...,2n−1
- 引理30.6(求和引理):对任意整数n ≥ 1 n\ge 1n≥1和不能被n整除的非负整数k,有
∑ j = 0 n − 1 ( ω n k ) j = 0 \sum_{j=0}^{n-1}(\omega_{n}^{k})^{j} = 0j=0∑n−1(ωnk)j=0
DFT
回顾一下,我们希望计算次数界为n的多项式
A ( x ) = ∑ j = 0 n − 1 a j x j A(x) = \sum_{j=0}^{n-1}a_{j}x^{j}A(x)=j=0∑n−1ajxj
在ω n 0 , ω n 1 , . . . , ω n n − 1 \omega_{n}^{0},\omega_{n}^{1},...,\omega_{n}^{n-1}ωn0,ωn1,...,ωnn−1处的值(即在n个n次单位复数根处)。假设A以系数形式给出:a = ( a 0 , a 1 , . . . , a n − 1 ) a = (a_{0},a_{1},...,a_{n-1})a=(a0,a1,...,an−1)。接下来对k=0,1,…,n-1,定义结果y k y_{k}yk
y k = A ( ω n k ) = ∑ j = 0 n − 1 a j ω n k j y_{k} = A(\omega^{k}_{n}) = \sum_{j=0}^{n-1}a_{j}\omega_{n}^{kj}yk=A(ωnk)=j=0∑n−1ajωnkj
向量 y = ( y 0 , y 1 , . . . , y n − 1 ) y = (y_{0},y_{1},...,y_{n-1})y=(y0,y1,...,yn−1)就是系数向量 a = ( a 0 , a 1 , . . . , a n − 1 ) a = (a_{0},a_{1},...,a_{n-1})a=(a0,a1,...,an−1) 的离散傅里叶变换(DFT)。记为y = D F T n ( a ) y = DFT_{n}(a)y=DFTn(a)。
FFT
通过使用一种称为**快速傅里叶变换(FFT)**的方法,利用复数单位根的特殊性质,可以在O ( n l o g n ) O(nlogn)O(nlogn)时间内计算出D F T n ( a ) DFT_{n}(a)DFTn(a)。
通篇假设n恰好是2的整数幂。
FFT利用了分治策略,采用A(x)中偶数下标的系数与奇数下标的系数,分别定义两个新的次数界为n/2的多项式 A [ 0 ] ( x ) A^{[0]}(x)A[0](x) 和 A [ 1 ] ( x ) A^{[1]}(x)A[1](x):
A [ 0 ] ( X ) = a 0 + a 2 x + a 4 x 2 + . . . + a n − 2 x n / 2 − 1 A^{[0](X)} = a_{0} + a_{2}x + a_{4}x^{2} + ... + a_{n-2}x^{n/2-1}A[0](X)=a0+a2x+a4x2+...+an−2xn/2−1
A [ 1 ] ( X ) = a 1 + a 3 x + a 5 x 2 + . . . + a n − 1 x n / 2 − 1 A^{[1](X)} = a_{1} + a_{3}x + a_{5}x^{2} + ... + a_{n-1}x^{n/2-1}A[1](X)=a1+a3x+a5x2+...+an−1xn/2−1
于是有:
A ( x ) = A [ 0 ] ( x 2 ) + x A [ 1 ] ( x 2 ) A(x) = A^{[0]}(x^{2}) + xA^{[1]}(x^{2})A(x)=A[0](x2)+xA[1](x2)
所以,求A(x)在ω n 1 , ω n 2 , . . . , ω n n − 1 \omega_{n}^{1},\omega_{n}^{2},...,\omega_{n}^{n-1}ωn1,ωn2,...,ωnn−1处的值的问题转换为:
求次数界为n/2的多项式 A [ 0 ] ( x ) A^{[0]}(x)A[0](x) 和 A [ 1 ] ( x ) A^{[1]}(x)A[1](x)在点
( ω n 0 ) 2 , ( ω n 1 ) 2 , . . . , ( ω n n − 1 ) 2 (\omega_{n}^{0})^{2},(\omega_{n}^{1})^{2},...,(\omega_{n}^{n-1})^{2}(ωn0)2,(ωn1)2,...,(ωnn−1)2
的取值。根据 A ( x ) = A [ 0 ] ( x 2 ) + x A [ 1 ] ( x 2 ) A(x) = A^{[0]}(x^{2}) + xA^{[1]}(x^{2})A(x)=A[0](x2)+xA[1](x2) 综合结果。
根据折半定理,式( ω n 0 ) 2 , ( ω n 1 ) 2 , . . . , ( ω n n − 1 ) 2 (\omega_{n}^{0})^{2},(\omega_{n}^{1})^{2},...,(\omega_{n}^{n-1})^{2}(ωn0)2,(ωn1)2,...,(ωnn−1)2并不是由n个不同的值组成,而是由n/2个n/2次单位复数根所组成,每个根正好出现2次。
因此,可以递归地对次数界为n/2的多项式 A [ 0 ] ( x ) A^{[0]}(x)A[0](x) 和 A [ 1 ] ( x ) A^{[1]}(x)A[1](x) 在n/2个n/2次单位复数根处进行求值。这些子问题与原始问题形式相同,但规模变为一半。
现在,我们以及成功地把一个n个元素的D F T n DFT_{n}DFTn 计算划分为两个规模为n/2的元素的D F T n / 2 DFT_{n/2}DFTn/2 计算。
这一分解是下面递归FFT算法的基础,此算法计算出一个由n个元素组成向量 a = ( a 0 , a 1 , . . . , a n − 1 ) a = (a_{0},a_{1},...,a_{n-1})a=(a0,a1,...,an−1) 的DFT,其中,n是2的幂。
2 ~ 3行:只有一个元素的DFT就是该元素本身
y 0 = a 0 ω 1 0 = a 0 ⋅ 1 = a 0 y_{0} = a_{0}\omega_{1}^{0} = a_{0} \cdot1 = a_{0}y0=a0ω10=a0⋅1=a0
6 ~ 7行:定义 A [ 0 ] ( x ) A^{[0]}(x)A[0](x) 和 A [ 1 ] ( x ) A^{[1]}(x)A[1](x) 的系数向量
4、5、13行:确保 ω \omegaω 的更新
8 ~ 9行:执行递归计算 D F T n / 2 DFT_{n/2}DFTn/2
对于k = 0,1,2,…,n/2-1,
y k [ 0 ] = A [ 0 ] ( ω n / 2 k ) y k [ 1 ] = A [ 1 ] ( ω n / 2 k ) y_{k}^{[0]} = A^{[0]}(\omega^{k}_{n/2})\\ y_{k}^{[1]} = A^{[1]}(\omega^{k}_{n/2})yk[0]=A[0](ωn/2k)yk[1]=A[1](ωn/2k)
或者,根据消去引理,有ω n / 2 k = ω n 2 k \omega^{k}_{n/2} = \omega^{2k}_{n}ωn/2k=ωn2k ,于是
y k [ 0 ] = A [ 0 ] ( ω n 2 k ) y k [ 1 ] = A [ 1 ] ( ω n 2 k ) y_{k}^{[0]} = A^{[0]}(\omega^{2k}_{n})\\ y_{k}^{[1]} = A^{[1]}(\omega^{2k}_{n})yk[0]=A[0](ωn2k)yk[1]=A[1](ωn2k)
对y 0 , y 1 , . . . , y n / 2 − 1 , y_{0},y_{1},...,y_{n/2-1},y0,y1,...,yn/2−1, 第11行推出:
y k = y k [ 0 ] + ω n k y k [ 1 ] = A [ 0 ] ( ω n 2 k ) + ω n k A [ 1 ] ( ω n 2 k ) = A ( ω n k ) ( x = ω n k ) y_{k} = y_{k}^{[0]} + \omega_{n}^{k}y_{k}^{[1]} = A^{[0]}(\omega_{n}^{2k}) + \omega_{n}^{k}A^{[1]}(\omega_{n}^{2k}) \\= A(\omega_{n}^{k})\\ (x = \omega_{n}^{k})yk=yk[0]+ωnkyk[1]=A[0](ωn2k)+ωnkA[1](ωn2k)=A(ωnk)(x=ωnk)
对 y n / 2 , y n / 2 + 1 , . . . , y n − 1 , y_{n/2},y_{n/2+1},...,y_{n-1},yn/2,yn/2+1,...,yn−1, 第12行推出:
y k + ( n / 2 ) = y k [ 0 ] − ω n k y k [ 1 ] = A [ 0 ] ( ω n 2 k ) + ω n k + ( n / 2 ) A [ 1 ] ( ω n 2 k ) = A [ 0 ] ( ω n 2 k + n ) + ω n k + ( n / 2 ) A [ 1 ] ( ω n 2 k + n ) = A ( ω n k + ( n / 2 ) ) ( x = ω n k + ( n / 2 ) ) y_{k+(n/2)} = y_{k}^{[0]} - \omega_{n}^{k}y_{k}^{[1]} \\= A^{[0]}(\omega_{n}^{2k}) + \omega_{n}^{k+(n/2)}A^{[1]}(\omega_{n}^{2k}) \\= A^{[0]}(\omega_{n}^{2k+n}) + \omega_{n}^{k+(n/2)}A^{[1]}(\omega_{n}^{2k+n}) \\= A(\omega_{n}^{k+(n/2)}) \\(x = \omega_{n}^{k+(n/2)})yk+(n/2)=yk[0]−ωnkyk[1]=A[0](ωn2k)+ωnk+(n/2)A[1](ωn2k)=A[0](ωn2k+n)+ωnk+(n/2)A[1](ωn2k+n)=A(ωnk+(n/2))(x=ωnk+(n/2))
其中,ω n k + ( n / 2 ) = − ω n k \omega_{n}^{k+(n/2)} = -\omega_{n}^{k}ωnk+(n/2)=−ωnk , ω n 2 k = ω n 2 k + n \omega_{n}^{2k} = \omega_{n}^{2k+n}ωn2k=ωn2k+n , ω n k \omega_{n}^{k}ωnk :旋转因子
因此,通过采用快速傅里叶变换,可以在 O ( n l g n ) O(nlgn)O(nlgn) 时间内,求出次数界为n的多项式在n个n次单位复数根处的值。
在单位复数根处插值
将多项式从点值表达转换回系数表达。
可以得到:对FFT算法进行 如下修改 就可以计算出 逆DFT :
- 把a与y互换
- 用 ω n − 1 \omega_{n}^{-1}ωn−1 替换 ω n \omega_{n}ωn
- 并将计算结果的每个元素除以n
因此,也可以在 O ( n l g n ) O(nlgn)O(nlgn) 时间内计算出 D F T n − 1 DFT^{-1}_{n}DFTn−1 。
- 定理30.8(卷积定理):对任意两个长度为n的向量a和b,其中n是2的幂,其中向量a和b用0填充,使其长度达到2n,并用” ⋅ \cdot⋅ “表示2个2n个元素组成向量的点乘。
a ⊗ b = D F T 2 n − 1 ( D F T 2 n ( a ) ⋅ D F T 2 n ( b ) ) a \otimes b = DFT_{2n}^{-1}(DFT_{2n}(a)\cdot DFT_{2n}(b))a⊗b=DFT2n−1(DFT2n(a)⋅DFT2n(b))
高效FFT实现
探究两种高效的FFT实现方法。
FFT的一种迭代实现
第10 ~ 13行的 for 循环中包含了 ω n k y k [ 1 ] \omega_{n}^{k}y_{k}^{[1]}ωnkyk[1] 的2次计算。称该值为 公用制表达式 。可以改变循环,使其只计算一次,将其值放在临时变量 t tt 中。
这一系列操作称为一个 蝴蝶操作 。
现在来说明如何使FFT算法采用迭代结构而不是递归结构。
树中每一个结点对应每次过程递归调用。由相应的输入向量标记。每次Recursive-FFT 调用产生两个递归调用,除非该调用已收到了1个元素的向量。第一次调用作为左孩子,第二次调用作为右孩子。
追踪过程,自底向上。
- 首先,成对取出元素,利用一次蝴蝶操作计算出每对的DFT,然后用其DFT取代这对元素。这样向量中就包含了n/2个二元素的DFT。(第四层到第三层)
- 下一步,按对取出这n/2个DFT,通过两次蝴蝶操作计算出具有四个元素向量的DFT,并用一个具有四个元素的DFT取代对应的两个二元素的DFT。于是,向量中包含n/4个四元素的DFT。
- 继续进行这一过程,直至向量包含两个具有n/2个元素的DFT,这时,综合应用n/2次蝴蝶操作,可以合成最终的具有n个元素的DFT。(第二层到第一层)
数组A[0,…,n-1]:初始时该数组包含输入向量a中的元素,其顺序为它们在上图中树叶出现的顺序(确定顺序:位逆序置换 )。
引入变量s:计算树的层次。
更精确的伪代码:
共三层循环:
- 第一层循环:s ss 确定所在的层数,m = 2 s m = 2^{s}m=2s ,要求 m个数的DFT ,ω m = e 2 π i / m \omega_{m} = e^{2\pi i/m}ωm=e2πi/m 为 m mm 次单位复数根(通过迭代乘以 ω \omegaω 得到 ω n k \omega_{n}^{k}ωnk (k次的))。
- 第二层循环:由于每组有m个数(k以m递增),故有n/m组。初始 ω n 0 = 1 \omega_{n}^{0} = 1ωn0=1。
- 第三层循环:蝴蝶操作 ,每组m个数,每组m/2次蝴蝶操作 (结合第二层循环:一共n/m组,共n/2次蝴蝶操作)
位逆序置换
让 r e v ( k ) rev(k)rev(k) 为 k kk 的二进制表示各位逆序所形成的 l g n lgnlgn 位的整数,把向量中的元素 a k a_{k}ak 放在数组的 A [ r e v ( k ) ] A[rev(k)]A[rev(k)] 位置上。
例如,在上上上图中树叶出现的次序为:0,4,2,6,1,5,3,7;这个序列二进制表示为:000,100,010,110,001,101,011,111;当把二进制表示各位逆序后,得到序列:000,001,010,011,100,101,110,111(递增)。
获得一般情况下的位逆序置换:在树的最底层,最低位为0的下标在左子树,以及最低位为1的下标在右子树中。在每一层去掉最低位后,沿着树继续这一过程,直到给出序列。
总体时间复杂度
并行FFT电路
SEE BOOK.
总结
T_T,感觉还是糊里糊涂的。
贴上FFT代码:
...