
"""
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()
