Clenshaw's Recurrence Research Articles

Using parity to accelerate Hermite function computations: Zeros of truncated Hermite series, Gaussian quadrature and Clenshaw summation

Although Hermite functions have been studied for over a century and have been useful for analytical and numerical solutions in a myriad of areas, the theory of Hermite functions has gaps. This article is a unified treatment of all the operations—quadrature, summation, differentiation, and rootfinding—that can be accelerated by exploiting parity. Any function u(y) can be decomposed into its parts that are symmetric and antisymmetric with respect to the origin. Suppose that the quadrature on y∈[−∞,∞] is symmetric in the sense that if yj is a quadrature abscissa with weight wj, then −yj is also a quadrature point with weight wj. The number of multiplications is halved by evaluating the quadrature as ∫−∞∞u(y)dy≈w0u(0)+∑n=1Mwnu(yn)+u(−yn). Parity is equally useful in computing the zeros, maxima and minima of a truncated Hermite series of degree N for the important special cases that the series terms are all either symmetric or antisymmetric. The zeros and critical points are the eigenvalues of a companion matrix whose dimension is N/2 instead of N. In addition, for Hermite functions, we show that parity exploitation halves the dimension of the Jacobi matrix (a special case of the companion matrix) whose eigenvalues are the abcissas of Hermite–Gauss quadrature. The number of floating point operations for the recursion for the weights can likewise be halved. Lastly, the same is true for Chenshaw summation of a Hermite series.

High Order and Degree Geopotential and Derivatives Computation Based on the Clenshaw Summation

An artificial satellite subject to the attraction of the Earth is disturbed due to nonsphericaldistribution and non-symmetrical Earth mass. This uneven distribution of mass is expressed by the socalledspherical harmonic coefficients of the Earth potential. For a faster computation, the accelerationderived from the potential is obtained by a series expansion in terms of these harmonics, the fullynormalized Legendre polynomials and their derivatives, and several recursions associated with thelongitude, geocentric latitude and altitude of the center of mass of the satellite. This paper analyzesthe detailed aspects of disturbances in artificial satellites, related with the modeling of the Earth'sgravitational potential as well as numerical implementation of a recursive algorithm to calculate theacceleration of the geopotential based on the Clenshaw summation. In general, one uses recursiveequations of high degree and order to calculate the Legendre polynomials in order to obtain fasterprocessing and numerical accuracy. However, the recursions can yield numerical errors at each stepof the recursion so that higher orders and degrees of harmonics, the accumulated error may be quitepronounced. The computational implementation of the algorithm is carried out by a PC computer.With the implementation of this algorithm it is possible to calculate the geopotential accelerationfor different orbits and different situations. Such approach aims at mitigating the numerical problemsarising from use of extended series expansion when computing recursively the Legendre polynomials.Once the favorable numerical properties are proven, the algorithm can be used in the solution ofpractical problems of orbital space mechanics and for the Brazilian Space Mission.

Accurate Evaluation of Polynomials in Legendre Basis

This paper presents a compensated algorithm for accurate evaluation of a polynomial in Legendre basis. Since the coefficients of the evaluated polynomial are fractions, we propose to store these coefficients in two floating point numbers, such as double-double format, to reduce the effect of the coefficients’ perturbation. The proposed algorithm is obtained by applying error-free transformation to improve the Clenshaw algorithm. It can yield a full working precision accuracy for the ill-conditioned polynomial evaluation. Forward error analysis and numerical experiments illustrate the accuracy and efficiency of the algorithm.

A General Condition Number for Polynomials

This paper presents a generic condition number for polynomials that is useful for polynomial evaluation of a finite series of polynomial basis defined by means of a linear recurrence. This expression extends the classical one for the power and Bernstein bases, but it also provides us a general framework for all the families of orthogonal polynomials like Chebyshev, Legendre, Gegenbauer, Jacobi, and Sobolev orthogonal polynomial bases. The standard algorithm for the evaluation of finite series in any of these polynomial bases is the extended Clenshaw algorithm. The use of this new condition number permits us to give a general theorem about the forward error for that evaluation algorithm. A running-error bound of the extended algorithm is also presented and all the bounds are compared in several numerical examples.

