Legendre functions and orthogonal polynomials

The Legendre polynomials $P_n(x)$, and more generally the classical orthogonal polynomials (Jacobi, Gegenbauer, Chebyshev, Laguerre, Hermite), are terminating hypergeometric series, for example

For integer degree they can be evaluated directly as polynomials, but for large $n$ this is slow and numerically unstable: the coefficients grow exponentially while $|P_n(x)| \le 1$ on $[-1,1]$. FLINT has a dedicated algorithm for $P_n(x)$ with real $x$ and word-size $n$, which also gives the nodes and weights of Gauss–Legendre quadrature; the other functions use the general ${}_2F_1$ and ${}_1F_1$ code with a few direct methods for integer degrees.

Legendre polynomials#

The algorithm of arb_hypgeom_legendre_p_ui is that of Johansson and Mezzarobba [JM2018]. It computes $P_n(x)$ and $P_n'(x)$ together for real $x$, and by $P_n(-x) = (-1)^n P_n(x)$ only $x \ge 0$ is needed. Four expansions are available.

• Expansion at $x = 1$ (arb_hypgeom_legendre_p_ui_one):

  $$P_n(x) = \sum_{k=0}^{n} \binom{n}{k}\binom{n+k}{k} \left(\frac{x-1}{2}\right)^k.$$

  The terms decrease quickly when $|1 - x|$ is small, and the series can be truncated after $K < n+1$ terms with a bound for the remainder.

• Expansion at $x = 0$ (arb_hypgeom_legendre_p_ui_zero): for even $n = 2m$,

  $$P_{2m}(x) = (-1)^m \frac{(1/2)_m}{m!}\, {}_2F_1\left(-m, m + \tfrac12; \tfrac12; x^2\right),$$

  and an analogous formula with $x\,{}_2F_1(-m, m + \frac32; \frac32; x^2)$ for odd $n$. The terms decrease when $|x|$ is small.

• Asymptotic expansion (arb_hypgeom_legendre_p_ui_asymp): with $x = \cos\theta$, $0 < \theta < \pi$, the Stieltjes series

  $$P_n(\cos\theta) = \frac{4}{\pi} \frac{\Gamma(n+1)}{\Gamma(n+\frac32)} \sum_{k=0}^{K-1} \frac{(\frac12)_k^2}{k!\,(n+\frac32)_k} \frac{\cos\left((n + k + \tfrac12)\theta - (k+\tfrac12)\tfrac{\pi}{2}\right)}{(2\sin\theta)^{k+1/2}} + R_K,$$

  summed as the real part of a complex series in $(1 - i\cot\theta)$, with the explicit bound $|R_K| \le 0.50795\, \frac{K!\, n!}{2^K (n+K)!\, y^K \sqrt{ny}}$, $y = \sin\theta$, which follows from the analysis of Bogaert, Michiels and Fostier [Bog2012]. The number of terms needed decreases as $n$ grows; for $n$ much larger than the precision a handful of terms suffice, independently of $n$.

• Recurrence (arb_hypgeom_legendre_p_ui_rec): the three-term recurrence $(k+1)P_{k+1} = (2k+1)xP_k - kP_{k-1}$, evaluated in fixed-point arithmetic with machine words and a rigorous bound for the accumulated rounding error. It is stable on $(-1, 1)$, costs $O(n)$ cheap operations, and is the fastest method for moderate $n$ and precision.

arb_hypgeom_legendre_p_ui selects the method from estimates of the cost of each expansion: the number of terms, the precision lost to cancellation in the series at 0 and at 1 (whose terms can be much larger than $P_n(x)$), and the size of the error term of the asymptotic expansion. The recurrence is used for $n$ below a few hundred at low precision, unless $x$ is very close to 0 or 1; for $|x| > 1$ only the expansion at 1 is used (with an early stop if it converges quickly). The working precision is raised by the estimated cancellation, and for inexact $x$ the function evaluates at the midpoint and adds the error propagated from the radius $r$, bounded by $r \cdot \max|P_n'|$ and $r \cdot \max|P_n''|$ on the ball (arb_hypgeom_legendre_p_ui_deriv_bound, which uses $|P_n'(x)| \le \min((n+1)^2/2,\ c\sqrt{n}\,(1 - x^2)^{-3/4})$ and $|P_n''(x)| \le \min((n+1)^4/8,\ 2c\, n^{3/2}(1-x^2)^{-5/4})$ for $|x| \le 1$, with $c \ge 2^{3/2}/\sqrt\pi$).

Gauss–Legendre nodes and weights#

