function ab=sr_jacobi(dig,N,a,b)
%SR_JACOBI Recurrence coefficients for monic Jacobi polynomials.
%   AB=SR_JACOBI(DIG,N,A,B) generates in DIG-digit arithmetic the 
%   Nx2 array AB of the first N recurrence coefficients for monic
%   Jacobi polynomials orthogonal on [-1,1] relative to the weight
%   function w(x)=(1-x)^A*(1+x)^B. The call AB=SR_JACOBI(DIG,N,A)
%   is the same as AB=SR_JACOBI(DIG,N,A,A) and AB=SR_JACOBI(DIG,N) 
%   the same as AB=SR_JACOBI(DIG,N,0,0).
%
%   Supplied by Dirk Laurie, 6-22-1998; edited by Walter
%   Gautschi, 4-4-2002.
%   Edited by Walter Leopardi 10-22-2006.

digits(dig);
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
sa=vpa(a); sb=vpa(b);
nu=(sb-sa)/(sa+sb+2);
mu=2^(sa+sb+1)*gamma(sa+1)*gamma(sb+1)/gamma(sa+sb+2);
if N==1, ab=[nu mu]; return, end 
N=N-1; n=1:N; nab=2*n+sa+sb;
A=[nu (sb^2-sa^2)*ones(1,N)./(nab.*(nab+2))];
n=2:N; nab=nab(n);
B1=4*(sa+1)*(sb+1)/((sa+sb+2)^2*(sa+sb+3));
B=4*(n+sa).*(n+sb).*n.*(n+sa+sb)./((nab.^2).*(nab+1).*(nab-1));
ab=[A' [mu; B1; B']];
