function A=invtoep(c,r)
% invtoep(c,r) == inv(toeplitz(c,r))

% Dirk Laurie 2000-08-11
% See notes to solvtoep.

if nargin==1, r=c; end
[y,w]=solvtoep(c,r,[]);
n=length(c); A=zeros(n);
A(:,n)=y; A(n,:)=w'; d=y(n);
i=1; ni=n-i;
while i<=ni,
  j=i:ni;
  A(i,j)=A(n+1-j,n+1-i)';
  A(j,i)=A(n+1-i,n+1-j)';
  if j==i, return; end
  j=i:ni-1;
  A(n-j,n-i)=A(i,j)'-(w(j)*y(i)-w(n-i)*y(n-j))/d;
  A(n-i,n-j)=A(j,i)'-(y(j)*w(i)-y(n-i)*w(n-j))'/d;
  i=i+1; ni=ni-1;
end

