function [beta,r,l1] = spefac2(right)
%
% function [beta,r,l1] = spefac2(right)      Mikael Sternad 1992, rev 950517
%
% Given the symmetric polynomial right(z,z^-1),
% the sequence beta(1)..beta(n) and the real r is computed such that
% r*[z**n+beta(1)*z**(n-1)+..+beta(n)]
%  *[z**(-n)+beta(1)*z**(-n+1)+..beta(n)] = right(z,z^-1)
% 'right' is assumed to be positive definite on the unit circle
%
% Method: calculates roots of right-hand side, and creates beta from
% the roots in |z|<1.
%
% The polynomial may have complex coefficients and roots close to |z|=1.
%
% Inaccurate results may occur for clustered roots.
%
% l1 is the 1-norm of the error rbetabeta* - right.

lr = length(right);
nr = fix(lr-1)/2;
if right(nr+1)<max(abs(right)),
  error('right-hand side not pos.definite')
  return
end
totalroots =roots(right);
betaroots = [];
for k=1:lr-1
  if abs(totalroots(k))<1,
    if length(betaroots)<nr,
      betaroots = [betaroots ; totalroots(k)];
    end
  end
end

if length(betaroots) < nr,
  error('Too few roots inside the unit circle!'),return, end
beta = poly(betaroots);
betabeta = abstar(beta,beta);
r = right(nr+1)/betabeta(nr+1);
beta = beta/beta(1);
if sum(abs(imag(beta)))<1e-4, beta = real(beta); end

fel  = r(1)*conv(beta,conj(beta(length(beta):-1:1))) - right;
l1 = norm(fel,1);