FFT在Python中使用Instole numpy.fft中的多项式乘法

Question

我想在Python中快速乘以两个多项式。由于我的多项式很大(> 100000(元素，因此我必须乘以很多。下面，您会找到我的方法

from numpy.random import seed, randint
from numpy import polymul, pad 
from numpy.fft import fft, ifft
from timeit import default_timer as timer
length=100
def test_mul(arr_a,arr_b):  #inbuilt python multiplication
    c=polymul(arr_a,arr_b)
    return c    
def sb_mul(arr_a,arr_b):    #my schoolbook multiplication
    c=[0]*(len(arr_a) + len(arr_b) - 1 )
    for i in range( len(arr_a) ):
        for j in range( len(arr_b) ):
            k=i+j
            c[k]=c[k]+arr_a[i]*arr_b[j]
    return c    

def fft_test(arr_a,arr_b):  #fft based polynomial multuplication
    arr_a1=pad(arr_a,(0,length),'constant')
    arr_b1=pad(arr_b,(0,length),'constant')
    a_f=fft(arr_a1)
    b_f=fft(arr_b1)
    c_f=[0]*(2*length)
    for i in range( len(a_f) ):
        c_f[i]=a_f[i]*b_f[i]
    return c_f

if __name__ == '__main__':
    seed(int(timer()))
    random=1
    if(random==1):
        x=randint(1,1000,length)
        y=randint(1,1000,length)
    else:
        x=[1]*length
        y=[1]*length
    start=timer()
    res=test_mul(x,y)
    end=timer()
    print("time for built in pol_mul", end-start)
    start=timer()
    res1=sb_mul(x,y)
    end=timer()
    print("time for schoolbook mult", end-start)
    res2=fft_test(x,y)
    print(res2)
    #########check############
    if( len(res)!=len(res1) ):
        print("ERROR");
    for i in range( len(res) ):
        if( res[i]!=res1[i] ):
            print("ERROR at pos ",i,"res[i]:",res[i],"res1[i]:",res1[i])

现在，这是我的详细方法，1.首先，我尝试了自己的复杂性o(n^2(的幼稚实施。但是，正如您可能期望的那样，事实证明它非常慢。

2. 第二，我在numpy库中了解了polymul。此功能比上一个功能快很多。但是我意识到这也是O(n^2(的复杂性。您可以看到，如果增加了长度k的时间，则时间增加了2次。

3. 我的第三种方法是使用内置的FFT函数尝试基于FFT的乘法。我遵循此处描述的众所周知的方法，但我无法使它起作用。

现在我的问题是

1. 在基于FFT的方法中我在哪里出错？你能告诉我如何修复它？

2. 我的观察是polymul功能具有O(n^2(的复杂性正确吗？

请让我知道您是否有任何疑问。预先感谢。

Answer 1

1. 我在基于FFT的方法中出错了哪里？你能告诉我如何解决吗？

主要问题是，在基于FFT的方法中，您应该在乘法后进行反变形，但是代码中缺少该步骤。在此缺失的步骤中，您的代码应该看起来如下：

def fft_test(arr_a,arr_b):  #fft based polynomial multiplication
    arr_a1=pad(arr_a,(0,length),'constant')
    arr_b1=pad(arr_b,(0,length),'constant')
    a_f=fft(arr_a1)
    b_f=fft(arr_b1)
    c_f=[0]*(2*length)
    for i in range( len(a_f) ):
        c_f[i]=a_f[i]*b_f[i]
    return ifft(c_f)

请注意，也可能有一些改进的机会：

• 可以通过将所需的FFT长度传递为第二个参数(例如a_f = fft(arr_a, length)(
来直接处理零填充物
• for循环中的系数乘法可以由numpy.multiply直接处理。
• 如果多项式系数是真实价值的，则可以使用numpy.fft.rfft和numpy.fft.irfft(而不是numpy.fft.fft和numpy.fft.ifft(进行一些额外的性能提升。

因此，实现现值输入的实现可能看起来像：

from numpy.fft import rfft, irfft
def fftrealpolymul(arr_a, arr_b):  #fft based real-valued polynomial multiplication
    L = len(arr_a) + len(arr_b)
    a_f = rfft(arr_a, L)
    b_f = rfft(arr_b, L)
    return irfft(a_f * b_f)

2. 我的观察是polymul功能具有O(N ² (的复杂性正确吗？

似乎也是我观察到的性能，并匹配我的numpy安装中的可用代码(版本1.15.4，并且在最近的1.16.1版本中似乎没有任何变化(。