Accurate evaluation of a polynomial in Chebyshev form

This paper presents a compensated algorithm to accurately evaluate a polynomial expressed in Chebyshev basis of the first and second kind with floating-point coefficients. The principle is to apply error-free transformations to improve the traditional Clenshaw algorithm. The new algorithm is as accurate as the Clenshaw algorithm performed in twice the working precision. Forward error analysis and numerical experiments illustrate the accuracy and properties of the proposed algorithm.

Accelerations of Zhao's methods for the numerical inversion of Laplace transform

Zhao recently presented two accurate and universal approaches to the numerical inversion of the Laplace transform. Both approaches are based on irregularly spaced intervals. Compared with the conventional approach by Durbin, the computational burden of Zhao's approaches is significantly heavier. In this report, the author proposed refinements by which the efficiency of Zhao's approaches can be improved significantly. This is achieved by subdividing the entire integrating range into several sub-bands, but the sampling interval in each sub-band is constant. In this way, the computation is accelerated by either applying Clenshaw's recurrence or the Chirp-Z transform. The efficiency of both accelerating approaches is verified numerically by three examples, which makes Zhao's two approaches practicable. In addition, the history of applying the fast Fourier transform to the Laplace transform is reviewed briefly. Copyright © 2009 John Wiley & Sons, Ltd.

Methods of harmonic synthesis for global geopotential models and their first-, second- and third-order gradients

Four widely used algorithms for the computation of the Earth’s gravitational potential and its first-, second- and third-order gradients are examined: the traditional increasing degree recursion in associated Legendre functions and its variant based on the Clenshaw summation, plus the methods of Pines and Cunningham–Metris, which are free from the singularities that distinguish the first two methods at the geographic poles. All four methods are reorganized with the lumped coefficients approach, which in the cases of Pines and Cunningham–Metris requires a complete revision of the algorithms. The characteristics of the four methods are studied and described, and numerical tests are performed to assess and compare their precision, accuracy, and efficiency. In general the performance levels of all four codes exhibit large improvements over previously published versions. From the point of view of numerical precision, away from the geographic poles Clenshaw and Legendre offer an overall better quality. Furthermore, Pines and Cunningham–Metris are affected by an intrinsic loss of precision at the equator and suffer from additional deterioration when the gravity gradients components are rotated into the East-North-Up topocentric reference system.

Evaluation of methods for spherical harmonic synthesis of the gravitational potential and its gradients

This work is concerned with the comparison of four of the best-known methods for the computation of the gravitational potential and its gradients: the traditional formulation in terms of Associated Legendre Functions in spherical coordinates; the non-singular method of Pines; the algorithm developed by Cunningham and extended by Metris and collaborators; and a variant of the first method based on the Clenshaw summation formula. Extensive numerical tests in double and quadruple floating point precision have been performed in order to assess and compare the efficiency and precision of these algorithms. Results show that when properly optimized the algorithm of Clenshaw is the most efficient, closely followed by the traditional Legendre formulation. All four methods are characterized by a high level of precision, although care should be taken when approaching the geographic poles due to the singularities which affect the methods of Legendre and Clenshaw. The methods of Cunningham–Metris and Pines are both characterized by some loss of relative precision at the equator, which is inherent in the choice of the coordinate system.

Recursive computation of Tchebichef moment and its inverse transform

Tchebichef moment is a novel set of orthogonal moment applied in the fields of image analysis and pattern recognition. Less work has been made for the computation of Tchebichef moment and its inverse moment transform. In this paper, both a direct recursive algorithm and a compact algorithm are developed for the computation of Tchebichef moment. The effective recursive algorithm for inverse Tchebichef moment transform is also presented. Clenshaw's recurrence formula was used in this paper to transform kernels of the forward and inverse Tchebichef moment transform. There is no need for the proposed algorithms to compute the Tchebichef polynomial values. The approaches presented are more efficient compared with the straightforward methods, and particularly suitable for parallel VLSI implementation due to their regular and simple filter structures.