The roots $1 > x_0 > \cdots > x_{n-1} > -1$ of $P_n$ and the weights $w_k = 2/((1-x_k^2)P_n'(x_k)^2)$ are computed by arb_hypgeom_legendre_p_ui_root: an initial approximation with an explicit error bound due to Petras [Pet1999], Newton iteration with a priori error bounds at low precision, and interval Newton iteration at doubling precisions above; each step evaluates $P_n$ and $P_n'$ with the method above, and the doubling of the precision means that the final evaluation dominates the cost. This is what makes high-degree Gauss–Legendre quadrature in acb_calc_integrate practical at high precision [JM2018]. See Roots of Legendre polynomials and Gauss–Legendre nodes.

Legendre functions#

For complex degree $n$, order $m$ and argument $z$, the associated Legendre function of the first kind is

with the regularized ${}_2F_1$ (type 0 has cuts on $(-\infty, -1]$ and $[1, \infty)$, type 1 on $(-\infty, 1]$), and the function of the second kind is

(and analogously for type 1). acb_hypgeom_legendre_p evaluates the ${}_2F_1$ representation with acb_hypgeom_2f1, which handles the transformations needed for large $|z|$. For integers $0 \le m \le n$ it uses instead the recurrence of acb_hypgeom_legendre_p_uiui_rec for $(1-z^2)^{-m/2}P_n^m(z)$, starting from $P_m^m$ and $P_{m+1}^m$, when the precision is low compared to $n$, $n - m$ is small, or $z$ is near 1 (where the transformations are singular); real $z$ with $m = 0$ goes to arb_hypgeom_legendre_p_ui. If $n - m > 2p$ the recurrence would lose all accuracy and the result is indeterminate.

For integer $m$, the defining formula for $Q_n^m$ is a $0 \cdot \infty$ limit, so acb_hypgeom_legendre_q uses connection formulas without the removable singularity. For type 0 it evaluates two ordinary ${}_2F_1$ series in $z^2$ with gamma and trigonometric prefactors [WQ3a]. For type 1 it evaluates one regularized ${}_2F_1$ in $1/z^2$ [WQ3b] when the ball $z$ is away from the real half-line $[-1, \infty)$ and from the square $|\operatorname{Re} z|, |\operatorname{Im} z| < 1/4$ (the formula is invalid on $(-1, 0)$), and otherwise combines the type 0 function $Q_n^m$ and $\mathcal{P}_n^m$ [WQ3c]. (The formula [WQ3d] would also be convenient, but it is incorrect in the lower half-plane.)

The spherical harmonics

are polynomials in $\cos\theta$ and $\sin\theta$, which acb_hypgeom_spherical_y evaluates with acb_hypgeom_legendre_p_uiui_rec.

Other orthogonal polynomials#

The other classical families are defined for complex degree through ${}_2F_1$ or ${}_1F_1$:

The direct methods cover the cases where the hypergeometric representation is singular although the polynomial is defined (the lower parameter $a + 1$, $m + \frac12$ or $m + 1$ is near a nonpositive integer). For non-integer degree, the Chebyshev functions use $T_\nu(\cos\theta) = \cos(\nu\theta)$ and $U_\nu(\cos\theta) = \sin((\nu+1)\theta)/\sin\theta$ (with $\cosh$ for $z$ outside $[-1,1]$) away from $z = 1$, and the ${}_2F_1$ series near it.

Use in FLINT#

• arb_hypgeom_legendre_p_ui and its methods arb_hypgeom_legendre_p_ui_zero, arb_hypgeom_legendre_p_ui_one, arb_hypgeom_legendre_p_ui_asymp, arb_hypgeom_legendre_p_ui_rec; arb_hypgeom_legendre_p_ui_root for quadrature, used by acb_calc_integrate and arb_fpwrap.
• acb_hypgeom_legendre_p, acb_hypgeom_legendre_q, acb_hypgeom_legendre_p_uiui_rec, acb_hypgeom_spherical_y and the real versions.
• acb_hypgeom_chebyshev_t, acb_hypgeom_chebyshev_u, acb_hypgeom_jacobi_p, acb_hypgeom_gegenbauer_c, acb_hypgeom_laguerre_l, acb_hypgeom_hermite_h and the real versions.

@misc{flint-doc-algorithms-special-functions-legendre,
  author       = {{The FLINT team}},
  title        = {{Legendre functions and orthogonal polynomials}},
  howpublished = {FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev},
  year         = {2026},
  url          = {https://flintlib.org/doc/algorithms/special-functions/legendre.html},
  note         = {Generated 2026-11-02}
}