% cyclotomicpolynomial_deconv.m
# recursive function to obtain cyclotomic polynomials
function retval = cyclotomicpolynomial_deconv (n)
  persistent keys = [];
  persistent values = [];
  
  if (~isscalar (n) || ~isreal(n) || n <= 0 || ceil(n) ~= floor(n))
    error ("cyclotomicpolynomial_deconv: 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);

    ## divide by each factor polynomial as we go    
    for i = d(end:-1:1)  ## dividing by bigger factors first is faster empirically
      if (any(keys==i))
        r = values(keys==i,:);
      else
        r = cyclotomicpolynomial_deconv (i);
      end
      r = r(find(r,1):end);
      retval = deconv (retval, r);
      retval = retval(find(retval,1):end);
    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

