function ab=sr_hahn(dig,N,a,b)
% SR_HAHN Recurrence coefficients for monic Hahn polynomials in
% variable-precision arithmetic.
%
%    AB=R_HAHN(N,A,B) generates the N+1 recurrence coefficients 
%    for the monic Hahn polynomials with parameters A and B. These
%    are orthogonal on the discrete set of N+1 points 0,1,2,...,N
%    with weights (A+k choose k)(B+N-k choose N-k), k=0,1,2,...,N.
%    The N+1 alpha-coefficients are stored in the first column,
%    the N+1 beta-coefficients in the second column, of the
%    (N+1)x2 array AB. The call AB=R_HAHN(N,A) is the same as AB=
%    R_HAHN(N,A,A), and AB=HAHN(N) the same as AB=R_HAHN(N,0,0),
%    which produces the recurrence coefficients of the discrete 
%    Chebyshev polynomials for the points 0,1,2,...,N.

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