Efficient computation of Laguerre polynomials

An efficient algorithm and a Fortran 90 module (LaguerrePol) for computing Laguerre polynomials $L^{(\alpha)}_n(z)$ are presented. The standard three-term recurrence relation satisfied by the polynomials and different types of asymptotic expansions valid for $n$ large and $\alpha$ small, are used depending on the parameter region. Based on tests of contiguous relations in the parameter $\alpha$ and the degree $n$ satisfied by the polynomials, we claim that a relative accuracy close or better than $10^{-12}$ can be obtained using the module LaguerrePol for computing the functions $L^{(\alpha)}_n(z)$ in the parameter range $z \ge 0$, $-1<\alpha \le 5$, $n \ge 0$.


Introduction
As is well known, Laguerre polynomials L (α) n (z) are involved in a vast number of applications in physics (quantum mechanics, plasma physics, etc) and engineering; for example see [1,2]. Also, the evaluation of Laguerre polynomials is central in the computation of nodes and weights in Gauss-Laguerre quadrature rules.
In this paper, we present an algorithm for computing Laguerre polynomials based on the use of the standard three-term recurrence relation satisfied by the polynomials and three types of asymptotic expansions valid for n large and small values of the parameter α: two Bessel-type expansions (used in the oscillatory region of the functions) and a uniform Airy-type expansion. This Airy expansion is specially suitable in the transition between the oscillatory and monotonic regimes of the Laguerre polynomials, although its domain of applicability extends to a large part of both the oscillatory and the monotonic regions.
The resulting algorithm, implemented in the Fortran 90 module Laguer-rePol is accurate and particularly efficient for large values of the parameter n.
An example of the behaviour of the Laguerre polynomials in the oscillatory region is shown in Figure 1.
Next, we are going to describe the theoretical expressions involved in the computation of Laguerre polynomials.

Asymptotic expansions for n large
In the asymptotic expansions for large degree we assume the α is fixed, which means small with respect to n. If we need results for which α is not small enough, we can use recursion with respect to α. That is, 150 (z) is plotted as an example. The value of the parameter ν (defined in (2)) is 610 in this particular case.
This follows from the corresponding c-recursion of the Kummer function U(−n, c, x) = (−1) n n! L (c−1) n (x).

An expansion in terms of Airy functions
We start with the representation 2 We have the relation The first coefficients of the expansions in (6) are where b = √ ζ if ζ ≥ 0 and b = i √ −ζ when ζ ≤ 0, and

A simple Bessel-type expansion
For the Laguerre polynomials we consider two types of asymptotic expansions in terms of Bessel functions, one for small values of the variable z of L (α) n (z) and one in which larger values are allowed. Next we give some details of the asymptotic expansion valid for small values of the variable z.
We use the expansion of the Kummer function 1 F 1 (a; c; x) for large negative values of a. First we mention Then, see [6, §10.3.4], This expansion of 1 F 1 (−a; c; z) is valid for bounded values of z and c, with a → ∞ inside the sector −π + δ ≤ ph a ≤ π − δ. This gives for the Laguerre polynomial (14) The coefficients a k (x) and b k (x) follow from the expansion of the function The function f is analytic in the strip |ℑs| < 2π and it can be expanded for The coefficients c k (x) are combinations of Bernoulli numbers and Bernoulli polynomials, the first ones being (with c = α + 1) The coefficients a k (x) and b k (x) are in terms of the c k (x) given by k = 0, 1, 2, . . ., and the first relations are again with c = α + 1.

A not so simple expansion in terms of Bessel functions
In this case we use the representation 3 with expansions Here, with ζ given by We have the relation The first coefficients are We give a few details about the coefficients A j (ζ) and B j (ζ) of the expansions in (21). The first ones are given in (25).
First we need coefficients c ± k of the expansions where, for 0 ≤ x < 1, b, s ± and the relation between s and u are defined by with ζ defined in (23), and s ± being the saddle points of the s-function and ±ib of the u-function. Because s(u) is an odd function of u, we have In the following we write c + k = c k . The first coefficients are .
(29) Next we consider the function with expansion where χ(ζ) is defined in (22). Because h(u) is even, an expansion at the other saddle point −ib has coefficients (−1) k d k ; also, d 0 = h(ib) = 1. This makes A 0 (ζ) = 1, see (25). The coefficient d 1 is given by where γ = α + 1. The coefficients in (21) follow from the following recursive scheme k = 0, 1, 2, . . ., with starting value h 0 (u) = h(u). The coefficients α k and β k follow from substituting u = ±ib. Because h(u) is even, β 0 = 0, and h 1 (u) is odd, giving α 1 = 0, and so on. Then, the coefficients in the expansions in (21) are given by This gives, again with γ = α + 1, For small values of x we need expansions. We can expand in terms of ζ or x. For example, we can write The first coefficients are A

