I'm trying to implement fast polynomial division using Fast Fourier Transform (fft). Here is what I have 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 same length, but if length is different then the result is wrong (I benchmark against RosettaCode's 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 weird is that it seems it's very close, but still a bit off. What did I do wrong? In other words: how to generalize fast polynomial division (using FFT) to vectors of different sizes. Also bonus if you can tell me how to compute the division quotient (currently I only have the remainder). Continue reading...