Fft for fast polynomial division

Question

I am trying to implement fast polynomial division using the fast Fourier transform (fft).

Here is what I got so far:

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]] 

This works well for polynomials of the same length, but if the length is different, then the result is incorrect (I compare with the RosettaCode extended_synthetic_division() function):

 # 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]) 

What is strange is that it seems very close, but still a bit. What have I done wrong? In other words: how to generalize fast polynomial division (using FFT) into vectors of different sizes .

Also a bonus if you can tell me how to calculate quotient division (currently I have a remainder ).

Answer 1

Here's a direct implementation of the fast polynomial division algorithm found in these lectures .

The separation is based on the fast / FFT multiplication of the dividend with the divisor by turnover. My implementation below strictly follows an algorithm that, as it turned out, has the time complexity of O(n*log(n)) (for polynomials with degrees of the same order of magnitude), but it is written with an emphasis on readability rather than efficiency.

 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([ab 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 

The function poly_divmod(u, v) returns a tuple (quotient, remainder) for the polynomials u and v (for example, the Python standard divmod for numbers).

For instance:

 >>> 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
• (2x^3 - 5x^2 + 10x + 8) / (x^2 + 2x -3) == 3x - 11 , with the remainder 41x - 25
• etc .. (The last two examples are yours.)