时间戳

首页 > 解决方案 > 使用内置 Numpy.fft 在 Python 中进行 FFT 多项式乘法

python - 使用内置 Numpy.fft 在 Python 中进行 FFT 多项式乘法

问题描述

我想在 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) 的 Schoolbook 的简单实现。但正如你所料，它变得非常缓慢。

  1. 其次，我是polymul在 Numpy 库中认识的。这个功能比前一个快很多。但我意识到这也是一个 O(n^2) 复杂度。你可以看到，如果你增加长度 k 时间增加了 k^2 倍。

  2. 我的第三种方法是使用内置的 FFT 函数尝试基于 FFT 的乘法。我遵循了这里也描述的众所周知的方法，但我无法让它工作。

现在我的问题是，

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

  2. 我观察到polymul函数具有 O(n^2) 复杂度是否正确？

如果您有任何问题，请告诉我。提前致谢。

标签： pythonfftmultiplicationpolynomials

解决方案

  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.rfftand numpy.fft.irfft（而不是numpy.fft.fftand 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)
  1. 我对polymul函数具有 O(n 2 ) 复杂度的观察是否正确？

这似乎也是我正在观察的性能，并且与我的 numpy 安装中的可用代码相匹配（版本 1.15.4，在最近的 1.16.1 版本中该部分似乎没有任何变化）。

推荐阅读