function xw=sgauss(dig,N,ab)
%SGAUSS Gauss quadrature rule.
%   XW=SGAUSS(DIG,N,AB) generates in DIG-digit arithmetic 
%   the Nx2 array XW of Gauss quadrature nodes and weights
%   for a given weight function w. The nodes, in increasing
%   order, are placed into the first column of XW and the
%   corresponding weights into the second column. The weight
%   function w is specified by the Nx2 input array AB of
%   recurrence coefficients for the polynomials orthogonal
%   with respect to the weight function w.

digits(dig)
N0=size(ab,1); if N0<N, error('input array ab too short'), end
sab=vpa(ab);
J=zeros(N); J=vpa(J);
for n=1:N, J(n,n)=sab(n,1); end
for n=2:N
  J(n,n-1)=sqrt(sab(n,2));
  J(n-1,n)=J(n,n-1);
end
[V,D]=eig(J);
D=diag(D);
[v,I]=sort(double(D));
D=D(I);
V=V(:,I);
xw=[D vpa(sab(1,2))*V(1,:)'.^2];
