首页 > Python算法 > 分治算法 阅读数:40

快速傅立叶变换—分治算法经典例题

快速傅氏变换简称 FFT,用于获得两个多项式的乘积,比如 f(x)=1+5x+3x2+2x3 与 g (x)=10+3x+x5 的乘积为 f(x)g(x)=10+53x+45x2+29x3+6x4+x5+5x6+3x7+2x8。用暴力法计算两个多项式的乘积的时间复杂度为 O(n2),但用 FFT 算法的时间复杂度仅为 O(nlgn),n 为乘积的次数。

我们正式的问题为:给定两个多项式 f(x) 与 g(x) 的系数,输出 f(x)g(x) 的系数。比如:
输入:

f(x) = [1,5,3,2],g(x) = [10,3,0,0,0,1]

输出:

[10,53,45,29,6,1,5,3,2]

请注意在输入和输出中,如果某个 xn 项为 0,则它的系数为 0。另外,系数的顺序从常数项开始,按照 x 的次数升序排序。

在理解 FFT 之前,我们需要先理解多项式的点值表示方式。一般情况下,我们都用系数表示多项式,比如 [10,3,0,0,0,1] 代表 10+3x+x5。按照这个式子,我们能够画出一条独一无二的曲线。但是,除了系数表示方式,多项式还有一种表达方式叫作点值表示方式。

如果多项式的次数为 n,那么多项式的点值表示为 n+1 个相互独立的点值对。例如,因为 10+3x+x5 的次数为 5,所以它的点值表示就是 6 个相互独立的点值对。点值对就是曲线上的 (x,y) 坐标。比如,(1,14) 是 10+3x+x5 的一个点,所以 (1,14) 是 10+3x+x5 的一个点值对。

下面我们将系数表示转化成点值表示:
  • 多项式:y=10+3x+x5
  • 取 6 个点:x=-2,-1,0,1,2,3;
  • 点值表示:(-2,-28),(-1,6),(0,10),(1,14),(2, 48),(3,262)。

以上 6 个点值对是 y=10+3x+x5 的点值表示方式,因为没有第二个多项式同时穿过这 6 个点。

当两个多项式相乘时,我们只需相乘它们的点值对就可以获得乘积的点值表示。例如:f(x)=1+5x+3x2+2x3,g(x)=10+3x+x5,求 f(x)g(x) 的点值表示。

1) 取 9 个点:x=-4,-3,-2,-1,0,1,2,3,4。

2) 计算点值对如下:
  • f(x):(-4,-99),(-3,-41),(-2,-13),(-1,-3),(0,1),(1,11),(2,39),(3,97),(4,197);
  • g(x):(-4,-1026),(-3,-245),(-2,-28),(-1,6),(0,10),(1,14),(2,48),(3,262),(4,1046)。

3) 点值乘积:(-4,101574),(-3,10045),(-2,364),(-1,-18),(0,10),(1,154),(2,1872),(3,25414),(4,202122)。

以上 9 个点值对就是乘积 10+5 5x2+23x+4 9x3+6x4+x5+5x6+3x7+2x8 的点值表示。

需要注意的是,我们在开始之前必须知道需要选择多少个点。也就是说,我们需要提前知道乘积的次数。多项式乘积的次数是两个多项式次数的和。比如如果 f(x) 和 g(x) 的次数分别为 3 和 5,那么它们的乘积的次数就是 3+5=8。

到目前为止,我们得到的结论是:计算两个多项式的点值表示,再将点值相乘,就能获得乘积的点值表示。我们只差一步就能得到答案:将乘积的点值表示转化为系数表示。

FFT 在解决多项式乘法问题中的作用有两个,一是将系数表示快速地转化为点值表示,二是将乘积的点值表示快速地转化为系数表示。FFT 算法是一个独立的算法,通过传入不同的参数,可以解决不同的问题。因此,解决第一部分和解决第二部分时我们需要向 FFT 方法传入不同的参数。

FFT 算法之所以快速,是因为它巧妙地选择了单位根为点。我们先来通过一个例子理解 FFT 怎样解决第一部分:将系数表示快速地转化为点值表示。

