फास्ट बहुपद विभाजन के लिए एफएफटी डिवीजन

Question

मैं फास्ट फूरियर ट्रांसफॉर्म (एफएफटी) का उपयोग कर तेजी से बहुपद विभाजन को कार्यान्वित करने की कोशिश कर रहा हूं।

यहाँ मैं अब तक क्या मिला है है:

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'sextended_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]) 

अजीब बात यह है कि ऐसा लगता है कि यह बहुत करीब है, लेकिन अभी भी थोड़ा सा है। मैंने गलत क्या किया? दूसरे शब्दों में: विभिन्न आकार के वैक्टरों को तेज़ बहुपद विभाजन (एफएफटी का उपयोग करके) को सामान्यीकृत करने के लिए कैसे करें।

इसके अलावा बोनस यदि आप मुझे बताओ कि गणना करने के लिए विभाजन भागफल कर सकते हैं (वर्तमान में मैं केवल शेष है)।

Answer 1

यहां इन lecture notes में पाए गए एक तेज बहुपद विभाजन एल्गोरिदम का प्रत्यक्ष कार्यान्वयन है।

विभाजन डिवीजन के पारस्परिक के साथ लाभांश के तेज़/एफएफटी गुणा पर आधारित है। नीचे मेरा कार्यान्वयन 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) समारोह बहुआयामी पद u और v के लिए एक (quotient, remainder) टपल (पायथन के मानक 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]) 

यानी:

उदाहरण के लिए

• (x^2 - 1)/(x - 1) == x + 1
• (2x^3 - 5x^2 + 10x + 8)/(x^2 + 2x -3) == 3x - 11, शेष के साथ 41x - 25
• आदि (पिछले दो उदाहरण तुम्हारा कर रहे हैं।)