vector polymul( vector a, vector b)
/*
   multiply two monic polynomials specified by their coefficients

   x^(m+n) + c1 x^(m+n-1) + ... + cm+n = 

     (x^m + a1 x^(m-1) + ... + am) * (x^n + b1 x^(n-1) + ... + bn)
*/
{
  int i, j, k, m = prod(sizeof(a)), n = prod(sizeof(b)), mn = m+n;

  vector c = zeros(mn,1);

  for( j = 1; j <= n; ++j) /* x^m bj x^(n-j) */
  {
    c(j) = c(j) + b(j);
  }

  for( i = 1; i <= m; ++i) /* x^n ai x^(m-i) */
  {
    c(i) = c(i) + a(i);
  }

  for( i = 1; i <= m; ++i) /* ai x^(m-i) bj x^(n-j) */
  {
    for( j = 1; j <= n; ++j)
    {
      k = i + j;
      c(k) = c(k) + a(i)*b(j);
    }
  }

  return c;
}
