Orthogonal systems in L 2 (R), once implemented in spectral methods, enjoy a number of important advantages if their differentiation matrix is skew-symmetric and highly structured. Such systems, where the differentiation matrix is skewsymmetric, tridiagonal and irreducible, have been recently fully characterised. In this paper we go a step further, imposing the extra requirement of fast computation: specifically, that the first N coefficients of the expansion can be computed to high accuracy in O (N log 2 N) operations. We consider two settings, one approximating a function f directly in (−∞, ∞) and the other approximatingIn each setting we prove that there is a single family, parametrised by α, β > −1, of orthogonal systems with a skew-symmetric, tridiagonal, irreducible differentiation matrix and whose coefficients can be computed as Jacobi polynomial coefficients of a modified function. The four special cases where α, β = ±1/2 are of particular interest, since coefficients can be computed using fast sine and cosine transforms. Banded, Toeplitz-plus-Hankel multiplication operators are also possible for representing variable coefficients in a spectral method. In Fourier space these orthogonal systems are related to an apparently new generalisation of the Carlitz polynomials.

Abstract. Many standard conversion matrices between coefficients in classical orthogonal polynomial expansions can be decomposed using diagonally-scaled Hadamard products involving Toeplitz and Hankel matrices. This allows us to derive O(N (log N ) 2 ) algorithms, based on the fast Fourier transform, for converting coefficients of a degree N polynomial in one polynomial basis to coefficients in another. Numerical results show that this approach is competitive with state-of-the-art techniques, requires no precomputational cost, can be implemented in a handful of lines of code, and is easily adapted to extended precision arithmetic.

In this paper we explore orthogonal systems in L2(R) which give rise to a skew-Hermitian, tridiagonal differentiation matrix. Surprisingly, allowing the differentiation matrix to be complex leads to a particular family of rational orthogonal functions with favourable properties: they form an orthonormal basis for L2(R), have a simple explicit formulae as rational functions, can be manipulated easily and the expansion coefficients are equal to classical Fourier coefficients of a modified function, hence can be calculated rapidly. We show that this family of functions is essentially the only orthonormal basis possessing a differentiation matrix of the above form and whose coefficients are equal to classical Fourier coefficients of a modified function though a monotone, differentiable change of variables. Examples of other orthogonal bases with skew-Hermitian, tridiagonal differentiation matrices are discussed as well.

The barycentric interpolation formula defines a stable algorithm for evaluation at points in [−1, 1] of polynomial interpolants through data on Chebyshev grids. Here it is shown that for evaluation at points in the complex plane outside [−1, 1], the algorithm becomes unstable and should be replaced by the alternative modified Lagrange or "first barycentric" formula dating to Jacobi in 1825. This difference in stability confirms the theory published by N. J. Higham in 2004 (IMA J. Numer. Anal., v. 24) and has practical consequences for computation with rational functions.

We introduce an algorithm for the least squares solution of a rectangular linear system Ax = b, in which A may be arbitrarily ill-conditioned. We assume that a complementary matrix Z is known such that A − AZ * A is numerically low rank. Loosely speaking, Z * acts like a generalized inverse of A up to a numerically low rank error. We give several examples of (A, Z) combinations in function approximation, where we can achieve high-order approximations in a number of nonstandard settings: the approximation of functions on domains with irregular shapes, weighted least squares problems with highly skewed weights, and the spectral approximation of functions with localized singularities. The algorithm is most efficient when A and Z * have fast matrix-vector multiplication and when the numerical rank of A − AZ * A is small.