Press here to get the full document in PostScript format.

Press here to get this subdocument in PostScript format.

Numerical Evaluation of Special Functions

D. W. Lozier and F. W. J. Olver

3.3. Comprehensive Libraries .

When new algorithms are developed they tend to appear first as subroutines in software packages (§ 3.1 above). Later they may be assimilated into more complete software products such as intermediate libraries (§ 3.2 above). Even more useful are comprehensive libraries that integrate evaluation of special functions with other essential elements of numerical computing and offer additional advantages such as uniform documentation, style of usage, and handling of error conditions. Corrections and improvements, particularly in orienting the programming toward particular computer architectures, are often made.

3.3.1. CERN Library.

The European Laboratory for Particle Physics maintains a comprehensive software library [ CER93] , mostly in Fortran but with a few routines in assembly language, in support of high-energy physics research. The coverage of special functions includes the error function of real and complex argument; the Dawson and Fresnel integrals; exponential, sine, cosine and arctangent integrals; the gamma and digamma functions of real and complex argument; incomplete gamma function; real dilogarithm and complex generalized polylogarithms; Bessel functions of real argument and orders 0, , , , , and 1; Bessel functions of real order and real or complex argument; Bessel functions of complex order and argument; zeros of the Bessel functions J and Y and of their derivatives; Coulomb wave functions of complex order, argument and parameter; complete and incomplete elliptic integrals; Jacobi's elliptic functions (real and complex); Jacobi's theta functions (real); Bose-Einstein and Fermi-Dirac integrals ; Legendre and associated Legendre functions; conical functions of the first kind; Struve functions. The library is distributed, with some restrictions, to organizations outside CERN.

3.3.2. IMSL Library.

International Mathematical and Statistical Libraries was incorporated in 1970 ``with the intent of providing high-quality, supported Fortran subroutine libraries in mathematics and statistics'' [ Air84] . In its first ten years it produced libraries tailored to twelve different computer lines, providing an alternative to manufacturer-supplied libraries. Currently the company offers a wide range of products for large-scale scientific computing. At the center of its product line is the IMSL numerical subroutine library for mathematics and statistics, which includes an extensive coverage of real and complex special functions [ IMS91] ; this reference includes a GAMS index [ BHK91] and a KWIC (keyword in context) index. The library is optimized, vectorized and parallelized where appropriate, depending on the target computer architecture, but it contains no vector or parallel support for special functions.

A subset of the IMSL library is offered also as a C library. This is an essential component of a powerful interactive system (§ 3.4 below) which has the capability of providing graphical and numerical computing with very large data sets. Optionally, Maple (§ 3.4.3 below) can be incorporated to provide for symbolic computing.

3.3.3. NAG Library.

An overview of the development, structure and contents of the NAG Fortran library [ NAG91] is given in [ FP84] . After originating as a cooperative project among several British computing centers in 1970 with the purpose of providing ``a balanced, general-purpose numerical algorithms library," the Numerical Algorithms Group formed a not-for-profit company in 1976 to provide for the wider distribution of the library. The library is organized around the ACM modification of the SHARE classification system (see § 3.1.1 above) and is available for a wide cross-section of computing systems. A KWIC (keyword in context) index and an index in the GAMS classification scheme [ BHK91] are provided in the library documentation.

Subsets of the full library are available in Ada, Algol 68, C, Fortran 90 and Pascal. NAG is actively developing and marketing an interactive system (§ 3.4 below) that integrates most of the numerical power of the full NAG library with online symbolic and graphical capabilities.

3.3.4. NSWC Library.

In 1976 the Naval Surface Warfare Center, Dahlgren, Virginia, began development of ``a [Fortran] library of general purpose subroutines that would provide a basic computational ability in a variety of mathematical activities'' [ Mor93] . The design goals stressed reliability, transportability, efficiency, ease of use, and generality. In 1993 the library contained 576 user-level subroutines, including ones in real arithmetic for the error function and its inverse; Dawson's integral and Fresnel integrals; exponential, sine and cosine integrals; gamma, psi and polygamma functions; dilogarithm; incomplete gamma function and its inverse; incomplete beta function; complete and incomplete elliptic integrals; Jacobi's and Weierstrass' elliptic functions; Bessel functions of real argument and order. It also contained Airy functions and complete elliptic integrals of complex argument, and Bessel functions of complex argument and integer or complex order.

3.3.5. NUMAL Library.

In 1973 the Mathematisch Centrum, Amsterdam, introduced this library of numerical procedures in Algol-60 with `` the aim ... to provide Algol-60 programmers with a high-level numerical library which contains a set of validated numerical procedures together with supporting documentation'' [ Hem81] . In 1979 it contained approximately 430 subroutines, including ones for the exponential, sine and cosine integrals; gamma function; incomplete gamma and beta functions; error and inverse error functions; Fresnel integrals; modified and unmodified Bessel functions of integer, half-integer or real order; Airy functions. All subroutines are for real variables.

3.3.6. Numerical Recipes.

