% cyclotomicpolynomial_conv.m
# recursive function to obtain cyclotomic polynomials
function retval = cyclotomicpolynomial_conv (n)
  persistent keys = [];
  persistent values = [];

  if (~isscalar (n) || ~isreal(n) || n <= 0 || ceil(n) ~= floor(n))
    error ("cyclotomicpolynomial_conv: N = %f is not a positive integer!", n)
  elseif (any(keys==n))     ## return memoized result
    retval = values(keys==n,:);
    return
  elseif (n == 1)  ## return (x - 1)
    retval = [+1 -1];
  elseif (isprime(n))  ## return 1 + x + ... + x^n-1
    retval = ones(1,n);
  else  ## composite n ==> need to factorize
    # start with x^n - 1
    retval = zeros(1, n+1);
    retval(1) = +1;
    retval(n+1) = -1;

    # divisors
    d = 1:(n/2);
    q = round(n ./ d);
    d = d(d.*q == n);

    # multiply all the factor polynomials together, then divide once
    denom = 1;
    for i = d
      if (any(keys==i))
        r = values(keys==i,:);
      else
        r = cyclotomicpolynomial_conv (i);
      end
      r = r(find(r,1):end);
      denom = conv2 (denom, r);
      denom = denom(find(denom,1):end);
    end
    retval = deconv (retval, denom);
    retval = retval(find(retval,1):end);
  end

  ## add to memoization table
  keys(end+1) = n;
  if (length(retval) > size(values,2))
    values = [zeros(size(values,1), length(retval)-size(values,2)), values];
  end
  values(end+1, (end+1-length(retval)):end) = retval;
end