Black-box modeling of passive systems by rational function approximation

In this paper, a rational interpolation approach is used to approximate the transfer function of passive systems characterized by sampled data. Orthogonal polynomials are used to improve the numerical stability of the ill-conditioned Vandermonde-like interpolation matrix associated with the ordinary power series. First, the poles of the system are obtained by efficiently and accurately transforming the coefficients of the orthogonal polynomials to the ordinary power series using Clenshaw's recurrence algorithm. Then, the residues are solved in real or in complex conjugate pairs to insure a physically realizable system. Finally, the passivity of the system is enforced by applying certain constraints on the poles and residues of the system. The performances of the three most common orthogonal polynomials, Legendre and Chebyshev of the first and second kinds, are also compared to that of the power series.

Parallel Recursive Computation of the Inverse Legendre Moment Transforms for Signal and Image Reconstruction

The inverse Legendre moment transform plays an important role in the field of image analysis applications. A recursive algorithm based on Clenshaw's recurrence formula is derived for the fast computation of the inverse Legendre moment transform for signal and image reconstruction purposes. The reconstruction can be implemented effectively using recursive equations. The presented algorithm is particularly suitable for parallel very large-scale integrated circuit (VLSIC) implementation due to its regular, modular, and simple filter structure.

Computation of forward and inverse MDCT using clenshaw's recurrence formula

The forward and inverse modified discrete cosine transform (MDCT) are two of the most computationally intensive operations in the MPEG audio coding standard. In this paper, we used Clenshaw's recurrence formula to transform kernels of the MDCT and IMDCT of the general length. Efficient implementations of MDCT and IMDCT are obtained. The proposed regular structures are particularly suitable for parallel VLSI realization. Our solution requires less hardware, and we achieved significant savings for a number of operations compared with existing related systems.

A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions

Spherical harmonic expansions form partial sums of fully normalised associated Legendre functions (ALFs). However, when evaluated increasingly close to the poles, the ultra-high degree and order (e.g. 2700) ALFs range over thousands of orders of magnitude. This causes existing recursion techniques for computing values of individual ALFs and their derivatives to fail. A common solution in geodesy is to evaluate these expansions using Clenshaw's method, which does not compute individual ALFs or their derivatives. Straightforward numerical principles govern the stability of this technique. Elementary algebra is employed to illustrate how these principles are implemented in Clenshaw's method. It is also demonstrated how existing recursion algorithms for computing ALFs and their first derivatives are easily modified to incorporate these same numerical principles. These modified recursions yield scaled ALFs and first derivatives, which can then be combined using Horner's scheme to compute partial sums, complete to degree and order 2700, for all latitudes (except at the poles for first derivatives). This exceeds any previously published result. Numerical tests suggest that this new approach is at least as precise and efficient as Clenshaw's method. However, the principal strength of the new techniques lies in their simplicity of formulation and implementation, since this quality should simplify the task of extending the approach to other uses, such as spherical harmonic analysis.

Backward Stability of Clenshaw's Algorithm

We study numerical properties of Clenshaw's algorithm for summing the series w = ∑n = 0Nbnpn where pn satisfy the linear three-term recurrence relation. We prove that under natural assumptions Clenshaw's algorithm is backward stable with respect to the data bn, n = 0,N.

Error Analysis of Clenshaw's Algorithm for Evaluating Derivatives of a Polynomial

A forward rounding error analysis is presented for the extended Clenshaw algorithm due to Skrzipek for evaluating the derivatives of a polynomial expanded in terms of orthogonal polynomials. Reformulating in matrix notation the three-term recurrence relation satisfied by orthogonal polynomials facilitates the estimate of the rounding error for the m-th derivative, which is recursively estimated in terms of the one for the (m − 1)-th derivative. The rounding errors in an important case of Chebyshev polynomial are discussed in some detail.

