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. 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
n=1:N; ab(1,2)=prod(1+(a+b+1)./n);
if a+b==0
  n=(1:N+1)';
  ab(n,1)=((2*n+a+b-1)*N+(b-a)*n+a)./(2*(2*n-1));
  n=(1:N)';
  ab(n+1,2)=.25*((N+1)^2)*(1+a./n).*(1+b./n).*(1-(n./(N+1)).^2)...
    ./(4-(1./n).^2);
elseif a+b+1==0
  n=(1:N+1)';
  ab(n,1)=((2*(n-1).^2+b)*N+(2*b+1)*(n-1).^2)./(4*(n-1).^2-1);
  n=(1:N)';
  ab(n+1,2)=.25*((N+1)^2)*(1+a./n).*(1+b./n).*(1-n./(N+1)).*...
    (1+(n-1)./(N+1))./(4-(1./n).^2);
else
  n=(1:N+1)';
  ab(n,1)=((n+a+b).*(n+a).*(N-n+1)./(2*n+a+b)+(n-1).*(n+b-1).*...
    (N+n+a+b)./(2*n+a+b-2))./(2*n+a+b-1);
  n=(1:N)';
  ab(n+1,2)=((N+1)^2)*(1+a./n).*(1+b./n).*(1+(a+b)./n).*...
    (1-n./(N+1)).*(1+(n+a+b)./(N+1))./(((2+(a+b)./n).^2).*...
    ((2+(a+b)./n).^2-(1./n).^2));
end
