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).

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).

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).

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).

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).

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)$).

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).

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).