I've been scanning across the web, and haven't found a good method to compute the Gauss Legendre abscissas and weights $\{ x\_j, w^j \} \_{j=1}^N$ for large $N\in\mathbb{N}$. My question is how to do it, and why should it work?


*To those who need some background:*


The goal is to approximate an integral by a discrete interpolating sum:


$$\int\limits\_{(-1,1)} f(x) dx = \sum\limits\_{j=-N}^N f(x\_j)w^j $$


The question is, how to choose $\{ x\_j, w^j \} \_j$ appropriately.The Gauss Legendre quadrature tells you (for good reasons) to choose $x\_j $ to be the roots of the $n$-th Legendre polynomial.


**Problem** : A straightforward computation of the Legendre polynomial for high $N$ is highly unstable, as it involves "big" coefficients of alternating signs.


**EDIT**: 
I've found a very simple code that computes weight and abscissas using eigenvalues of a symmetric matrix, but doesn't seem use Golub Welsch. The matrix is


$$\forall 1\leq j\leq N-1 \, A\_{j,j+1} = A\_{j+1,j} = \frac{j}{\sqrt{4j^2 -1}}$$ with all other entries are zero. The discussion about it was [split to another post.](https://mathoverflow.net/questions/236087/symmetric-matrix-formula-for-gaus-legendre-quadrature)



Computing Gauss Legendre quadrature for large $N$
There are asymptotic methods that essentially give you $N$ nodes and weights in $O(N)$ time if the precision is assumed to be fixed (e.g. at double precision).


See Nicholas Hale and Alex Townsend, "Fast and Accurate Computation of Gauss-Legendre and Gauss-Jacobi Quadrature Nodes and Weights", SIAM J. Sci. Comput., 35(2) (a PDF is available at <http://eprints.maths.ox.ac.uk/1629/1/finalOR79.pdf>, [Wayback Machine](https://web.archive.org/web/20181102000941/http://eprints.maths.ox.ac.uk/1629/1/finalOR79.pdf)).


They claim that their algorithm achieves double precision accuracy for $N \ge 100$.


For $N < 100$, you may as well precompute a big table with perfect accuracy using a computer algebra system or arbitrary precision library of your choice (or look up tables that others have published).


As to computing Legendre polynomials in a numerically stable way, use the three-term recurrence $(n+1) P\_{n+1}(x) = (2n+1) x P\_n(x) - n P\_{n-1}(x)$ to evaluate $P(x)$ directly instead of computing the coefficients of the polynomial and using Horner's rule (similarly for $P'(x)$).


**Update (2019)**: the issue of efficient computation for large N *with arbitrary precision* (and also with rigorous error bounds) is addressed in new work by myself and Marc Mezzarobba: SIAM Journal on Scientific Computing, 2018, Vol. 40, No. 6 : pp. C726-C747
<https://doi.org/10.1137/18M1170133> (<https://arxiv.org/abs/1802.03948>).


