python– 用于快速多项式除法的fft除法

概述我正在尝试使用快速傅立叶变换(fft)实现快速多项式除法.这是我到目前为止所得到的：from numpy.fft import fft, ifft def fft_div(C1, C2): # fft expects right-most for significant coefficients C1 = C1[::-1] C2 =

我正在尝试使用快速傅立叶变换(fft)实现快速多项式除法.

这是我到目前为止所得到的：

@H_404_8@from numpy.fft import fft,ifftdef fft_div(C1,C2):    # fft expects right-most for significant coefficIEnts    C1 = C1[::-1]    C2 = C2[::-1]    d = len(C1)+len(C2)-1    c1 = fft(List(C1) + [0] * (d-len(C1)))    c2 = fft(List(C2) + [0] * (d-len(C2)))    res = List(ifft(c1-c2)[:d].real)    # Reorder back to left-most and round to integer    return [int(round(x)) for x in res[::-1]]

这适用于相同长度的多项式,但如果长度不同则结果是错误的(我对RosettaCode’s extended_synthetic_division()函数进行基准测试)：

@H_404_8@# Most signficant coefficIEnt is leftN = [1,-11,-22,1]D = [1,-3,1,2]# OK case,same length for both polynomialsfft_div(N,D)>> [0,-8,-23,-1]extended_synthetic_division(N,D)>> ([1],[-8,-1])# NOT OK case,D is longer than N (also happens if shorter)D = [1,2,20]fft_div(N,-1,4,-24,-19]extended_synthetic_division(N,D)>> ([],[1,1])

奇怪的是它看起来非常接近,但仍然有点偏离.我做错了什么？换句话说：如何将快速多项式除法(使用FFT)推广到不同大小的向量.

如果你能告诉我如何计算除法商(目前我只有余数),也可以获得奖励.

最佳答案这是在这些lecture notes中找到的快速多项式除法算法的直接实现.

除法基于除数的快速/ FFT乘法与除数的倒数.我在下面的实现严格遵循经证明具有O(n * log(n))时间复杂度的算法(对于具有相同数量级的多项式),但是它的重点在于可读性而非效率.

@H_404_8@from math import ceil,logfrom numpy.fft import fft,ifftdef poly_deg(p):    return len(p) - 1def poly_scale(p,n):    """Multiply polynomial ``p(x)`` with ``x^n``.    If n is negative,poly ``p(x)`` is divIDed with ``x^n``,and remainder is    discarded (truncated division).    """    if n >= 0:        return List(p) + [0] * n    else:        return List(p)[:n]def poly_scalar_mul(a,p):    """Multiply polynomial ``p(x)`` with scalar (constant) ``a``."""    return [a*pi for pi in p]def poly_extend(p,d):    """Extend List ``p`` representing a polynomial ``p(x)`` to    match polynomials of degree ``d-1``.    """    return [0] * (d-len(p)) + List(p)def poly_norm(p):    """normalize the polynomial ``p(x)`` to have a non-zero most significant    coefficIEnt.    """    for i,a in enumerate(p):        if a != 0:            return p[i:]    return []def poly_add(u,v):    """Add polynomials ``u(x)`` and ``v(x)``."""    d = max(len(u),len(v))    return [a+b for a,b in zip(poly_extend(u,d),poly_extend(v,d))]def poly_sub(u,v):    """Subtract polynomials ``u(x)`` and ``v(x)``."""    d = max(len(u),len(v))    return poly_norm([a-b for a,d))])def poly_mul(u,v):    """Multiply polynomials ``u(x)`` and ``v(x)`` with FFT."""    if not u or not v:        return []    d = poly_deg(u) + poly_deg(v) + 1    U = fft(poly_extend(u,d)[::-1])    V = fft(poly_extend(v,d)[::-1])    res = List(ifft(U*V).real)    return [int(round(x)) for x in res[::-1]]def poly_recip(p):    """Calculate the reciprocal of polynomial ``p(x)`` with degree ``k-1``,defined as: ``x^(2k-2) / p(x)``,where ``k`` is a power of 2.    """    k = poly_deg(p) + 1    assert k>0 and p[0] != 0 and 2**round(log(k,2)) == k    if k == 1:        return [1 / p[0]]    q = poly_recip(p[:k/2])    r = poly_sub(poly_scale(poly_scalar_mul(2,q),3*k/2-2),poly_mul(poly_mul(q,p))    return poly_scale(r,-k+2)def poly_divmod(u,v):    """Fast polynomial division ``u(x)`` / ``v(x)`` of polynomials with degrees    m and n. Time complexity is ``O(n*log(n))`` if ``m`` is of the same order    as ``n``.    """    if not u or not v:        return []    m = poly_deg(u)    n = poly_deg(v)    # ensure deg(v) is one less than some power of 2    # by extending v -> ve,u -> ue (mult by x^nd)    nd = int(2**ceil(log(n+1,2))) - 1 - n    ue = poly_scale(u,nd)    ve = poly_scale(v,nd)    me = m + nd    ne = n + nd    s = poly_recip(ve)    q = poly_scale(poly_mul(ue,s),-2*ne)    # handle the case when m>2n    if me > 2*ne:        # t = x^2n - s*v        t = poly_sub(poly_scale([1],2*ne),poly_mul(s,ve))        q2,r2 = poly_divmod(poly_scale(poly_mul(ue,t),-2*ne),ve)        q = poly_add(q,q2)    # remainder,r = u - v*q    r = poly_sub(u,poly_mul(v,q))    return q,r

poly_divmod(u,v)函数返回多项式u和v的(商,余数)元组(如Python的标准divmod数字).

例如：

@H_404_8@>>> print poly_divmod([1,-1],-1])([1,1],[])>>> print poly_divmod([3,-5,10,8],-3])([3,-11],[41,-25])>>> print poly_divmod([1,2])([1],-1])>>> print poly_divmod([1,20])([],1])

即：

>(x ^ 2 – 1)/(x – 1)== x 1
>(2x ^ 3 – 5x ^ 2 10x 8)/(x ^ 2 2x -3)== 3x – 11,其余为41x – 25
>等等(最后两个例子是你的.) 总结

以上是内存溢出为你收集整理的python – 用于快速多项式除法的fft除法全部内容，希望文章能够帮你解决python – 用于快速多项式除法的fft除法所遇到的程序开发问题。

如果觉得内存溢出网站内容还不错，欢迎将内存溢出网站推荐给程序员好友。