Asymptotic approximations to the nodes and weights of Gauss-Hermite and Gauss-Laguerre quadratures

Asymptotic approximations to the zeros of Hermite and Laguerre polynomials are given, together with methods for obtaining the coefficients in the expansions. These approximations can be used as a standalone method of computation of Gaussian quadratures for high enough degrees, with Gaussian weights computed from asymptotic approximations for the orthogonal polynomials. We provide numerical evidence showing that for degrees greater than $100$ the asymptotic methods are enough for a double precision accuracy computation ($15$-$16$ digits) of the nodes and weights of the Gauss--Hermite and Gauss--Laguerre quadratures.

Iterative algorithms are interesting methods of computation of Gaussian nodes and weights, very clearly outperforming matrix methods (Golub-Welsch [1]) for high degrees. They are based on the computation of the roots of the orthogonal polynomial by an iterative method and the subsequent computation of the weights by using function relations such as those in Eqs. (2)- (4). Most iterative methods for the computation of the Gaussian nodes (with the exception of [2]) require accurate enough first approximations to ensure the convergence of the iterative method (typically the Newton method); for two recent examples, see [3,4]. An alternative approach [5], although less efficient for high degrees than iterative methods with asymptotic first approximations [3,4], consists in guessing these first approximations by integrating a Prufer-transformed ordinary differential equation (ODE) with a Runge-Kutta method, and then refining these guesses by the Newton method (however, asymptotic approximations were also used in this reference for the particular case of Gauss-Legendre quadrature). More recently, noniterative methods based on asymptotic approximations for the computation of Gauss-Legendre nodes and weights were developed in [6], which were shown to outperform iterative approaches.
In this paper, our aim is to provide asymptotic approximations for the accurate computation of the nodes and weights of Gauss-Hermite and Gauss-Laguerre quadrature. These approximations provide a fast and accurate method of computation, which can be used for arbitrarily large degree, but which also provide accurate results for not so large degrees (n ≥ 100). The methods are able to compute both the nodes and the weights with nearly double precision accuracy, improving the accuracy of the available fixed precision iterative methods.
As we will discuss in a subsequent paper, a fully noniterative approach is also possible for the case of Gauss-Jacobi quadrature [7], similarly as was shown for the particular case of Legendre polynomials [6].

Hermite polynomials
In [4], first estimates of the zeros of Hermite polynomials are based on work of Tricomi for the middle zeros; these first guesses follow from expansions in terms of elementary functions. For the remaining zeros near the positive endpoint √ 2n + 1 of the zeros interval, the first estimates are taken from the work of Gatteschi, and are in terms of the zeros of the Airy functions.
In this section, we give an expansion of the zeros based on the asymptotic expansion in terms of Airy functions. The expansion can be used for all positive zeros, however, the approximations are less accurate for the small zeros. For these we give an approximation based on an asymptotic expansion in terms of elementary functions. We start discussing this expansion.

Expansions in terms of elementary functions
An expansion in terms of elementary functions for the Hermite polynomials is given in [8, §18.15(v)] with a limited number of coefficients. However, we prefer an expansion for the parabolic cylinder function derived in [9]; these results are summarized in [10, §12.10(iv)] and [11, §30.2.3].
We use the notations and we have the asymptotic representation with expansions where δ is an arbitrary small positive number.
The first few coefficients are and more u s (t) follow from the recurrence relations The quantity g(μ) is only known in the form of an asymptotic expansion where the coefficients γ k are defined by The first ones are .
2.1.1. Expansions of the zeros. Next we discuss expansions for the zeros of H n (x), x k , 1 ≤ k ≤ n (x 1 < x 2 < · · · < x n ). We introduce a function W (η) (see (7)) and try to solve the equation W (η) = 0 for large values of n. We define a first approximation η 0 such that the cosine term vanishes and η 0 and the corresponding t and x-values are (in first-order approximation) related to a zero of H n (x). The small zeros are around x = 0 and t = 0, that is, for η near η(0) = 1 4 π . We define In this way, cos(μ 2 η 0 − 1 4 π ) = 0, and this choice of η 0 follows from the location of the zeros of the cosine function and those of H n (x). Observe that, when n is odd and k = 1 2 (n + 1), that is, x k = 0, it follows that η 0 = 1 4 π . If η = 1 4 π we have t = 0 and x = 0. We assume that the equation W (η) = 0 has a solution η that can be expanded in the form and consider the Taylor expansion and equation where W (η) and its derivatives are taken at η = η 0 . Because the expansions in (8) are in terms of t, we need dt/dη = −1/ √ 1 − t 2 . When we have found η, the corresponding t-value is obtained by inverting the relation for η(t) in (6). For this purpose we use the expansion It is also possible to invert the relation (19) by using an iterative method. For this purpose it is convenient to write t = sin 1 2 θ . Then the equation to be solved for θ ∈ (−π, π) reads A Newton or related procedure can be used to solve this equation, but in our algorithms we prefer to use the series shown in (19), which is faster (and of more restricted applicability, but sufficient for our purposes).
After a few symbolic manipulations, we find that η 2k+1 = 0, k = 0, 1, 2, . . ., and that the first nonzero coefficients are Because we have a recurrence relation for the coefficients u s (t) in (10), it is quite easy to generate many u s (t) and also much more coefficients η j than given in (21).
Algorithm: For the computation of the approximations of the zeros x k , we summarize the procedure as follows.

