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