Methods for designing efficient parallel-recursive filter structures for computing discrete transforms

Analytical and numerical approaches are presented for the design of first?order and second?order recursive digital filter banks for computing linear, discrete transforms. This subject has been studied extensively for the case of trigonometric transforms. The focus of this paper is on discrete polynomial transforms, and Clenshaw's recurrence formulae are used to design the second?order filters. The efficiency of the implementation for a specific transform is dependent upon the characteristics of recurrence relations for the transform basis vectors. Efficient implementations are derived for the discrete cosine transform and the inverse discrete Legendre transform from analytical expressions for basis vector recurrence relations. A numerical procedure is presented to search for the existence and parameters of an efficient implementation when analytical expressions for the basis vector recurrence relations are unknown.

Polynomial evaluation and associated polynomials

We show a connection between the Clenshaw algorithm for evaluating a polynomial \(q_n\), expanded in terms of a system of orthogonal polynomials, and special linear combinations of associated polynomials. These results enable us to get the derivatives of \(q_n(z)\) analogously to the Horner algorithm for evaluating polynomials in monomial representations. Furthermore we show how a polynomial \(\hat q_n(z)\) given in monomial (!) representation can be evaluated for \(z\in{\Bbb C}\) using the Clenshaw algorithm without complex arithmetic. From this we get a connection between zeros of polynomials expanded in terms of Chebyshev polynomials and the corresponding polynomials in monomial representation with the same coefficients.

Fast algorithms for discrete polynomial transforms

Consider the Vandermonde-like matrix P := ( P k ( cos ⁡ j π N ) ) j , k = 0 N {\mathbf {P}}:=(P_k(\cos \frac {j\pi }{N}))_{j,k=0}^N , where the polynomials P k P_k satisfy a three-term recurrence relation. If P k P_k are the Chebyshev polynomials T k T_k , then P {\mathbf {P}} coincides with C N + 1 := ( cos ⁡ j k π N ) j , k = 0 N {\mathbf {C}}_{N+1}:= (\cos \frac {jk\pi }{N})_{j,k=0}^N . This paper presents a new fast algorithm for the computation of the matrix-vector product P a {\mathbf {Pa}} in O ( N log 2 ⁡ N ) O(N \log ^2N) arithmetical operations. The algorithm divides into a fast transform which replaces P a {\mathbf {Pa}} with C N + 1 a ~ {\mathbf {C}}_{N+1} {\mathbf {\tilde a}} and a subsequent fast cosine transform. The first and central part of the algorithm is realized by a straightforward cascade summation based on properties of associated polynomials and by fast polynomial multiplications. Numerical tests demonstrate that our fast polynomial transform realizes P a {\mathbf {Pa}} with almost the same precision as the Clenshaw algorithm, but is much faster for N ≥ 128 N\ge 128 .

Computation of discrete cosine transform using Clenshaw's recurrence formula

Clenshaw's recurrence formula is used to derive recursive algorithms for the discrete cosine transform (DCT) and the inverse discrete cosine transform (IDCT). The recursive DCT algorithm presented requires one fewer delay element per coefficient and one fewer multiply operation per coefficient compared with two other proposed methods. Clenshaw's recurrence formula provides a unified development for the recursive DCT and IDCT algorithms. The recursive algorithms apply to arbitrary length algorithms and are appropriate for VLSI implementation. >

An efficient algorithm for choosing scattering directions in Monte Carlo work with arbitrary phase functions

This paper describes an efficient Monte Carlo algorithm for choosing a new direction of a photon after a scattering interaction. The algorithm chooses a scattering angle by linear interpolation in a table of the inverse cumulative scattering probability. A Legendre expansion of the phase function makes it easy to apply Clenshaw's algorithm to build the interpolation table. The points in the table are close enough together that linear interpolation is accurate. With a table of 100,000 entries, we can keep the absolute and relative errors in matching the probability distribution below 10 −5.