For a given positive measure on a fixed domain, Gaussian quadrature routines can be defined via their property of integrating an arbitrary polynomial of degree $2n+1$ exactly using only $n+1$ quadrature nodes. In the special case of Gauss--Jacobi quadrature, this means that $$\int_{-1}^1 (1+x)^\alpha(1-x)^\beta f(x) dx = \sum_{j=0}^{n} w_j f(x_j), \quad \alpha, \beta > -1, $$ whenever $f(x)$ is a polynomial of degree at most $2n+1$. When $f$ is not a polynomial, but a function analytic in a neighbourhood of $[-1,1]$, the above is not an equality but an approximation that converges exponentially fast as $n$ is increased.
An undergraduate mathematician taking a numerical analysis course could tell you that the nodes $x_j$ are roots of the Jacobi polynomial $P^{\alpha,\beta}_{n+1}(x)$, the degree $n+1$ polynomial orthogonal with respect to the given weight, and that the quadrature weight at each node is related to the derivative $P'^{\alpha,\beta}_{n+1}(x_j)$. However, these values are not generally known in closed form, and we must compute them numerically... but how?
Traditional approaches involve applying the recurrence relation satisfied by the orthogonal polynomials, or solving the Jacobi matrix eigenvalue problem in the algorithm of Golub and Welsch, but these methods are inherently limited by a minimal complexity of $O(n^2)$. The current state-of-the-art is the $O(n)$ algorithm developed by Glasier, Liu, and Rokhlin, which hops from root to root using a Newton iteration evaluated with a Taylor series defined by the ODE satisfied by $P^{\alpha,\beta}_{n+1}$.
We propose an alternative approach, whereby the function and derivative evaluations required in the Newton iteration are computed independently in $O(1)$ operations per point using certain well-known asymptotic expansions. We shall review these expansions, outline the new algorithm, and demonstrate improvements in both accuracy and efficiency.