function [numn,denn]=onestepr(numnp1,dennp1);

% Best rational approximation of degree n of rational function of degree n+1
% numnp1 and dennp1 contain coeff. of x^0,1,2,... of given numerator and
% denominator
[aux,degn]=size(numnp1);[aux,degd]=size(dennp1);
np1=max(degn,degd)-1;n=np1-1;
disp(' degree ');disp(np1);
xres=[];xroots=[];zres=[];zroots=[];
if degd==1, xpolpart=numnp1(degn:-1:1)/dennp1(1);
else, [xres,xroots,xpolpart]=residue(numnp1(degn:-1:1),dennp1(degd:-1:1));end;
[pold1,pold2]=size(xpolpart);if pold1==0, pold1=1;pold2=1;xpolpart(1)=0;end;
zpolpart=zeros(1,2*pold2-1);
for k=1:pold2;zpolpart=([zpolpart(2:2*pold2-1),0]+[0,zpolpart(1:2*pold2-2)])/2;
     zpolpart(pold2)=zpolpart(pold2)+xpolpart(k);
end;
zpolpart=zpolpart(1:pold2);if pold2>1, zpolpart(1:pold2-1)=2*zpolpart(1:pold2-1);end;
if degd==1,  znumnp1=zpolpart;zdennp1=zeros(1,pold2);zdennp1(pold2)=1;
else,
  zroots=xroots+sqrt(xroots.*xroots-1);
  for k=1:degd-1;if abs(zroots(k))<1, zroots(k)=1/zroots(k);end;end;
  [sr1,sr2]=size(xres);onone=ones(sr1,sr2);
  zres=4*xres./(onone - onone./(zroots.*zroots));
  zpolpart(pold2)=zpolpart(pold2)+sum(zres./zroots)/2;
  [znumnp1,zdennp1]=residue(zres,zroots,zpolpart);
end;
  xres,xroots,xpolpart,
  zres,zroots,zpolpart,
  znumnp1,zdennp1,
[aux,degzn]=size(znumnp1);[aux,degzd]=size(zdennp1);
smatr=[toeplitz([zdennp1,zeros(1,2*n+2-degzd)],zeros(1,np1)), ...
        hankel([zeros(1,2*np1-degzd),zdennp1(degzd:-1:1)],zeros(1,np1))],
nmatr=[hankel([zeros(1,2*np1-degzn),znumnp1(degzn:-1:1)],zeros(1,np1)),zeros(2*np1,np1)],
smatr\nmatr,
nmatr/smatr,
[vecp,lambd]=eig(nmatr,smatr),
[lambs,ksort]=sort(-abs(diag(lambd)));
for k=1:np1,
  ks=ksort(k);
  disp(ks);format long;disp(lambd(ks,ks));format short;
  rv=roots(vecp(1:np1,ks));
  [rv,abs(rv)],
  if k==np1, kmin=ks; zroots2=rv;lamb2=lambd(ks,ks);end;
end;
numnz=lamb2*vecp(np1+1:2*np1,kmin),
dennz=vecp(1:np1,kmin),
% in x=(z+1/z)/2  :
denn=dennz(1);if n>0,denn=poly( (zroots2+ones(n,1)./zroots2)/2 );end;
% (n(z)/d(z)+n(1/z)/d(1/z))/2 :
numn=numnz(1);if n>0,
 chebp=exp(i*pi*(0:n)/n);chebp2=ones(1,np1)./chebp;
 zval=0.5*( polyval(numnz,chebp)./polyval(dennz,chebp)    ...
              +polyval(numnz,chebp2)./polyval(dennz,chebp2));
 numn=polyfit((chebp+chebp2)/2,polyval(denn,(chebp+chebp2)/2).*zval,n);
end;
denn=denn(np1:-1:1),
numn=numn(np1:-1:1),