This partly pedagogical series of books offers ``for each topic considered, a certain amount of general discussion, a certain amount of analytical mathematics, a certain amount of discussion of algorithmics, and (most important) actual implementations of these ideas in the form of working computer routines" [ PTVF92] . Besides being listed fully in the text, the computer routines are available for purchase on diskette or under a variety of licensing arrangements, one of which is tailored to the needs of classroom instructors. Example books with test programs and diskettes are available also. Standard fields of numerical computation are covered, with approximation of functions and evaluation of special functions included. The book is published in four versions with the software coded in Basic, C, Fortran or Pascal; another volume for Basic is [ Spr91] .

3.3.7. NUMPAC Library.

The Nagoya University Mathematical Package is used widely in Japan. It is a comprehensive Fortran library oriented toward Japanese computers, including vector supercomputers. Coverage includes Airy functions; error and inverse error functions; Dawson and Fresnel integrals; exponential, sine and cosine integrals; complex gamma function; digamma function; dilogarithm; Riemann's zeta function; Bessel functions of integer or real order and real or complex argument; zeros and integrals of Bessel functions; complete and incomplete elliptic integrals; Jacobi's elliptic functions; incomplete beta and gamma functions; Legendre polynomials and associated Legendre functions; classical orthogonal polynomials; Struve functions; Abramowitz, Debye and elliptic theta functions; solutions of the Blasius and Thomas-Fermi equations. Information can be obtained from Ichizo Ninomiya, Chubu University, Kasugai, Aichi, 487 Japan or Yasuyo Hatano, Chukyo University, Yagoto, Nagoya, 466 Japan.

3.3.8. PORT Library.

This library [ FHS78b,Fox84] is mentioned here because it provides a framework [ FHS78a] for constructing portable Fortran libraries that has proven its utility. The framework supplies computer arithmetic parameters via Fortran function calls. Algorithms are coded so as to be valid for a range of values of the arithmetic parameters; actual values are substituted at run time. The PORT framework is used, for example, in the SLATEC library (§ 3.3.10 below). It is particularly valuable in the special function routines because of their sensitivity to precision, underflow and overflow. The PORT framework is available as ACM Algorithm 528 (§ 3.1.1 above).

3.3.9. Scientific Desk Library.

C. Abaci offers the following products: (i) the Scientific Desk Library, a Fortran-based collection of numerical software; (ii) the Scientific Desk Analysis System, an interactive system (§ 3.4 below); (iii) software produced by others, including the ACM algorithms (§ 3.1.1 above). The library is available in object code for personal computers under a variety of Fortran compilers and in Fortran source code for other computers. The Analysis System, which is strongly oriented toward statistics, simplifies the programming burden and provides for simple graphical output. C. Abaci distributes the SPECFUN collection [ Cod93b] of Fortran programs for special functions and the ELEFUNT, INTFUNT and CELEFUNT tests [ Cod93a,CW80] for elementary functions. Inquire at C. Abaci, Inc., P. O. Box 2626, Raleigh, NC 27602.

3.3.10. SLATEC Library.

The acronym stands for Sandia, Los Alamos, Air Force Weapons Laboratory Technical Exchange Committeegif, formed in 1974 to ``foster the exchange of technical information among the three computing departments''. In 1977 a subcommittee undertook the development of a complete, noncommercial Fortran library for numerical supercomputing [ Buz84] . The primary motivation was that the suppliers of commercial libraries regarded the supercomputing market as too small. The library subcommittee admitted subsequently five additional U. S. Government agencies (the Lawrence Livermore, Oak Ridge and Sandia Livermore National Laboratories, the National Energy Supercomputer Center at Lawrence Livermore, and the National Institute of Standards and Technology). SLATEC Version 1.0 appeared in 1981. Version 4.0, the third major revision and expansion, was released in December 1992. The initial coverage of special functions coincided with FNLIB [ Ful77] , FUNPACK [ Cod75,Cod84a] and AMOSLIB [ AD79] . Subroutines from [ ADW77a,ADW77b,Amo80a,Amo83a,Amo83b,Amo86,CN81,LS81,OS83] were added later. Available from Energy Science and Technology Center, P. O. Box 1020, Oak Ridge, TN 37831 and from Netlib (§ 3.1.1 above).


D. E. Amos and S. L. Daniel, AMOSLIB, A special function library, version 9/77, Report 77-1390, Sandia Laboratories, August 1979.

D. E. Amos, S. L. Daniel, and M. K. Weston, Algorithm 511. CDC 6600 subroutines IBESS and JBESS for Bessel functions and , , , ACM Trans. Math. Software 3 (1977), 93--95, for erratum see same journal v. 4 (1978), p. 411.

D. E. Amos, S. L. Daniel, and M. K. Weston, CDC 6600 subroutines IBESS and JBESS for Bessel functions and , , , ACM Trans. Math. Software 3 (1977), 76--92.

T. J. Aird, The IMSL library, Sources and Development of Mathematical Software (W. R. Cowell, ed.), Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1984, pp. 264--301.

D. E. Amos, Algorithm 556. Exponential integrals, ACM Trans. Math. Software 6 (1980), 420--428, for remark see same journal v. 9 (1983), p. 525.

