continues to be one of the best-selling mathematical
books of all time.
Press here to get the full document in PostScript format.
Press here to get this subdocument in PostScript format.
D. W. Lozier and F. W. J. Olver
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:
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 
---including of course 
---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]
and also available online.
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.
68
--function, H. M.
Stationery Office, London, 1922, reissued by Biometrika, 1934.
The original abstract follows.
Higher transcendental functions continue to play varied and important roles in investigations by engineers, mathematicians, scientists and statisticians. The purpose of this paper is to assist in locating useful approximations and software for the numerical generation of these functions, and to offer some suggestions for future developments in this field.
Applied and Computational Mathematics Division, National Institute of Standards and Technology, Gaithersburg, Md 20899
E-mail address: dlozier@nist.gov
Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742
E-mail address: olver@bessel.umd.edu
The research of the second author has been supported by NSF Grant CCR 89-14933.
1991 Mathematics Subject Classification. Primary 65D20; Secondary 33-00.