function ab=sr_hahn_ops(dig,N,a,b)
%SR_HAHN_OPS Recurrence coefficients for the monic Hahn 
%orthogonal polynomials.
%   AB=SR_HAHN_OPS(DIG,N,A,B) generates in Dig-digit arithmetic 
%   the (N+1)x2 array AB of the N+1 recurrence coefficients for
%   the weight function w(k)=(A+k choose k)*(B+N-k choose N-k),
%   k=0:N, A>=0, B>=0. The call AB=SR_HAHN_OPS(DIG,N,A) is the 
%   same as AB=SR_HAHN_OPS(DIG,N,A,A), and AB=SR_HAHN_OPS(DIG,N) 
%   the same as AB=SR_HAHN_OPS(DIG,N,0,0). which produces the
%   recurrence coefficients of the discrete Chebyshev polynomials
%   for the points 0,1,2,...,N.

digits(dig); 
if nargin<3, a=0; end; if nargin<4, b=a; end
if((N<=0)|(a<-1)|(b<-1)), error('parameter(s) out of range'), end
ab=vpa(zeros(N+1,2));
sN=vpa(N); sa=vpa(a); sb=vpa(b); s1=vpa(1); s2=vpa(2); s4=vpa(4);
k=1:N; sk=vpa(k); ab(1,2)=prod(s1+(sa+sb+s1)./sk);
if double(a+b)==0
  k=(1:N+1)'; sk=vpa(k);
  ab(k,1)=((2*sk+sa+sb-s1)*sN+(sb-sa)*sk+sa)./(s2*(s2*sk-s1));
  k=(1:N)'; sk=vpa(k);
  ab(k+1,2)=vpa(1/4)*((sN+s1)^2)*(s1+sa./sk).*(s1+sb./sk) ...
    .*(s1-(sk./(sN+s1)).^2)./(s4-(s1./sk).^2);
elseif double(a+b+1)==0
  k=(1:N+1)'; sk=vpa(k);
  ab(k,1)=((s2*(sk-s1).^2+sb)*sN+(s2*sb+s1)*(sk-1).^2)./(s4*(sk-s1).^2-s1);
  k=(1:N)'; sk=vpa(k);
  ab(k+1,2)=vpa(1/4)*((sN+s1)^2)*(s1+sa./sk).*(s1+sb./sk).*(s1-sk./(sN+s1)).* ...
    (s1+(sk-s1)./(sN+s1))./(s4-(s1./sk).^2);
else
  k=(1:N+1)'; sk=vpa(k);
  ab(k,1)=((sk+sa+sb).*(sk+sa).*(sN-sk+s1)./(s2*sk+sa+sb)+(sk-s1).*(sk+sb-s1).* ...
    (sN+sk+sa+sb)./(s2*sk+sa+sb-s2))./(s2*sk+sa+sb-s1);
  k=(1:N)'; sk=vpa(k);
  ab(k+1,2)=((sN+s1)^2)*(s1+sa./sk).*(s1+sb./sk).*(s1+(sa+sb)./sk).* ...
    (s1-sk./(sN+s1)).*(s1+(sk+sa+sb)./(sN+s1))./(((s2+(sa+sb)./sk).^2).* ...
    ((s2+(sa+sb)./sk).^2-(s1./sk).^2));
end