Expansions in terms of Airy functions
For the large zeros we use the Airy-type expansion of the Hermite polynomials. We write (see [10, section 12.10(vii) with expansions where μ = √ 2n + 1, t = x/μ, and g(μ) is the function with asymptotic expansion given in (11). For ζ we have the definition where η(t) is defined in (6); χ (ζ ) is defined by The variable ζ is analytic in a neighborhood of t = 1. We have the differential equation and we have the following expansions in powers of t − 1 and ζ the expansions The relation between t and ζ is singular at t = −1, ζ (−1) = −(3π/4) 2/3 = −1.770 · · · , and the series in the second line converges for |ζ | < 1.770 · · · . The coefficients are given by 2s m=0 β m (χ (ζ )) 6(2s−m) u 2s−m (t), where u s (t) are as in Section 2.1, and A recursion for α m reads The first few coefficients of the expansions in (23) are given by: (31) Here, χ = χ (ζ ) is given by (25). To avoid numerical cancellations when ζ is small in the above representations, we can expand the coefficients, which are analytic at ζ = 0, in powers of ζ .

Expansions of the zeros.
An expansion for the zeros is obtained as follows. First, we determine the zeros in terms of ζ .
For the first-order approximation of a zero x n−k+1 of H n (x) we compute where a k is a zero of the Airy function Ai(x). Because of the symmetry of the Hermite polynomial, we assume that 1 ≤ k ≤ 1 2 n . We introduce an expansion of ζ corresponding to the zero of H n (x) by writing and we try to obtain the ζ j , j ≥ 1. We introduce a function W (ζ ) by writing (see (22)) and expand W (ζ ) at ζ = ζ 0 , writing ζ = ζ 0 + ε, which gives In this equation we substitute the expansion given in (32) and those in (23), compare equal powers of μ and obtain the first few coefficients where the derivative is with respect to ζ and the coefficients are given in (28).
To obtain the derivative of B 0 (ζ ) we need which follow from (25) and (26). This gives For small values of ζ we have expansions of the form Algorithm: When we have obtained a value ζ that corresponds to a zero of the Hermite polynomial, the corresponding t-value should be obtained from the second equation in (24). This equation has to be solved by a numerical procedure. A first estimate, when ζ is small, can be obtained from the second line in (27), and more terms of that expansion can easily be obtained by a symbolic package.
For an iterative procedure it is convenient to substitute t = cos 1 2 θ , with θ ∈ [0, 2π ). Then the equation to be solved for θ reads 8 3 (−ζ ) 3 2 = θ − sin θ and we can use, for instance, the Newton method for this purpose. However, in our algorithms we prefer to invert using enough terms in (27), which is a faster method.
We proceed as follows for computing approximations for the zeros.
(1) To approximate the zero x n−k+1 , define the starting value ζ 0 = μ − 4 3 a k , 1 ≤ k ≤ 1 2 n, where a k is a zero of the Airy function. (2) Compute t from the second line of (27).

Numerical performance of the expansions
The approximation (17) (obtained from the expansion in terms of elementary functions) is accurate for large n and particularly for the small zeros. As a first numerical example of the accuracy, even for quite small n, we take n = 11, k = 7 (the smallest positive zero). Then, η 0 = 0.
where δ is an arbitrary small positive number. Hence, for the large zeros this method is not reliable, and we need to restrict the number of zeros that we can compute. For example, we can request that |t| ≤ 1 2 , the corresponding η-value satisfies |η − 1 4 π | ≤ 1 12 π + 1 8 √ 3 = 0.478. When we use the first estimate η 0 given in (16) in the equation |η 0 − 1 4 π | ≤ 0.478, we find for k the bound | 1 2 n − k| 0.478 π (2n + 1) = 0.304n + 0.152. This says that roughly 0.3n of the positive zeros can be computed by using the asymptotic approximations of Section 2.1, when we request |t| ≤ 1 2 . In practice, as we will see later, the expansions in terms of elementary functions can be used for larger values of |t| and when they are accurate, they are preferable to the expansions in terms of Airy functions because the algorithm is faster.
More extensive tests of the expansions have been performed using finite precision implementations coded in Fortran 90. In these implementations only noniterative methods (power series) are used for the inversion of the variables. Figure 1 shows the performance of the expansion in terms of elementary functions. In this figure, the relative accuracy obtained for computing the positive zeros of H n (x) for n = 100, 1, 000, 10, 000 is plotted. The label i in the abscissa represents the order of the zero (starting from i = 1 for the smallest positive zero). The algorithm for testing the accuracy of the zeros has been implemented in finite precision arithmetic using the first six nonzero terms in the expansion. We compare the asymptotic expansions against an extended precision accuracy (close to 32 digits) iterative algorithm, which uses the global fixed point method of [2], with orthogonal polynomials computed by local Taylor series. As can be seen, a very large number of the zeros for the three values of n tested can be computed with the expansion with a relative accuracy near full double precision. Actually, the points not shown in the plot correspond to values with all digits correct in double precision accuracy. However, the expansion fails for the largest zeros, as expected.
As for the asymptotic expansion in terms of the zeros of Airy functions (32), the situation is the reverse: the further we are from the turning point at t = 1 (ζ = 0), the larger the relative errors become. Therefore, for n fixed the maximum errors in the computation are obtained for the small zeros. For example, using Maple with digits = 16, we take n = 11 and six coefficients in (32). Then we have for the zero x 6 at the origin ζ 0 = μ − 4 3 a 6 = −1.115618210110694, t = −0.1668495251592333 × 10 −3 , and the better values ζ = −1.115460237225190 and t = 1.746192313216916 × 10 −13 . This gives x 6 . = 8.374444141492045 × 10 −13 and for the largest zero x 11 the relative accuracy is 10 −15 . A test of the expansion for very large values of n using a finite precision arithmetic implementation is shown in Fig. 2. In this figure, we show the relative accuracy obtained with the asymptotic expansion (32) for computing the largest 1, 000 positive zeros of H n (x) for n = 10, 000, 100, 000, 1, 000, 000. As can be seen, an accuracy near 10 −16 can be obtained in all cases. The zeros a k of the Airy function have been computed using a k = −T ( 3 8 π (4k − 1)), where T (t) has the Poincaré's expansion (see [12, §9.9(iv)])  This expansion is valid for moderate/large values of k. In our implementation we use precomputed values for the first 10 zeros of the Airy function and the Poincaré's expansion for the rest.
The accuracy of the two expansions (17) and (32) for approximating the zeros of Hermite polynomials for n = 100 is compared in Fig. 3. As can be seen, the combined use of both expansions allow the computation of all the zeros with a double precision accuracy of 15-16 digits.
In Table 1, we illustrate the efficiency of the expansions for approximating the zeros of Hermite polynomials for n = 100, 10, 000. In particular, the first 0.6n zeros of the Hermite polynomials have been computed with the asymptotic expansion in terms of elementary functions and the last 0.4n zeros with the asymptotic expansion in terms of the zeros of Airy functions. With this splitting and by taking enough terms, it is possible to use the series (19) and (27) for computing the t-values in the expansions instead of using an iterative method for solving the nonlinear equations. In the table we show average CPU times (obtained using an Intel Core i54310U 2.6 GHz processor under Windows) per node. The second column shows the CPU times when the number of terms required (no more than five or six depending on the expansion) for a double precision accuracy for the zeros is considered, while the first column shows the CPU times for only two terms. For n = 10, 000, this is the number of terms needed in the expansions to obtain double precision accuracy. For n = 100 we observe that there is not much difference in speed between the more simple (two terms) and the more accurate approximation; this favors the use of accurate asymptotic approximations with no ulterior iterative refinements. The table also shows that the computation of the expansion in terms of elementary functions is more efficient than the expansion in terms of zeros of Airy functions although for n = 10, 000 the difference in speed is not very significant.
Once the nodes (the zeros of H n (x)) of the Gauss-Hermite quadrature have been computed, approximations to the weights given in (2) can be also obtained by using the asymptotic results in Section 2.1 (elementary functions) and Section 2.2 (Airy functions).
For the computation of the weights, one needs to be careful to avoid overflows in the computation both as a function of n and as a function of the values of the nodes. With respect to the dependence on n, we observe that the large factor 2 n n! in (2) can be cancelled out by the factors in front of the expansions (7) and (22). This is as expected because using the first approximations from the elementary asymptotic expansions as n → ∞ we obtain the estimate for the weights: This estimation shows that underflow may occur for computing the large zeros. In this case the range of computation of the weights can be enlarged by scaling the factor e x 2 /2 in the asymptotic approximations and computing scaled weights given byw With this, the overflow/underflow limitations are eliminated. Using (2) this scaled weight can be written as This expression does not have overflow/underflow limitations neither with respect to x nor with respect to n. Using (7) or (22) we observe that the dominant factors e −x 2 /2 and 2 n+1 n! can be explicitly cancelled out.
Another interesting property of this expression is that it is well conditioned with respect to the values of the nodes. Indeed, we havew i = W (x i ), where we define the function W (x) = √ π2 n+1 n!/y (x) 2 . Now, it is straightforward to check that W (x i ) = 0, which means that, at the nodes x = x i , the value of the weight is little affected by variations on the actual value of the node. This, as we will show, will allow us to compute scaled weights with nearly full double precision in all the range.
For computing the scaled weights in this way, we need to compute y (x) from the asymptotic expansions (7) or (22). This is a straightforward computation and, for instance, starting from (7) we have that and The dots mean differentiation with respect to t. Two examples of computation of the scaled weights (for n = 1, 000, 10, 000) using the expansion in terms of elementary functions are shown in Fig. 4. As can be seen, most of the scaled weights can be computed with almost double precision accuracy. Also, as expected, there is some loss of accuracy for the weights corresponding to the largest nodes (as discussed, for these values one has to use the expansion for the Hermite polynomials in terms of Airy functions).
Typically, the additional computation of the weights requires about 70% more CPU time than when using the asymptotic expansion in terms of elementary functions and about 133% more CPU time than when using the asymptotic expansion in terms of Airy functions (due to the computation of these functions). This shows that, when possible, the direct computation of nodes and weights using asymptotics will be more efficient than computing more crude first approximations and then refining with an iterative method which uses values of the orthogonal polynomial. Each time the function (and its derivative when we use Newton's method) is computed, the CPU time increases by this same amount, and only when one iteration is needed the speed would be comparable.

Laguerre polynomials
We consider asymptotic expansions for the Laguerre polynomials L (α) n (x) in terms of Bessel functions, Airy functions and Hermite polynomials. Some of these expansions have been used to build an efficient scheme for computing the Laguerre polynomials for large values of n and small values of α (−1 < α ≤ 5) [13]. We discuss how to use the expansions to obtain approximations to the zeros of Laguerre polynomials. Later, in Section 3.5 we give expansions valid for large n and α.
For a survey of the work of several authors on inequalities and asymptotic formulas for the zeros of L (α) n (x) as n or α or ν = 4n + 2α + 2 → ∞, we refer to [14]. See also [15], were an alternative method, based on nonlinear steepest descent analysis of Riemann-Hilbert problems, is given for Laguerre-type Gaussian quadrature (and in particular Gauss-Laguerre).

A simple Bessel-type expansion
We have the following representation 1 with expansions valid for bounded values of x and α.
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) (50) 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.

Expansions of the zeros. Approximations of the zeros of L (α)
n (x) can be obtained from (46) and expressed in terms of zeros of the Bessel function J α (x). Because the expansion is valid for bounded values of x, the approximation can only be used for the small zeros. For example, in Table 2  Table 2 Relative Errors in the Computation of the Zeros x k (See (54) we show the results for the first 10 zeros when n = 100, and for these early zeros the approximations are satisfactory. We write (see (46)) A first approximation to the zero x k of L (α) n (x) follows from writing 2 √ nx k = j k , where j k is the kth zero of J α (x). A further approximation will be obtained by writing By expanding W (x) at the zero x = ξ + ε, assuming that ε is small, we find and substituting an expansion of the form we find the following first few values ξ 1 = ξ 12 (ξ − 6(α + 1)), where ξ is defined in (54).
Algorithm and first numerical examples for the zeros: The algorithm for computing the asymptotic approximation of the zeros runs in the same way as described for the Hermite polynomials, but is quite simple now. First compute ξ from (54) and the ξ j given in (57), then compute ε from (56), and finally x k from (54).
In Table 2, we show the results of a first numerical verification for the expansion. We take n = 100, α = 1 3 , and compute the first 10 zeros by using Maple with digits = 32. We show the relative errors in our approximations when we take 2, 4, and 6 terms in the expansion (56). As can be seen in the table, it is possible to obtain an accuracy near double precision (10 −16 ) in the computation of the first two zeros of L (1/3) 100 (x) using just the expansion with six terms.

An expansion in terms of Airy functions
We start with the representation 2 L (α) n (νσ ) = (−1) n e 1 2 νσ χ (ζ ) 2 α ν with expansions uniformly for bounded α and σ ∈ (σ 0 , ∞], where σ 0 ∈ (0, 1), a fixed number. Here, and We have the relation For the derivative, we can use the relation The first coefficients of the expansions in (59) are , (65) More coefficients can be obtained by the method described in [11, §23.2]. Starting point in this case is the integral (see [17,§VII.5,(5.11) where L is an Airy-type contour and f (u) is given by The relation between z and u follows in this case from the cubical transformation The function f (u) can be expanded in a two-point Taylor series in which the coefficients can be expressed in terms of the derivatives of f (u) at u = ± √ ζ . An integration by parts procedure then gives the coefficients α 2 j and β 2 j+1 of (59).
In Section 3.3.1 we describe in detail this method for a Bessel-type expansion.

Expansions of the zeros. We write
where A(ζ ) and B(ζ ) have the expansions shown in (59). Similarly as in Section 2 we write the zeros x j of L (α) n (x) in terms of the zeros a k of the Airy function. These zeros are negative, and a 1 will correspond the nth zero of L (α) n (x), a 2 with the (n − 1)th zero, and so on. A zero of L (α) n (x) is a zero of W (ζ ) and it can be written in terms of ζ in the form and we assume that we can expand By expanding W (ζ ) at ζ 0 we have and substituting the expansions shown in (59) we can obtain the coefficients ζ j . We obtain where β 1 is given in (64). The coefficients are evaluated at ζ 0 .
Algorithm and first numerical examples for the zeros: In Section 2.2.1 we have described the algorithm for computing the asymptotic approximation of the zeros for the Airy case. The present algorithm runs in the same way. For the zero x n+1− j , j = 1, 2, . . ., first compute ζ 0 , from (71). Then compute σ 0 by inverting the first relation in (61). This is done by using the expansion An alternative would be to use an iterative method. In that case it is convenient to write σ = cos 2 θ , and the equation to be solved for θ becomes With σ = σ 0 we compute the coefficients in (74), then ε and ζ from (72) and (71). A final inversion of the relation in the first line of (61) gives the σ , and then x n+1− j ∼ νσ .

Another expansion in terms of Bessel functions
After substituting t = e −s in the integral representation 3 we obtain the representation where ν = 2n + α + 1 and h(s, ρ) = s − ρ coth s. The contour starts at −∞ with ph u = −π , encircles the origin anticlockwise, and returns to −∞ with ph u = π . The transformation to a standard form for this case is where By using an integration by parts procedure (see Section 3.3.1), we can obtain the representation with expansions uniformly for ρ ≤ 1 − δ, where δ ∈ (0, 1) is a fixed number. Here, with ζ given by We have the relation The first coefficients are A 0 (ζ ) = 1, where More coefficients can be obtained by using the method described in Section 3.3.1.
To remove in (81) the singularities due to the Bessel functions at ζ = 0, it is convenient to use the function E ν (z) introduced by Tricomi; see [18, p. 34]. We have It is an analytic function of z. In terms of the modified Bessel function, we can write The representation in (81) can be written in the form and we can use this representation also for ζ < 0, i.e., ρ < 0.
For more details about the coefficients A j (ζ ) and B j (ζ ) of the expansions in (82), see [13].

A general method for the coefficients in Bessel-type expansions.
We describe a general method for evaluating the coefficients A k (ζ ) and B k (ζ ) used in (82).
We consider the standard form where the contour C starts at −∞ with ph u = −π , encircles the origin anti-clockwise, and returns to −∞ with ph u = π . The f (u) is assumed to be analytic in a neighborhood of C, and in particular in a domain that contains the saddle points ±ib, where b = √ ζ . When we replace f by unity, we obtain the Bessel function: The coefficients of the expansions in (82) follow from the recursive scheme with f 0 (u) = f (u), the coefficient function.
Using this scheme and integration by parts, we can obtain the asymptotic expansion The coefficients A j (ζ ) and B j (ζ ) can all be expressed in terms of the derivatives f (k) (±ib) of f (u) at the saddle points ±ib; we will need these for 0 ≤ k ≤ 2 j (see (98)).
We expand the functions f j (u) in two-point Taylor expansions Using (79), we derive the following recursive scheme for the coefficients for j, k = 0, 1, 2, . . ., and the coefficients A j and B j follow from In the present case of the Laguerre polynomials, the functions f 2 j are even and f 2 j+1 are odd, and we have A 2 j+1 (ζ ) = 0 and B 2 j (ζ ) = 0. A few nonvanishing coefficients are To have A 0 (ζ ) = 1 in the first expansion in (82) we have scaled all A and B coefficients with respect to A 0 (ζ ) = χ (ζ ); see (83).
Remark 1. The main step for obtaining the coefficients A j (ζ ) and B j (ζ ) is the evaluation of those for j = 0 in (95) and we summarize the method described in [19]. We rewrite the two-point Taylor expansion in the form where, in the present case, u 1 = −b and u 2 = b. Then, We have a 0 (u 1 , u 2 ) = f (b)/(2b) and a 0 (u 2 , u 1 ) = − f (−b)/(2b), and, for k = 1, 2, 3, . . ., a k (u 2 , u 1 ) follows from a k (u 1 , u 2 ) by replacing b by −b.

Expansions of the zeros.
From the Bessel-type expansion, we derive expansions of the first half of the zeros of the Laguerre polynomial. We write A zero of L (α) n (2νx) is a zero of W (ζ ) and it can be written in terms of ζ in the form where j k is a zero of J α (z). By expanding W (ζ ) we have with the zero ζ in this form We assume that ε can be expanded in the form and substituting this expansion, we obtain ζ 1 = −B 1 (ζ ) (see (86)) and (106) In the algorithm we use ζ = ζ 0 .
Algorithm and first numerical examples for the zeros: As in the previous cases we describe how the asymptotic approximations for the zeros can be obtained. For the zero x k , k = 1, 2, . . ., first compute ζ 0 , from (103). Then compute ρ 0 by inverting the second relation in (84). This is done by using the expansion An alternative is solving with an iterative method. In that case it is convenient to write ρ = sin 2 1 2 θ , and the equation to be solved for θ becomes 8 √ ζ = θ + sin θ , 0 ≤ θ < π. With ρ = ρ 0 we compute the coefficients ζ j in (105), see also (86). Compute ζ from (103) and perform a final inversion of the relation in the second line of (84). This gives the ρ, and then x k ∼ 2νρ.
Because the expansions in (82) become useless when ρ → 1, we should use the present result for a limited number of zeros, say, only for k = 1, 2, 3, . . . , 1 2 n; the remaining zeros can be obtained by using the Airy-type expansion.
In the next section, we analyze in more detail the performance of the different expansions for the zeros and we also discuss the stable computation of the weights. An implementation of the expansions in finite precision arithmetic (coded in Fortran 90) has been considered for testing. As for the Hermite case, in these implementations only noniterative methods (power series) are used for the inversion of the variables. For computing the first zeros of Bessel functions, we use the algorithm describe in [20]. For large zeros we use the MacMahon's expansion (see [21, §10.21(vi) where μ = 4ν 2 , a = (m + ν/2 − 1/4)π .
As can be seen in Fig. 5, the validity of the first asymptotic expansion in terms of zeros of Bessel functions (56) is limited to the first zeros. On the contrary, Fig. 7 shows that the other Bessel expansion (103) works very  well for approximating a large number of zeros of the Laguerre polynomial but fails for the last zeros. For these zeros, the Airy expansion (75) should be used. The accuracy of the Bessel and Airy expansions for n = 100 is illustrated in Fig. 8. As in the case of the Hermite approximations, the combined use of the expansions allow the computation of the zeros of Laguerre polynomials for n = 100 with an accuracy of 15-16 digits. The efficiency of the expansions is compared in Table 3. As in the case of the zeros of Hermite polynomials, to improve the speed of the methods we apply the expansions only in the regions where the inversion of the variables can be done accurately by using the series expansions (75) and (107) in the case of the Airy expansion and the second Bessel expansion (103), respectively: the first 0.75n zeros for the Bessel expansion and the last 0.25n zeros for the Airy expansion. For these two expansions, we observe in the table that there is no much difference in speed between using two terms and the more accurate approximation (for clarity the table  includes the number of terms needed for the three different expansions). With respect to the comparison between the different expansions, we observe that the computation of the first expansion in terms of Bessel functions is, as expected, extremely efficient in its range of validity. On the other hand, the expansion (103) in terms of Bessel functions is slightly more efficient than the Airy expansion.
As in the case of the Gauss-Hermite quadrature, overflow/underflow limitations in the computation of the weights can be eliminated by balancing the large terms as a function of n in the expressions and by scaling out the dependence on the weights. A first estimation of the weights as n → ∞ is given by The range of computation of the weights of Gauss-Laguerre quadrature can be enlarged by simply scaling out the dominant factor in the asymptotic expansions for Laguerre polynomials. When α is small, this factor is given by e x/2 . With this, one can define the scaled weights bỹ These normalized weights do not overflow/underflow as a function of n, α and x i . In addition, similarly as we did for the Hermite case, we can compute this scaled weights in a numerically stable way. We notice that the weights (3) can be written as where z = √ x, and therefore z i = √ x i . Now, in the new variable z, the scaled weights can be expressed as where the dots mean differentiation with respect to z and Now, we define W (z) = 4 (n + α + 1)/(n!(ẏ(z)) 2 ) and with this we have that w i = W (z i ), and it is straightforward to check that we have again the desirable property d dz W (z i ) = 0. This means that the computation is well conditioned in the sense that the error for the weights will be approximately proportional to the square of the error for the nodes. As a consequence, as we will shown, the weights can can be computed with almost no accuracy loss.
All that is left for computing the nodes is to use the expansions for the Laguerre polynomials to computeẏ(z) by differentiation. In particular, starting from (81) we havė where in this expression x is the variable defined in Section (3.3) and and in these equations prime denotes the derivative with respect to ρ. Similarly as we did for the Hermite case, we show in Fig. 9 two examples of computation of the scaled weights (112) for n = 1, 000, 10, 000 (with α = 1/4). We use the expansion in terms of Bessel functions (81). As can be seen, the accuracy for the scaled weights is better than 10 −15 in most cases. There is some loss of accuracy for the weights corresponding to the largest nodes (as discussed, for these values one has to use the expansion for the Laguerre polynomials in terms of Airy functions).

An expansion for large values of n and α.
In [24] we have given expansions for large n in which α = O(n) is allowed; for a summary see [25]. The results follow also from uniform expansions of Whittaker functions obtained by using differential equations; see [26]. 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. In this section, explicit expressions for the first few coefficients of the expansion are given.
By using an integral we can derive the following asymptotic representation with expansions where κ = n + 1 2 (α + 1), We assume that τ < 1. The quantity b is a function of x and follows from the relation where and The relation in (125) can be used for x ∈ [x 1 , x 2 ], in which case b ≥ 1 2 τ . In this interval the zeros of L , B 1 = P R 3 + QW 3 192R 3 W 4 (τ 2 − 1) , P = 4(2τ 2 + 12b 2 )(1 − τ 2 ), where the relation between b and x is given in (125). We write a zero in terms of b in the form where j k is a zero of the Bessel function J α (z). We assume for ε an expansion in the form By expanding U (b) at b 0 we have Using the representation of U (b) given in (129), substituting the expansion of ε, those of A(b) and B(b) given in (123), and comparing equal powers of κ, we can obtain the coefficients b j of (133). The first coefficients are For example, when we take n = 100, α = 75, then we obtain for the first zero b 0 = 0.1504907582034649. With this value for b we find from (125) a first approximation x = 0.0231157462791716, with a relative error 2.45 × 10 −5 . We compute with this x and b = b 0 the coefficient b 2 and find from b ∼ b 0 + b 2 /κ 2 the value b = 0.1504905751793771. Again inverting