gaussq
c accuracy
c
c the routine was tested up to n = 512 for legendre quadrature,
c up to n = 136 for hermite, up to n = 68 for laguerre, and up
c to n = 10 or 20 in other cases. in all but two instances,
c comparison with tables in ref. 3 showed 12 or more significant
c digits of accuracy. the two exceptions were the weights for
c hermite and laguerre quadrature, where underflow caused some
c very small weights to be set to zero. this is, of course,
c completely harmless.
c
c method
c
c the coefficients of the three-term recurrence relation
c for the corresponding set of orthogonal polynomials are
c used to form a symmetric tridiagonal matrix, whose
c eigenvalues (determined by the implicit ql-method with
c shifts) are just the desired nodes. the first components of
c the orthonormalized eigenvectors, when properly scaled,
c yield the weights. this technique is much faster than using a
c root-finder to locate the zeroes of the orthogonal polynomial.
c for further details, see ref. 1. ref. 2 contains details of
c gauss-radau and gauss-lobatto quadrature only.
c