D. E. Amos, Algorithm 609. A portable Fortran subroutine for the Bickley functions , ACM Trans. Math. Software 9 (1983), 480--493.

D. E. Amos, Algorithm 610. A portable Fortran subroutine for derivatives of the psi function, ACM Trans. Math. Software 9 (1983), 494--502.

D. E. Amos, Algorithm 644. A portable package for Bessel functions of a complex argument and nonnegative order, ACM Trans. Math. Software 12 (1986), 265--273, for remark see same journal v. 16 (1990), p. 404.

R. F. Boisvert, S. E. Howe, and D. K. Kahaner, The guide to available mathematical software problem classification system, Comm. Statist. B---Simulation Comput. 20 (1991), 811--842.

B. L. Buzbee, The SLATEC common mathematical library, Sources and Development of Mathematical Software (W. R. Cowell, ed.), Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1984, pp. 302--320.

CERNLIB short writeups, CERN Program Library Office, CERN--CN Division, CH--1211 Geneva 23, Switzerland, 1993, electronic mail address is cernlib@

B. C. Carlson and E. M. Notis, Algorithm 577. Algorithms for incomplete elliptic integrals, ACM Trans. Math. Software 7 (1981), 398--403.

W. J. Cody, The FUNPACK package of special function routines, ACM Trans. Math. Software 1 (1975), 13--25.

W. J. Cody, FUNPACK---A package of special function routines, Sources and Development of Mathematical Software (W. R. Cowell, ed.), Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1984, pp. 49--67.

W. J. Cody, Algorithm 714. CELEFUNT : A portable test package for complex elementary functions, ACM Trans. Math. Software 19 (1993), 1--21.

W. J. Cody, Algorithm 715. SPECFUN : A portable Fortran package of special function routines and test drivers, ACM Trans. Math. Software 19 (1993), 22--32.

W. J. Cody and W. Waite, Software manual for the elementary functions, Prentice Hall, Englewood Cliffs, New Jersey, 1980.

P. A. Fox, A. D. Hall, and N. L. Schryer, Algorithm 528. Framework for a portable library, ACM Trans. Math. Software 4 (1978), 177--188.

P. A. Fox, A. D. Hall, and N. L. Schryer, The PORT mathematical subroutine library, ACM Trans. Math. Software 4 (1978), 104--126.

P. A. Fox, The PORT mathematical subroutine library, Sources and Development of Mathematical Software (W. R. Cowell, ed.), Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1984, pp. 346--374.

B. Ford and J. C. T. Pool, The evolving NAG library service, Sources and Development of Mathematical Software (W. R. Cowell, ed.), Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1984, pp. 375--397.

L. W. Fullerton, Portable special function routines, Lecture Notes in Computer Science 57: Portability of Numerical Software (Oak Brook 1976) (W. R. Cowell, ed.), Springer-Verlag, Berlin, 1977, pp. 452--483.

P. W. Hemker (ed.), NUMAL numerical procedures in Algol 60, Mathematisch Centrum Syllabus 47.1--7, Centrum voor Wiskunde en Informatica, Kruislaan 413, 1098 SJ, Amsterdam, 1981.

MATH/LIBRARY special functions, version 2.0, IMSL, Inc., 2500 CityWest Boulevard, Houston, Texas 77042, September 1991, current address is Visual Numerics, Inc., 9990 Richmond Ave, Suite 400, Houston, Texas 77042-4548.

D. W. Lozier and J. M. Smith, Algorithm 567. Extended-range arithmetic and normalized Legendre polynomials, ACM Trans. Math. Software 7 (1981), 141--146.

A. H. Morris, Jr., NSWC library of mathematics subroutines, NSWCDD/TR--92/425, Naval Surface Warfare Center, Dahlgren Division, Dahlgren, Virginia, 22448, January 1993.

NAG Fortran library manual, mark 15, vol. 10, NAG Inc., 1400 Opus Place, Suite 200, Downers Grove, Illinois 60515-5702, June 1991.

F. W. J. Olver and J. M. Smith, Associated Legendre functions on the cut, J. Comput. Phys. 51 (1983), 502--518.

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical recipes. The art of scientific computing, second ed., Cambridge University Press, 1992, diskettes and example books available. Editions exist in Basic (1991), C (1992), Fortran (1992), Macintosh Fortran (1988) and Pascal (1989).

J. C. Sprott, Numerical recipes : routines and examples in BASIC, Cambridge University Press, 1991.


This document is an excerpt from the current hypertext version of an article that appeared in Walter Gautschi (ed.), Mathematics of Computation 1943--1993: A Half-Century of Computational Mathematics, Proceedings of Symposia in Applied Mathematics 48, American Mathematical Society, Providence, RI 02940, 1994. The symposium was held at the University of British Columbia August 9--13, 1993, in honor of the fiftieth anniversary of the journal Mathematics of Computation.

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:

Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742

E-mail address:

The research of the second author has been supported by NSF Grant CCR 89-14933.

1991 Mathematics Subject Classification. Primary 65D20; Secondary 33-00.

The Air Force Weapons Laboratory has been renamed the Phillips Laboratory.

Daniel W Lozier
Fri Apr 7 13:39:24 EDT 1995