FFT除法用于快速多项式除法

Question

我正尝试使用快速傅里叶变换（fft）实现快速多项式除法。

这是我迄今为止得到的东西：

from numpy.fft import fft, ifft
def 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 的 extended_synthetic_division() 函数进行了基准测试）：

# Most signficant coefficient is left
N = [1, -11, 0, -22, 1]
D = [1, -3, 0, 1, 2]
# OK case, same length for both polynomials
fft_div(N, D)
>> [0, 0, 0, 0, 0, -8, 0, -23, -1]
extended_synthetic_division(N, D)
>> ([1], [-8, 0, -23, -1])

# NOT OK case, D is longer than N (also happens if shorter)
D = [1, -3, 0, 1, 2, 20]
fft_div(N, D)
>> [0, 0, 0, 0, -1, 4, -11, -1, -24, -19]
extended_synthetic_division(N, D)
>> ([], [1, -11, 0, -22, 1])

奇怪的是它似乎很接近，但还是有点偏差。我做错了什么？换句话说：如何将快速多项式除法（使用FFT）推广到不同大小的向量。

如果您能告诉我如何计算商，那就更好了（目前只有余数）。

- gaborous

---

你从哪里得到这个算法的？ - Alex

以下链接中第33页的算法看起来不像是你的：http://www.diag.uniroma1.it/sankowski/lecture4.pdf - Alex

这是反卷积的一个简单常见实现：ifft(fft(A).*fft(B))，因为我读到多项式除法等于反卷积。但实际上，在你提供的链接中这是不同的算法，我会研究一下（除非你能编写 Python 代码来实现它，那你就可以得到赏金！:D） - gaborous

Answer 1

这是一个快速多项式除法算法的直接实现，该算法在这些讲义中找到。

该除法基于将被除数与除数的倒数进行快速/FFT乘法。我下面的实现严格遵循已证明具有O(n*log(n))时间复杂度（对于阶数相同数量级的多项式）的算法，但它更注重可读性而非效率。

from math import ceil, log
from numpy.fft import fft, ifft

def poly_deg(p):
    return len(p) - 1


def 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,b in zip(poly_extend(u, d), poly_extend(v, 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, 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)函数返回多项式和的商和余数的元组（类似于Python中数字的标准divmod函数）。

例如：

>>> print poly_divmod([1,0,-1], [1,-1])
([1, 1], [])
>>> print poly_divmod([3,-5,10,8], [1,2,-3])
([3, -11], [41, -25])
>>> print poly_divmod([1, -11, 0, -22, 1], [1, -3, 0, 1, 2])
([1], [-8, 0, -23, -1])
>>> print poly_divmod([1, -11, 0, -22, 1], [1, -3, 0, 1, 2, 20])
([], [1, -11, 0, -22, 1])

I.e:

• (x^2 - 1) / (x - 1) == x + 1（其中x不等于1）
• (2x^3 - 5x^2 + 10x + 8) / (x^2 + 2x -3) == 3x - 11，余数为41x - 25
• 等等。（最后两个例子是您的。）

- randomir

---

非常欢迎。顺便说一下，我相信这里采用的反转多项式的新方法链接1和链接2可以导致更简单的实现，但我还没有尝试过。 - randomir