I have written a package FastGauss.jl available here: https://github.com/ajt60gaibb/FastGauss.jl
to compute Gauss quadrature rules to 16-digit precision (so far Legendre, Jacobi, Lobatto, Radau), aiming to be the fastest implementation in Julia. For example, this is how long it takes to compute 1,000,000 Gauss-Legendre nodes and weights: tic(), GaussLegendre( 1000000 ); toc() elapsed time: 0.336489122 seconds In my comparisons, GaussLegendre() is faster than Base.gauss() for n>60 (nothing in it for small n), but my implementation right now does not allow for higher precision. I couldn't find a Gauss-Jacobi code in Base. I am a Julia beginner (only been learning for 2 weeks) so I am assuming the code can be improved in a million and one ways. Please tell me if I've done something that Julia does not like. I am not sure if it is appropriate to make this an official package.