Roots of Legendre polynomials and Gauss–Legendre nodes

arb_hypgeom_legendre_p_ui_root computes an enclosure of the k-th root x_k of the Legendre polynomial P_n and, optionally, the weight w_k = 2/((1-x_k^2) P_n'(x_k)^2), so that

is the n-point Gauss–Legendre quadrature rule, exact for polynomials of degree less than 2n. The roots are numbered in decreasing order, 1 > x_0 > x_1 > \cdots > x_{n-1} > -1, which is the order of increasing \theta for the roots of P_n(\cos\theta). Only the roots with k < \lceil n/2 \rceil need to be computed, since x_{n-1-k} = -x_k and w_{n-1-k} = w_k; the function uses this symmetry itself. An index k \ge n raises an exception. The function is the one used by acb_calc_integrate to generate its quadrature rules (Rigorous numerical integration). The method is that of Johansson and Mezzarobba [JM2018].

Algorithm#

1. Initial approximation. With \phi = \frac{4k+3}{4n+2}\pi (that is, (4k'-1)\pi/(4n+2) for the 1-based index k' = k + 1), the approximation

   x \approx \cos\psi, \qquad \psi = \phi + \frac{2\cos\phi}{(4n+2)^2 \sin\phi}\left(1 - \frac{11}{(4n+2)^2 \sin^2\phi}\right)

   of Petras [Pet1999] has an explicit error bound, of the order of \phi^2 \cdot 22 y^4 (1 + 2y^2) with y = 1/((4k'-1)\pi).

2. Newton iteration with a priori bounds. Petras also shows that Newton's method converges quadratically from this approximation, so the error bound is squared at each step; the iteration runs, at a precision about 20% above the target, until the bound is below 2^{-\mathit{prec}}. This gives a rigorous enclosure for modest precision.
3. Interval Newton refinement. At high precision (when 40 bits plus a padding of 8 + 2\log_2 n bits is less than half of prec), the steps above are carried out only at a low starting precision, and the enclosure is refined by interval Newton steps x \mapsto m - P_n(m)/P_n'(X) at doubling precisions, where m is the midpoint of the current enclosure X and P_n'(X) is bounded by P_n'(m) plus a bound for P_n'' from arb_hypgeom_legendre_p_ui_deriv_bound times the radius.
4. Weight. P_n'(x_k) is computed with arb_hypgeom_legendre_p_ui at the enclosure of x_k.

All evaluations of P_n and P_n' go through arb_hypgeom_legendre_p_ui, which chooses between the hypergeometric series at 0 and at 1, the asymptotic expansion in n and the recurrence (see Legendre functions and orthogonal polynomials).

Example#

The 20-point rule integrates e^x over [-1, 1] to more than 30 digits; the error of the rule itself is about 10^{-60}:

#include <flint/arb.h>
#include <flint/arb_hypgeom.h>

int main(void)
{
    arb_t x, w, s, t, wsum;
    ulong n = 20, k;
    slong prec = 128;

    arb_init(x);
    arb_init(w);
    arb_init(s);
    arb_init(t);
    arb_init(wsum);

    for (k = 0; k < n; k++)
    {
        arb_hypgeom_legendre_p_ui_root(x, w, n, k, prec);
        if (k == 0)
        {
            flint_printf("x_0 = "); arb_printn(x, 30, 0); flint_printf("\n");
            flint_printf("w_0 = "); arb_printn(w, 30, 0); flint_printf("\n");
        }
        /* s += w_k exp(x_k) */
        arb_exp(t, x, prec);
        arb_addmul(s, w, t, prec);
        arb_add(wsum, wsum, w, prec);
    }

    flint_printf("sum of weights = "); arb_printn(wsum, 30, 0); flint_printf("\n");
    flint_printf("quadrature     = "); arb_printn(s, 30, 0); flint_printf("\n");

    /* exact value e - 1/e */
    arb_const_e(t, prec);
    arb_inv(x, t, prec);
    arb_sub(t, t, x, prec);
    flint_printf("e - 1/e        = "); arb_printn(t, 30, 0); flint_printf("\n");

    arb_clear(x);
    arb_clear(w);
    arb_clear(s);
    arb_clear(t);
    arb_clear(wsum);
    flint_cleanup();
    return 0;
}

x_0 = [0.993128599185094924786122388471 +/- 3.21e-31]
w_0 = [0.0176140071391521183118619623519 +/- 4.72e-32]
sum of weights = [2.00000000000000000000000000000 +/- 1.08e-37]
quadrature     = [2.35040238728760291376476370119 +/- 1.21e-30]
e - 1/e        = [2.35040238728760291376476370119 +/- 1.21e-30]

The enclosure of the quadrature sum is rigorous for the sum, but it does not include the error of the quadrature rule. For an integral with a proved error bound, use acb_calc_integrate, which bounds the quadrature error using the analyticity of the integrand [Joh2018a].