""" Naive and FFT polynomial multiplication. Connelly Barnes 2005-11-12. Public domain. """ import FFT # Requires Numeric. def roundup_pow2(x): """ Round up to power of 2 (obfuscated and unintentionally faster :). """ while x&(x-1): x = (x|(x>>1))+1 return max(x,1) def fft_mul(coeff1, coeff2): """ FFT multiply two real polynomials (lowest power coeffs are first). """ d = roundup_pow2(len(coeff1)+len(coeff2)-1) c1 = FFT.fft(list(coeff1) + [0] * (d-len(coeff1))) c2 = FFT.fft(list(coeff2) + [0] * (d-len(coeff2))) return list(FFT.inverse_fft(c1 * c2)[:len(coeff1)+len(coeff2)-1].real) def naive_mul(coeff1, coeff2): """ Naively multiply two polynomials (lowest power coeffs are first). """ d = len(coeff1)+len(coeff2)-1 ans = [0] * d for i in range(d): for j1 in range(len(coeff1)): j2 = i-j1 if 0 <= j2 < len(coeff2): ans[i] += coeff1[j1] * coeff2[j2] return ans def eval_poly(coeff, x): """ Evaluate polynomial (lowest power coeffs are first). """ return reduce(lambda a, b: a*x+b, coeff[::-1]) def test(): """ Unit test. """ import random for i in range(1000): c1 = [random.random() for i in range(random.randrange(1,20))] c2 = [random.random() for i in range(random.randrange(1,20))] p1 = naive_mul(c1, c2) p2 = fft_mul(c1, c2) assert [abs(p1[i]-p2[i])<1e-8 for i in range(len(p1))].count(0)==0 if __name__ == '__main__': test()