h(x)=a0+a1x+a2x2+a3x3+a4x4+a5x5+a6x6+a7x7,h (x) 的系数表示为 [a0,a1,a2,a3,a4,a5,a6,a7]。我们的目的是得到 h(x) 的点值表示。最简单的方法是选择 8 个实数 x 值,代入 x 值得到 y 值。下面的方法是更加快速的 FFT 方法。

声明以下两个数组:

[a0,a2,a4,a6] 包括偶数项的系数,也就是说对应的 xn 项的 n 为偶数。
[a1,a3,a5,a7] 包括奇数项的系数。


令 h1(x) 和 h2(x) 为次数为 3 的多项式:

h1(x)=a0+a2x+a4x2+a6x3
h2(x)=a1+a3x+a5x2+a7x3


注意到:

h(x)=h1(x2)+xh2(x2)
h(-x)=h1(x2)-xh2(x2)


这意味着,只要我们得到 h1(x2) 和 xh2(x2),就能同时得到 h(x) 和 h(-x)。换句话说,我们为了得到 h(x) 而计算的信息,可以被二次利用,用来计算 h(-x)。这样的话,我们付出一份时间能够得到两个值。

因此,我们理所当然地应该选择正负对称的点,比如,如果乘积的次数为 6,我们应选 [-3,-2,-1,0,1,2,3] 为点。这样,计算 f(1)、f(2)、f(3) 的同时即可得到 f( 1)-、f( 2)- 、f( 3)- 。

读者可能会有疑问,我们并不是花了一份时间,而是两份时间,因为我们计算了 h1(x2) 和 h2(x2),这不是两次运算吗?FFT 算法巧妙的地方就在于,我们可以递归地计算 h1(x2) 和 h2(x2),这导致它们的计算时间也被优化。而唯一能够满足重复递归条件的,只有单位根。平方根能够满足在每一次递归时,传入的值都是正负对称的。

例如,令 8 个点为单位的 8 次根:[1,w,i, w3,-1,-w,-i,-w3]。ω 是第一个单位根,为 

在第二次递归的时候,我们会计算 h1(x2),也就是 h1(1)、h1(i)、h1(1)、h1(i) 的值。对我们有利的是,[1,i,-1,-i] 是正负对称的。因此:

h1(x)=a0+a2x+a4x2+a6x3
h11(x)=a0+a4x
h12(x)=a2+a6x


并且:

h1(x)=h11(x2)+xh12(x2)
h1(-x)=h11(x2)-xh12(x2)


在第三次递归的时候,我们会计算 h11(x2),也就是 h11(1) 和 h11(−1) 的值。同样的,[1,-1] 是两个正负对称的数字。在第四次递归的时候,因为只需要计算一个值,所以直接计算并输出答案即可。接下来我们将子答案层层合并,得到 h(x) 的点值表示。

在解决多项式乘法问题时,我们需要传入第二个多项式,将步骤重复一遍,得到第二个多项式的点值表示。

读者可能会觉得这样计算点值表示麻烦,不如随便选择几个点,将 x 值直接代入。但是,如果输入的多项式的项数更大,我们就可以清晰地发现用单元根的优势。更重要的是,我们还有另外一项工作需要单元根的帮助:将点值表示转化为系数表示。

在解释 FFT 算法怎样做第二项工作之前,我们先通过一个例子解释为什么实数不满足重复递归条件。

假设乘积的次数为 7,按照以上原则,我们选择的 8 个点应为 [-4,-3,-2,-1,1,2,3,4]。在第二次递归的时候,也就是计算 h1(x2)的 时候(或者 h1(-x2) 的时候,但是我们只看前者的情况就够了),我们想要得到 h1(1)、h1(4)、h1(9)、h1(16) 的值。但是,我们没有办法继续优化这个过程,因为 1、4、9、16 全都是正数。毕竟传入的是实数的平方,所以我们不可能得到负数。

在做第一项工作的时候,我们用 n 个系数值求 n 个点值。现在做第二项工作,用 n 个点值求 n 个系数值。这两项工作看上去需要用到不同的方法,其实不然。如方程所示,左边第一个矩阵代表单位根,第二个矩阵代表系数,右边的矩阵代表我们想要得到的点值。在求点值表示的时候,我们相当于在求两个矩阵的乘积。

现在,做第二项工作的时候,我们在两边同时乘以单元根矩阵的逆矩阵。得到:

