Fast and accurate computation of classical Gaussian quadratures

Algorithms for computing the classical Gaussian quadrature rules (Gauss--Jacobi, Gauss--Laguerre, and Gauss--Hermite) are presented, based on globally convergent fourth-order iterative methods combined with asymptotic approximations, which are applied in complementary regions of the parameter space. This approach yields methods that improve upon existing algorithms in speed, accuracy, and computational range. The MATLAB algorithm for Gauss--Jacobi is faster than previous methods and lifts the upper restrictions on the parameters imposed by those methods ($α,β\le 5$); for example, for degrees up to $10^6$ all nodes and weights can be computed within the underflow limit for $-1<α,β\le 30$, and the computable range of parameters is much larger for smaller degrees, limited only by intrinsic overflow/underflow constraints. For the particular case of Gauss--Legendre quadrature ($α=β=0$), a specific asymptotic approach is considered, which yields the most efficient MATLAB implementation available so far. The Gauss--Laguerre and Gauss--Hermite algorithms incorporate subsampling, and scaling is also available in order to extend the computational range. Gauss--Radau and Gauss--Lobatto variants are also considered, along with the computation of the associated barycentric weights. Additionally, arbitrary-precision algorithms (in Maple) are offered for the symmetric cases (Gauss--Gegenbauer and Gauss--Hermite), which can be used to compute thousands of nodes with hundreds of digits in a matter of seconds.