function  Aa=mpolflipr(A,z0)
%function  Aa=mpolflipr(A,z0)
% A is an s-1 order square matrix polynomial with a root at z0.
% Aa has a root at 1/z0 instead, and satisfies Aa(z)*Aa'(1/z)=A(z)*A'(1/z).
% The algorithm follows Rozanov's Stationary Random Processes, p.46-47.
realsmall=1e-12;
[n,m,s]=size(A);
zv=z0.^[0:s-1];
Az=prodt(A,zv,3,2);
[u,d,v]=svd(Az);
if ~(abs(d(n,n))<realsmall)
   error(sprintf('p(z0)=%g~=0',d(n,n)))
else
   Aa=prodt(A,v,2,1);
   Aa=permute(Aa,[1 3 2]);
   Aa(:,n,:)=prodt(Aa(:,n,:),toeplitz([1 zeros(1,s-1)],1../zv),3,1);
   Aa(:,n,:)=(Aa(:,n,:)-conj(z0)*cat(3,zeros(n,1,1),Aa(:,n,1:end-1 )))/abs(z0);
end
% Steps below make Aa real, if possible (i.e. if Aa(:,:,1)*Aa(:,:,1)' is real.)
% These steps are never invoked if A and z0 are both real.
if norm(imag(Aa(:,:,1)))>n*realsmall
   sig=Aa(:,:,1)*Aa(:,:,1)';
   if ~(norm(imag(sig))>n*realsmall)
      sig=real(sig);
      [w,err]=chol(sig);
      if err==0
         Aa=prodt(Aa,Aa(:,:,1)'/w,2,1);
         Aa=permute(Aa,[1 3 2]);
         Aa=real(Aa);
      end
   end
end