右边的第一个矩阵里面是单位根的倒数,第二个矩阵里是已知的点值,左边的矩阵是我们想要得到的系数。

这意味着,求系数表示的时候,我们可以传入多项式:

f(x)=y0+y1x+y2x2+…+yn-1xn-1

同时,让 1,w-1,w-2,…,w-(n-1) 除以 n 就是最终答案。为点即可得到系数表示。再将得到的系数

用 FFT 解决多项式乘法的过程如下:
  1. 根据乘积的次数计算对应的单位根;
  2. 将单位根和系数表示传入FFT方法,递归地得到两个多项式的点值表示;
  3. 将两个多项式的点值表示相乘,得到乘积的点值表示;
  4. 将单位根的倒数和乘积的点值表示传入FFT方法,递归地得到乘积的系数表示。

有一点我们在写最终代码的时候需要特别注意。

因为我们需要递归地调用 FFT 方法,并且 FFT 方法要求传入的系数数组正负对称,所以我们需要保证每一次传入 FFT 方法的系数数组的长度为 2 的倍数。这意味着一开始的系数数组的长度应为 2 的平方。如果原数组的长度不是 2 的平方,我们应该在它的尾部填 0,直到它的长度变成 2 的平方。

比如,第一个多项式是 f(x)=1+4x−9x2,第二个多项式是 g(x)=−7+4x−9x2+9x6,它们的次数分别是 2 和 6,所以乘积的次数会是 8。乘积的点值数组的长度会是 9。

我们希望乘积的点值数组的长度是 2 的平方,也就是 16。因为在做第二项工作的时候,我们会将乘积的点值当作系数传入 FFT 方法中,而 FFT 方法要求传入的数组经过多次递归获得长度仍然是 2 的倍数。

因此,我们需要得到 16 个乘积的点值对。这意味着我们需要 16 个 f(x) 的点值对和 16 个 f(x) 的点值对,它们相乘的结果就是乘积的 16 个点值对。所以,f(x) 的系数数组从 [1,4,-9] 变成 [1,4,-9,0,0,0,0,0,0,0,0,0,0,0,0,0],同时 g(x) 的系数数组从长度 7 的 [-7,4,-9,0,0,0,9] 变成长度 16 的 [-7,4,-9,0,0,0,9,0,0,0,0,0,0,0,0,0]。

代码如下:
from cmath import pi
from cmath import exp
import numpy as np
#输出系数表示A所对应的点值表示
#w为A数组长度所对应的单位根
def FFT(A,w):
    length = len(A)
    if length==1:           #边界条件
        return [A[0]]         #直接返回常数项
    else:
        A1 =[]             #偶数次数项的系数
        A2 = []            #奇数次数项的系数
        for i in range(0,length//2):  #length肯定是2的倍数
            A1.append(A[2*i])
            A2.append(A[2*i+1])
        F1 = FFT(A1, w**2)
        F2 = FFT(A2,w**2)       #计算A1与A2的点值表示
        x=1
        #默认values第i坐标上的是第i个单位根对应的y值
        values = [None for _ in range(length)]
        #通过A1和A2的点值表示计算A的点值表示
        for i in range(0,length//2):
            values[i]=F1 [i] + x*F2[i]
            values[i+length//2]= F1[i] - x*F2[i]
            x=x*w
        return values

#输出多项式A、B的乘积
def solver(A,B):
    #填充0,将A、B的长度扩展到2的次数
    length = len(A)+len(B)-1
    n = 1
    while 2**n < length:
        n+=1
    length = 2**n
    A.extend([0]*(length-len(A)))
    B.extend([0]*(length-len(B)))
    w = exp(2*pi*1j/length)      #n次的第一个单位根
    #通过A、B的点值表示,计算A×B的点值表示
    A_values = FFT(A,w)
    B_values = FFT(B,w)
    AB_values = [A_values[i]*B_values[i] for i in range(length)]
    #将点值表示变成系数表示
    result = [round((x/length).real) for x in FFT(AB_values,w**-1)]
    while result[-1] == 0:     #将result尾部不必要的0删除
        del result[-1]
    print(result)
solver([1,5,3,2],[10,3,0,0,0,1])
输出:

[10,53,45,29,6,1,5,3,2]