1 Introduction

When automatic computers began to appear in the 1950s various confident, and often incorrect, predictions were made concerning the impact of these devices on applied mathematics, science and engineering. One of these predictions was that the need for special functions, or higher transcendental functions (as they are also known), would disappear entirely. This was based on the observation that the main use of these functions in those days was to approximate the solutions of classical partial differential equations: with automatic computers it would become possible to solve these equations by direct numerical methods. This observation is indeed correct; nevertheless, a perusal of current computational journals in the sciences reveals a persistent need for numerical algorithms to generate Airy functions, Bessel functions, Coulomb wave functions, error functions and exponential integrals--to name but a few of the classical special functions. Equally significantly, the National Bureau of Standards' Handbook of Mathematical Functions [AS64]³continues to be one of the best-selling mathematical books of all time⁴.

The purpose of the present paper is to provide some assistance to those mathematicians, engineers, scientists, and statisticians who discover that they need to generate numerical values of the special functions in the course of solving their problems. ``Generate'' is the operative word here: we are thinking primarily of either software or numerical approximations that can be programmed fairly easily. Numerical tables are not covered in this survey. Furthermore, for the most part we shall concentrate on the functions themselves; only in certain cases do we include, for example, zeros, inverse functions or indefinite integrals. Elementary functions, also, are excluded⁵. Lastly, we believe that the majority of readers would prefer us to emphasize the more useful algorithms rather than make an attempt to be encyclopedic: algorithms or approximations that have clearly been superseded are omitted.

We identify three stages in the development of computational procedures for the special functions:

1. Derivation of relevant mathematical properties.

2. Development of numerical approximations and algorithms.

3. Construction and testing of robust software.

Included in Stage (i) are asymptotic expansions, continued fractions, difference and differential equations, functional identities, integral representations, and Taylor-series expansions. Included in Stage (ii) are expansions in series of Chebyshev polynomials (``Chebyshev series''), minimax polynomial and rational approximations, Padé approximations, numerical quadrature, and numerical solution of difference and differential equations. In this paper the emphasis will be on Stages (ii) and (iii), but in §2 we supply some general references for Stage (i).

In §3 we make a general survey of software libraries and packages⁶that include collections of special functions.

In §§4 and 5 the functions are treated individually. We list software that is already available and readily programmable numerical approximations. In §6 we also include references to articles (or books) that may be useful for the testing or comparison of existing software, or in the construction of new libraries. We do not attempt what would be a herculean task of testing and comparing everything that is available. Our philosophy in this survey has to be that of caveat emptor: no algorithm, approximation or package that we mention should be relied upon to produce accurate output in the absence of evidence of independent and systematic checks.

As in the progress of other branches of numerical analysis, procedures used to evaluate special functions are influenced heavily by the computing equipment available at the time. In the era of desk-calculating machines the medium was a printed table of numerical values of the wanted function, or functions, generally for equispaced values of the arguments. Nontabular values were calculated by means of Lagrange's interpolation formula or central-difference interpolation formulas [Fox56]. In other words, local polynomial approximations were employed. These interpolation procedures were reasonably successful for functions of a single variable, but two-dimensional interpolation on desk-calculating machines was often a laborious computation that was prone to error. The daunting task that faced a (human) computer is summed up in the following quotation from the Introduction to Karl Pearson's tables of the incomplete gamma function [Pea22]:

``As a matter of fact, supposing the use of a machine, which every modern computer has at his command, no interpolation suggested ought to take more than an hour's work and many much less. If the user of these tables groans under that hour, let him compute de novo a function value, say $I(6.86877, 47.1813)$ --including of course $\Gamma (48.1813)$--to seven-figure accuracy, and when he has completed the task, we believe his feelings towards those who have provided him with these tables will be very sensibly modified.''

With electronic computers, the number of arithmetic operations that could be contemplated for the generation of a single function value increased considerably. In consequence, high-degree global approximations appeared for functions of a single variable in the form of minimax polynomial or rational approximations, or truncated Chebyshev-series expansions; see, for example, [Cle62,HCL+68] and [Luk77b].

Chebyshev series in two dimensions also became feasible [CP66,Luk71a,Luk71b,Luk72a]. However, because of their more complicated asymptotic behavior, special functions of two variables cannot be covered comprehensively simply by use of polynomial or rational approximations or Chebyshev series. The situation is cloudier still when the variables or parameters are complex, or of course when they are more than two in number. For this reason, to achieve maximum speed, a comprehensive software package for generating a function of two (or more) variables typically employs several different algorithms in addition to, or quite commonly in place of, polynomial or rational approximations or Chebyshev series. The construction and testing of such a package invariably entails prodigious effort.

Computers continue to increase in sophistication and power; in consequence we should expect further changes in the algorithms used to generate the special functions. So far, the potential offered by the introduction of vector and parallel computing machines has not been exploited to any great extent. It might well lead to simplifications in the algorithms needed for many functions, as well as to an increase in execution speeds. We refer again to this possibility in the concluding section (§7).

In assembling the bibliography of this paper we have been assisted by the references collected and classified by the late Dr. L. W. Fullerton in his 1980 Bell Laboratories report [Ful80], by access to Dr. N. M. Temme's private collection of references, and by GAMS, the Guide to Available Mathematical Software prepared by the National Institute of Science and Technology [BHKS90]. Accessible at http://gams.nist.gov/, GAMS is a convenient, free source for documentation and nonproprietary source code.

We have searched through issues of the following journals of the past twenty-five years for relevant references:

Applied Statistics (Appl. Statist.), Association for Computing Machinery Transactions on Mathematical Software (ACM Trans. Math. Software), BIT, Collected Algorithms from the Association for Computing Machinery (CALGO), Communications of the Association for Computing Machinery (Comm. ACM), Computer Physics Communications (Comput. Phys. Comm.), Computing, Journal of Computational and Applied Mathematics (J. Comput. Appl. Math.), Journal of Computational Physics (J. Comput. Phys.), Mathematical Reviews (Math. Rev.), Mathematics of Computation (Math. Comp.), Numerische Mathematik (Numer. Math.), U.S.S.R. Computational Mathematics and Mathematical Physics (U.S.S.R. Comput. Math. and Math. Phys.)⁷, Zeitschrift für Angewandte Mathematik und Physik (Z. Angew. Math. Phys.), Zentralblatt für Mathematik und ihre Grenzgebiete (Zbl.).

However, because of the sheer volume and diversity of publications on special functions it is almost inevitable that we have overlooked some useful algorithms and important articles. In this event we tender, in advance, our apologies to the authors.