Expansions for large values of n and α
In [7] we have given expansions for large n in which α = O(n) is allowed; for a summary see [8]. These results can be obtained by using an integral representation, but they follow also from uniform expansions of Whittaker functions obtained by using differential equations; see [9]. These expansions include the J-Bessel function, and are valid in the parameter domain where order and argument of the Bessel function are equal, that is, in the turning point domain. Because no explicit forms of the coefficients in the expansions are available, we omit further details.

An algorithm for computing the Bessel functions J ν (z)
The algorithm for computing the Bessel function J ν (z) in the expansions (14) and (20) is based in the following methods of approximation: Power series.
The power series given in Eq.(10.2.2) of [10, §10.19(ii)] is used for computing J ν (z) when z is small: .

Debye's asymptotic expansions.
Debye's asymptotic expansions are also used in the algorithm. When ν < z, we use The coefficients U k (p) are polynomials in p of degree 3k given by U 0 (p) = 1 and Asymptotic expansions for large z.
For large values of the argument z, we use the Hankel's expansion given in [10, §10.17(i)]: The coefficients a k (ν) are given by Airy-type expansions.
An important ingredient in our algorithm for computing Bessel functions are Airy-type expansions. We use the representation given in [11,Chapter 8] The variable ζ is written in terms of the variable x as Three-term recurrence relation using Miller's algorithm.
The standard three-term recurrence relation for the cylinder functions is computed backwards (starting from large values of ν) using Miller's algorithm.

Three-term recurrence relation
The generalized Laguerre polynomials satisfy the following three-term recurrence relation This recurrence relation is not ill conditioned in both backward and forward directions. Therefore it can be used with starting values L (α) 0 (z) = 1 and L (α) 1 (z) = 1 + α − z, to compute the generalized Laguerre polynomials when n is small/moderate. As n increases, it is more efficient, as we later discuss, to use the asymptotic expansions described in the previous sections.

Overview of the software structure
The Fortran 90 package includes the main module LaguerrePol, which includes as public routine the function laguerre.
In the module LaguerrePol, the auxiliary modules Someconstants (a module for the computation of the main constants used in the different routines), BesselJY (for the computation of Bessel functions) and AiryFunction (for the computation of Airy functions) are used.

Description of the individual software components
The calling sequence of this routine is laguerre(a,n,z,lagp,ierr) where the input data are: a, n and z (arguments of the Laguerre polynomial). The outputs of the function are error flag ierr and the value of the Laguerre polynomial value lagp. The possible values of the error flag are: ierr = 0, successful computation; ierr = 1, computation failed due to overflow/underflow; ierr = 2, arguments out of range.

Testing the algorithms
The performance of the asymptotic expansions for the Laguerre polynomials has been tested by considering the relation given in Eq.(18.9.13) of [12] written in the form This chek fails close to the zeros of L (α+1) n ; in this case, we can consider the alternative test Notice that, because the zeros of L (α) n and L (α+1) n are interlaced, both tests will not fail simultaneously. We can therefore take A test for the accuracy obtained using the Bessel-type and Airy-type expansions for n large is shown in Figures 2, 3 and 4, respectively. In these figures, the points where the value of ǫ in Eq. (40) is greater than 5 10 −12 are shown. Parameter values have been randomly generated in the oscillatory region of the Laguerre polynomials with (α, n) ∈ (−1, 10) × (200, 10000).  As can be seen, the Bessel-type expansions (Figures 2, 3) are valid for small values of the variable z/ν; in particular, the simple Bessel-type expansion works well when the argument of the Laguerre polynomials is very close to the origin. On the other hand, the domain of applicability of the Airy-type expansion extends well beyond the transition between the oscillatory and the monotonic regions of the functions, as Figure 4 shows for the oscillatory region. The validity of the three asymptotic expansions is, in all cases, limited by the value of the parameter α, as commented in §(2.1.4). In the Fortran 90 module, we restrict the values of this parameter to the interval (−1, 5) in order to avoid the use of the recursion relation (37) for large n.
We have also tested the efficiency of using the asymptotic expansions in their region of applicability in comparison with the use of the three-term recurrence relation given in Eq. (37) for computing the functions. Our tests show that for n > 200 it is more efficient to use the asymptotic expansions than the recurrence relation. As an example, Table 1 shows few of the test values for particular choices of the parameters.

Test run description
The Fortran 90 test program testlag.f90 includes the computation of 25 function values and their comparison with the corresponding pre-computed results. Also, the relation given in (38) is tested for several values of the parameters (z, α, n).