A Proposal to Add Mathematical Special Functions to the C++ Standard Library

Document number:   WG21/N1422 = J16/03-0004
Date:   February 24, 2003
Project:   Programming Language C++
Reference:   ISO/IEC IS 14882:1998(E)
Reply to:   Walter E. Brown <wb@fnal.gov>
  Fermi National Accelerator Laboratory
  Batavia, Illinois, USA


Contents

I. Background and Motivation
II. Impact On the C++ Standard
III. Design Decisions
IV. Proposed Text
V. Relationship to Earlier Proposal
VI. Acknowledgements
VII. References

I. Background and Motivation

Why is this important? What kinds of problems does it address, and what kinds of programmers is it intended to support? Is it based on existing practice?

Compared to C++ [ISO:14882], C99 [ISO:9899] provides an extended <math.h> header and library. Among the additions introduced by C99 are selected mathematical functions from categories of particular interest to the numerical computing communities: exponential and logarithmic functions, circular and hyperbolic functions, and special functions.

In particular, C99 specifies the following traditional and extended functions of numerical interest in <math.h>, each with variants to accomodate arguments of types float, double, and long double:

All these functions either (a) are already part of the C++ standard library, or (b) have already been proposed [Plauger] and discussed [Dawes, Item 1] for incorporation into the forthcoming C++ Library Technical Report [Josuttis]. We here propose to augment the C++ standard library with additional functions from the above Special Functions category.

Mathematical Special Functions are appropriate for this TR because they plug an obvious hole ("Filling Gaps," as [Austern] phrases it) in the existing standard library. While these functions are clearly numerical in nature and will likely be most heavily used by the scientific and engineering communities, other communities of programmers also have needs, ranging from frequent to intermittent, for these functions.

This Special Functions proposal additionally falls into the "Standards Coordination and Infrastructure" categories identified in [Austern] as targets for the TR, for this proposal is based on an existing standard, Quantities and units [ISO:31] Part 11: Mathematical signs and symbols for use in the physical sciences and technology. We draw particular attention to the tables constituting ISO:31-11 paragraphs 8 ("Exponential and logarithmic functions"), 9 ("Circular and hyperbolic functions"), 10 ("Complex numbers"), and 14 ("Special functions").

There is a long history of implementation experience with these functions, as evidenced, for example, by section C of the Fortran-based SLATEC Common Mathematical Library [SLATEC]. Today, special functions constitute important subsets of such well-respected add-on libraries as the NAG C Library, [NAG], the IMSL C Numerical Library [IMSL], and the GNU Scientific Library [GSL]. Further, some standard C libraries such as the SGI C library [SGI] and the GNU C library [GNU C] also provide, as extensions, a few of the proposed special functions. The benefits of incorporating them into the C++ Standard Library include predictability of interface and behavior across a broad spectrum of implementations, leading to improved portability and interoperability for applications that make use of these functions.

Finally, we believe that adoption of this proposal would send a clear message to the various numeric computing communities that, contrary to significant popular belief within these communities, C++ is an eminently suitable programming language for their problem domain, too.

II. Impact On the C++ Standard

What does it depend on, and what depends on it? Is it a pure extension, or does it require changes to standard components? Does it require core language changes?

This proposal is a pure extension. It does not require any changes in the core language. It does not require changes to any standard classes or functions. It does not require changes to any of the standard requirement tables. All the functions are mathematically well-understood, all have proven their utility in practice over a considerable period of time, and all have been previously implemented in C and C++.

This proposal does not depend on any other C++ library extensions. This proposal potentially overlaps slightly with another proposal [Plauger] that would incorporate the bulk of C99's library additions into C++. The potential commonality between the two proposals is, however, limited to a rather tiny part of <math.h> that is essentially identical in the two proposals; see section V for details.

III. Design Decisions

Why did you choose the specific design that you did? What alternatives did you consider, and what are the tradeoffs? What are the consequences of your choice, for users and implementors? What decisions are left up to implementors? If there are any similar libraries in use, how do their design decisions compare to yours?

A. How to Package the Additional Declarations?

Following the precedent set by C99, this proposal recommends that declarations for all the proposed special functions be incorporated into <math.h> and thence extended into <cmath> in the obvious way.

An alternative design would present these additional declarations in a new header. Obvious names for this header, e.g., <special_functions.h>, seem unwieldy, and no suitably descriptive shorter names have come to mind. Further, it seems likely that implemention of some of the special functions can make advantageous use of extant functionality in <math.h> and so it seemed reasonable to avoid the separation.

B. Which Special Functions to Incorporate?

Because the set of mathematical functions that can be considered Special Functions is potentially unbounded, we considered several options in selecting our list of candidates for standardization.

This proposal recommends adoption of the list of Special Functions specified in the ISO standard Quantities and units [ISO:31, Part 11, paragraph 14]. This list has already received international scrutiny and endorsement via the standardization process. Further consultations with respected scientific colleagues have confirmed that these Special Functions would, if incorporated into the C++ standard library, constitute a significant contribution to the numeric community.

A second possibility was the adoption of all the Special Functions listed in the (exhaustive!) Handbook of Mathematical Functions [Abramowitz & Stegun], the generally-accepted standard reference for this domain. However, a careful inspection of this work's extensive contents strongly suggests that its scope may be overly broad. More importantly, many of the listed functions appear to be very difficult to implement. (Indeed, a colleague suggested, not entirely in jest, that "several Ph.D. dissertations could result" from such efforts!)

We considered, third, adopting a list taken from an existing library in this domain. The Special Functions portion of the GNU Scientific Library [GSL], for example, seemed to present a reasonable set for consideration. Indeed, these functions largely constitute a superset of the list recommended above, and are all clearly implementable. However, we felt it advantageous to require only the functions in the above list, in order to allow implementers the freedom to add value by providing additional functions.

A final possibility was the construction of an ad hoc list of Special Functions. We rejected this as the least defensible of the alternatives, since there are no obvious criteria for accepting some and rejecting others, other than a feel for general utility.

We note, in passing, that we have omitted both the cylindrical and the spherical Hankel functions from our list of candidates for standardization. Our motivation for this is straightforward: we have restricted ourselves to functions taking real-valued arguments, and producing real-valued results. The Hankel functions do not qualify.

C. Function Templates or Overloaded Functions?

It is possible to declare the desired additional <cmath> functionality in either of two ways: as (specialized) function templates or as families of overloaded functions. We recommend overloading.

This recommendation based in significant part on arguments favoring consistency of form with the existing contents of the affected headers: There are no templates today in <math.h> or in <cmath>. Further, we are unaware of any current Special Functions implementation that is based on template technology.

Additional reasons take into account the possible future extension of the new functions to additional headers (such as, for example, <complex>). To do so in the presence of function templates would raise such issues as the location of the primary template and the concomitant need to coordinate multiple cooperating headers. We prefer to avoid such entanglements.

D. Exceptions or Error Codes?

Most of the proposed functions must advise their caller of domain and/or range errors. In the context of C++, this would arguably be best done by throwing appropriate exceptions.

However, no standard functions in <cmath> today throw exceptions. Rather, mimicking the behavior of the functions in <math.h>, each sets a global errno variable to a suitable code (i.e., EDOM or ERANGE) defined in <cerrno>.

We recommend that this existing behavior be preserved with respect to the proposed new functions. Not only is this a matter of consistency, but it preserves the possibility that a compatible version of this proposal might be incorporated into a future revision of C99.

E. Traditional or Extended Error Codes?

Having recommended, in the previous subsection, the continued use of error codes, a further question arises: Are the existing EDOM, ERANGE error codes sufficient to the needs of the proposed Special Functions?

We note that several current implementations of some of the proposed Special Functions do define their own codes to supplement the codes mandated by the standard. They incorporate codes for such situations as overflow, underflow, loss of precision, singularities, and the like, and make these codes accessible via additional global variables analogous to errno, or via other means.

Nothing in this proposal should be construed as preventing implementors from such optional extensions. However, this proposal recommends against requiring such behavior. We base this recommendation on consistency with current standards. In particular, we call attention to §7.12.1 of C99 [ISO:9899].

F. Traditional or Descriptive Function Names?

In selecting names for the proposed new functions, we were moved to retain their traditional (mathematical) names. For example, we kept the names of a few functions which are customarily denoted using Greek letters (spelled out, of course: beta, zeta [but tgamma for compatibility with C99]). We also kept the function names erf and erfc because they are both traditional and descriptive, as well as for compatibility with C99.

In the remaining (majority of) cases, the traditional mathematical names are mostly single letters (e.g., J and N). We judged such names to be too brief and insufficiently descriptive for programming purposes. Also, such one-character names are frequently reserved by coding standards to signify local variables. For these reasons, we recommend against use of the traditional names. Instead, we chose such names as bessel_J and neumann_N, combining a descriptive prefix with a traditional suffix.

We note in passing that this policy resulted in a few pairs of our names that differ only in the case of a single character (e.g., bessel_J and bessel_j). In each instance, this is an artifact of the traditional mathematical naming convention that we preserved on the basis of prior art as the accepted canonical mathematical nomenclature. While some coding standards recommend or require avoidance of multiple identifiers that differ only in case, we believe it appropriate in the present context to embrace case-sensitivity in distinguishing otherwise-identical names.

G. Real- or Complex-valued Domains and Results?

Many of the proposed special functions have definitions over some or all of the complex plane as well as over some or all of the real numbers. Further, some of these functions can produce complex results, even over real-valued arguments. The present proposal restricts itself by considering only real-valued arguments and (correspondingly) real-valued results.

Our investigation of the alternative led us to realize that the complex landscape for the special functions is figuratively dotted with land mines. In coming to our recommendation, we gave weight to the statement from a respected colleague that "Several Ph.D. dissertations would [or could] result from efforts to implement this set of functions over the complex domain." This led us to take the position that there is insufficient prior art in this area to serve as a basis for standardization, and that such standardization would be therefore premature. While we could perhaps consider standardizing some subset of the special functions over the complex domain, we far prefer to treat this set of special functions as a unit.

We further ruled out (via domain errors, for example) the possibility that these special functions could return complex results. In making this recommendation, we followed the past practice of C++ and of C99: functions taking real arguments always return real results; only functions taking complex arguments return complex results. Perhaps the best example of this is the sqrt function: it always returns a value of type matching its parameter's type. To do otherwise opens the door to a number of small, but bothersome technical issues. As one example, which header (<cmath> or <complex>) would declare such a function whose domain and range are different?

Finally, none of our colleagues or reviewers has presented any compelling need or rationale for the extension to the complex domain or range. While there would certainly be some segments of the user community that could take advantage of such functionality (and we certainly don't mean to prohibit vendors from providing such as extensions), there seems to be insufficient demand to require such at present.

H. Which Conventions?

Among the functions in this proposal, there are several (e.g., spherical harmonics: sph_Y()) for which multiple phase and normalization conventions exist. Because the choices do not easily co-exist, we have opted to resolve any such ambiguities by appealing to a common source, the aforementioned standard Handbook of Mathematical Functions by [Abramowitz & Stegun].

I. Relationship to LIA?

It has been suggested that the functions constituting the subject of this proposal be reviewed within the framework of one or more parts of the International Standard for Language Independent Arithmetic [ISO:10967-1, ISO:10967-2, ISO:10967-3]. While we applaud the motivation underlying the suggestion, we recommend against such a perspective, believing it to be premature and not yet feasible.

Only LIA Part 1 has been formally adopted as an International Standard: LIA Part 2 is a Final Draft, while LIA Part 3 is only a Committee Draft. Further, LIA Part 1 does not speak to the subject of the present proposal, while Parts 2 and 3 restrict themselves to coverage of "elementary" numerical functions such as those already part of the C++ standard library. Thus, none of the LIA documents addresses (or even mentions) any of the functions comprising the present proposal.

Finally, we note that the C++ standards body has to date not adopted any "conformity statement" regarding the LIA "binding standard" to be used by a conforming C++ implementation. In the absence of such guidance, it is unclear how to apply to C++ the principles (let alone the details) of the LIA documents.

In our view, the present situation is best summarized as constituting a lack of prior art: C++ has not determined how LIA is to apply to an implementation and, in any event, none of the functions comprising the present proposal are in the scope of the LIA documents as currently drafted. For these reasons, we believe there is today an insufficient basis on which to evaluate the present proposal with respect to LIA.

IV. Proposed Text

Note: the wording presented in these subsections describes the contents of the <cmath> header. The changes to describe analogous contents in <math.h> are straightforward and consist primarily of function renaming to avoid overloading.

A. To be inserted into Table 80 (Clause 26)

bessel_I
bessel_J
bessel_K
bessel_j
beta
ei
ellint_E
ellint_E2
ellint_F
ellint_K
ellint_P
ellint_P2


hermite
hyperg_1F1
hyperg_2F1
laguerre_0
laguerre_m
legendre_Pl
legendre_Plm
neumann_N
neumann_n
sph_Y

zeta



B. New section to be added to Clause 26

26.x.1 Synopsis

  // (26.x.2) cylindrical Bessel functions (of the first kind):
  double       bessel_J( double l, double x );
  float        bessel_J( float l, float  x );
  long double  bessel_J( long double l, long double x );

  // (26.x.3) cylindrical Neumann functions;
  // cylindrical Bessel functions (of the second kind):
  double       neumann_N( double l, double x );
  float        neumann_N( float l, float  x );
  long double  neumann_N( long double l, long double x );

  // (26.x.4.1) regular modified cylindrical Bessel functions:
  double       bessel_I( double l, double x );
  float        bessel_I( float l, float  x );
  long double  bessel_I( long double l, long double x );

  // (26.x.4.2) irregular modified cylindrical Bessel functions:
  double       bessel_K( double l, double x );
  float        bessel_K( float l, float  x );
  long double  bessel_K( long double l, long double x );

  // (26.x.5) spherical Bessel functions (of the first kind):
  double       bessel_j( double l, double x );
  float        bessel_j( float l, float  x );
  long double  bessel_j( long double l, long double x );

  // (26.x.6) spherical Neumann functions;
  // spherical Bessel functions (of the second kind):
  double       neumann_n( double l, double x );
  float        neumann_n( float l, float  x );
  long double  neumann_n( long double l, long double x );

  // (26.x.7) Legendre polynomials:
  double       legendre_Pl( unsigned l, double x )
  float        legendre_Pl( unsigned l, float x )
  long double  legendre_Pl( unsigned l, long double x )

  // (26.x.8) associated Legendre functions:
  double       legendre_Plm( unsigned l, unsigned m, double x )
  float        legendre_Plm( unsigned l, unsigned m, float x )
  long double  legendre_Plm( unsigned l, unsigned m, long double x )

  // (26.x.9) spherical harmonics:
  double       sph_Y( unsigned l, int m, double theta, double phi )
  float        sph_Y( unsigned l, int m, float theta, float phi )
  long double  sph_Y( unsigned l, int m, long double theta, long double phi )

  // (26.x.10) Hermite polynomials:
  double       hermite( unsigned n, double x );
  float        hermite( unsigned n, float x );
  long double  hermite( unsigned n, long double x );

  // (26.x.11) Laguerre polynomials:
  double       laguerre_0(unsigned n, double x);
  float        laguerre_0(unsigned n, float x);
  long double  laguerre_0(unsigned n, long double x);

  // (26.x.12) associated Laguerre polynomials:
  double       laguerre_m(unsigned n, unsigned m, double x);
  float        laguerre_m(unsigned n, unsigned m, float x);
  long double  laguerre_m(unsigned n, unsigned m, long double x);

  // (26.x.13) hypergeometric functions:
  double       hyperg_2F1(double a, double b, double c, double x) ;
  float        hyperg_2F1(float a, float b, float c, float x) ;
  long double  hyperg_2F1(long double a, long double b, long double c, long double x) ;

  // (26.x.14) confluent hypergeometric functions:
  double       hyperg_1F1(double a, double c, double x) ;
  float        hyperg_1F1(float a, float c, float x) ;
  long double  hyperg_1F1(long double a, long double c, long double x) ;

  // (26.x.15.1) (incomplete) elliptic integral of the first kind:
  double       ellint_F( double k, double phi );
  float        ellint_F( float k, float phi );
  long double  ellint_F( long double k, long double phi );

  // (26.x.15.2) (complete) elliptic integral of the first kind:
  double       ellint_K( double k );
  float        ellint_K( float k );
  long double  ellint_K( long double k );

  // (26.x.16.1) (incomplete) elliptic integral of the second kind:
  double       ellint_E( double k, double phi );
  float        ellint_E( float k, float phi );
  long double  ellint_E( long double k, long double phi );

  // (26.x.16.2) (complete) elliptic integral of the second kind:
  double       ellint_E2( double k );
  float        ellint_E2( float k );
  long double  ellint_E2( long double k );

  // (26.x.17.1) (incomplete) elliptic integral of the third kind:
  double       ellint_P( double k, double n, double phi );
  float        ellint_P( float k, float n, float phi );
  long double  ellint_P( long double k, long double n, long double phi );

  // (26.x.17.2) (complete) elliptic integral of the third kind:
  double       ellint_P2( double k, double n, double phi );
  float        ellint_P2( float k, float n, float phi );
  long double  ellint_P2( long double k, long double n, long double phi );

  // (26.x.19) beta function:
  double       beta( double x, double y );
  float        beta( float x, float y );
  long double  beta( long double x, long double y );

  // (26.x.20) exponential integral:
  double       ei( double );
  float        ei( float );
  long double  ei( long double );

  // (26.x.22) Riemann zeta function:
  double       zeta( double );
  float        zeta( float );
  long double  zeta( long double );

26.x.2 cylindrical Bessel functions (of the first kind)

  double       bessel_J( double l, double x );
  float        bessel_J( float l, float  x );
  long double  bessel_J( long double l, long double x );

Effects:
These functions compute the cylindrical Bessel functions of the first kind (as defined by [Abramowitz & Stegun, §9.1.10, etc.]) of their respective arguments l and x. A domain error occurs if l is less than zero. A range error (due to underflow) occurs if x is too close to any of the roots of the bessel_J function.

Returns:
The bessel_J functions return Jl(x) as denoted in [ISO:31, Item No. 11-14.1].

26.x.3 cylindrical Neumann functions; cylindrical Bessel functions (of the second kind)

  double       neumann_N( double l, double x );
  float        neumann_N( float l, float  x );
  long double  neumann_N( long double l, long double x );

Effects:
These functions compute the cylindrical Neumann functions (as defined by [Abramowitz & Stegun, §9.1.2, etc.]), also known as the cylindrical Bessel functions of the second kind, of their respective arguments l and x. A domain error occurs if x is less than or equal to zero. A range error (due to overflow) occurs (a) if the magnitude of l is too large, or (b) if x is too small. A range error (due to underflow) occurs if x is too close to any of the roots of the neumann_N function.

Returns:
The neumann_N functions return Nl(x) as denoted in [ISO:31, Item No. 11-14.2].

26.x.4.1 regular modified cylindrical Bessel functions

  double       bessel_I( double l, double x );
  float        bessel_I( float l, float  x );
  long double  bessel_I( long double l, long double x );

Effects:
These functions compute the regular modified cylindrical Bessel functions (as defined by [Abramowitz & Stegun, §9.6.3, etc.) of their respective arguments l and x. A domain error occurs if l is not an integer and either (a) l is less than or equal to zero or (b) x is less than or equal to zero. A range error occurs if x is too large.

Returns:
The bessell_I functions return Il(x) as denoted in [ISO:31, Item No. 11-14.4].

26.x.4.2 irregular modified cylindrical Bessel functions

  double       bessel_K( double l, double x );
  float        bessel_K( float l, float  x );
  long double  bessel_K( long double l, long double x );

Effects:
These functions compute the irregular modified cylindrical Bessel functions (as defined by [Abramowitz & Stegun, §9.6.4, etc.]) of their respective arguments l and x. A domain error occurs if l is not an integer and either (a) l is less than or equal to zero or (b) x is less than or equal to zero. A range error occurs (a) if the magnitude of l is too large, or (b) if x is too small.

Returns:
The bessell_K functions return Kl(x) as denoted in [ISO:31, Item No. 11-14.4].

26.x.5 spherical Bessel functions (of the first kind)

  double       bessel_j( double l, double x );
  float        bessel_j( float l, float  x );
  long double  bessel_j( long double l, long double x );

Effects:
These functions compute the spherical Bessel functions of the first kind (as defined by [Abramowitz & Stegun, §10.1.1, etc.) of their respective arguments l and x. A domain error occurs if l is less than zero. A range error (due to overflow) occurs if l equals zero and the magnitude of x is too small. A range error (due to underflow) occurs if x is too close to any of the roots of the bessel_j function.

Returns:
The bessell_j functions return jl(x) as denoted in [ISO:31, Item No. 11-14.5].

26.x.6 spherical Neumann functions; spherical Bessel functions (of the second kind)

  double       neumann_n( double l, double x );
  float        neumann_n( float l, float  x );
  long double  neumann_n( long double l, long double x );

Effects:
These functions compute the spherical Neumann functions (as defined by [Abramowitz & Stegun, §10.1.1, etc.]), also known as the spherical Bessel functions of the second kind, of their respective arguments x. A domain error occurs if x is less than or equal to zero. A range error (due to overflow) occurs (a) if the magnitude of l is too large, or (b) if the magnitude of x is too small. A range error (due to underflow) occurs if x is too close to any of the roots of the neumann_n function.

Returns:
The neumann_n functions return nl(x) as denoted in [ISO:31, Item No. 11-14.6].

26.x.7 Legendre polynomials

  double       legendre_Pl( unsigned l, double x )
  float        legendre_Pl( unsigned l, float x )
  long double  legendre_Pl( unsigned l, long double x )

Effects:
These functions compute the Legendre polynomials (as defined by [Abramowitz & Stegun, §8.1.1, etc.]) of their respective arguments l and x. A domain error occurs if x is greater than one. A range error (due to underflow) occurs if x is too close to any of the roots of the legendre_Pl function.

Returns:
The legendre_Pl functions return Pl(x) as denoted in [ISO:31, Item No. 11-14.8].

26.x.8 associated Legendre functions

  double       legendre_Plm( unsigned l, unsigned m, double x )
  float        legendre_Plm( unsigned l, unsigned m, float x )
  long double  legendre_Plm( unsigned l, unsigned m, long double x )

Effects:
These functions compute the associated Legendre functions (as defined by [Abramowitz & Stegun, §8.1.1, etc.]) of their respective arguments l, m, and x. A domain error occurs (a) if m is greater than l, or (b) if x is greater than one. A range error (due to overflow) occurs if l is too large. A range error (due to underflow) occurs if x is too close to any of the roots of the legendre_Plm function.

Returns:
The legendre_Plm functions return Plm(x) as denoted in [ISO:31, Item No. 11-14.9].

26.x.9 spherical harmonics

  double       sph_Y( unsigned l, int m, double theta, double phi )
  float        sph_Y( unsigned l, int m, float theta, float phi )
  long double  sph_Y( unsigned l, int m, long double theta, long double phi )

Effects:
These functions compute spherical harmonic functions (as defined by [Abramowitz & Stegun, §8.1.1 footnote 2, etc.) of their respective arguments l, m, theta, and phi. A domain error occurs if the magnitude of m is greater than l. A range error (due to overflow) occurs if l is too large. A range error (due to underflow) occurs if either theta or phi is too close to any of the roots of the sph_Y function.

Returns:
The sph_Y functions return Ylm(φ,θ) as denoted in [ISO:31, Item No. 11-14.10].

26.x.10 Hermite polynomials

  double       hermite( unsigned n, double x );
  float        hermite( unsigned n, float x );
  long double  hermite( unsigned n, long double x );

Effects:
These functions compute the Hermite polynomials (as defined by [Abramowitz & Stegun, §13.6.17 and -.18, etc.) of their respective arguments n and x. A range error (due to overflow) occurs if the magnitude of x is too large. A range error (due to underflow) occurs if x is too close to any of the roots of the hermite function.

Returns:
The hermite functions return Hn(x) as denoted in [ISO:31, Item No. 11-14.11].

26.x.11 Laguerre polynomials

  double       laguerre_0(unsigned n, double x);
  float        laguerre_0(unsigned n, float x);
  long double  laguerre_0(unsigned n, long double x);

Effects:
These functions compute the Laguerre polynomials (as defined by [Abramowitz & Stegun, §13.6.9, etc.]) of their respective arguments n and x. A range error occurs (due to overflow) if the magnitude of x is too large. A range error (due to underflow) occurs if x is too close to any of the roots of the laguerre_0 function.

Returns:
The laguerre_0 functions return Ln(x) as denoted in [ISO:31, Item No. 11-14.12].

26.x.12 associated Laguerre polynomials

  double       laguerre_m(unsigned n, unsigned m, double x);
  float        laguerre_m(unsigned n, unsigned m, float x);
  long double  laguerre_m(unsigned n, unsigned m, long double x);

Effects:
These functions compute the associated Laguerre polynomials (as defined by [Abramowitz & Stegun, §13.6.9, etc.]) of their respective arguments n, m, and x. A domain error occurs if m is greater than n. A range error (due to overflow) occurs if the magnitude of x is too large. A range error (due to underflow) occurs if x is too close to any of the roots of the laguerre_m function.

Returns:
The laguerre_m functions return Lnm(x) as denoted in [ISO:31, Item No. 11-14.13].

26.x.13 hypergeometric functions

  double       hyperg_2F1(double a, double b, double c, double x) ;
  float        hyperg_2F1(float a, float b, float c, float x) ;
  long double  hyperg_2F1(long double a, long double b, long double c, long double x) ;

Effects:
These functions compute the hypergeometric functions (as defined by [Abramowitz & Stegun, §15.1.1, etc.]) of their respective arguments a, b, c, and x. A domain error occurs if the magnitude of x is greater than or equal to one. A range error (due to overflow) occurs if the magnitude of x is too close to one at the same time that c-a-b is an integer. A range error (due to underflow) occurs if x is too close to any of the roots of the hyperg_2F1 function.

Returns:
The hyperg_2F1 functions return F(a,b;c;x) as denoted in [ISO:31, Item No. 11-14.14].

26.x.14 confluent hypergeometric functions

  double       hyperg_1F1(double a, double c, double x) ;
  float        hyperg_1F1(float a, float c, float x) ;
  long double  hyperg_1F1(long double a, long double c, long double x) ;

Effects:
These functions compute the confluent hypergeometric functions (as defined by [Abramowitz & Stegun, §13.1.2, etc.]) of their respective arguments a, c, and x. A domain error occurs (a) if c is a negative integer, or (b) if c is zero. A range error (due to overflow) occurs if the magnitude of x is too large. A range error (due to underflow) occurs if x is too close to any of the roots of the hyperg_1F1 function.

Returns:
The hyperg_1F1 functions return F(a;c;x) as denoted in [ISO:31, Item No. 11-14.15].

26.x.15.1 (incomplete) elliptic integral of the first kind

  double       ellint_F( double k, double phi );
  float        ellint_F( float k, float phi );
  long double  ellint_F( long double k, long double phi );

Effects:
These functions compute the incomplete elliptic integral of the first kind (as defined by [Abramowitz & Stegun, §17.2.6, etc.) of their respective arguments k and phi. A domain error occurs (a) if k is less than or equal to zero, or (b) if k is greater than or equal to one, or (c) if phi is negative. A range error occurs if the magnitude of k is too close to one.

Returns:
The ellint_F functions return F(k,φ) as denoted in [ISO:31, Item No. 11-14.16].

26.x.15.2 (complete) elliptic integral of the first kind

  double       ellint_K( double k );
  float        ellint_K( float k);
  long double  ellint_K( long double k);

Effects:
These functions compute the complete elliptic integral of the first kind (as defined by [Abramowitz & Stegun, §17.3.1, etc.]) of their respective arguments k. A domain error occurs (a) if k is less than or equal to zero, or (b) if k is greater than or equal to one. A range error occurs if the magnitude of k is too close to one.

Returns:
The ellint_K functions return K(k) as denoted in [ISO:31, Item No. 11-14.16].

26.x.16.1 (incomplete) elliptic integral of the second kind

  double       ellint_E( double k, double phi );
  float        ellint_E( float k, float phi );
  long double  ellint_E( long double k, long double phi );

Effects:
These functions compute the incomplete elliptic integral of the second kind (as defined by [Abramowitz & Stegun, §17.2.9, etc.]) of their respective arguments k and phi. A domain error occurs (a) if k is less than or equal to zero, or (b) if k is greater than or equal to one, or (c) if phi is negative. A range error occurs if the magnitude of k is too close to one.

Returns:
The ellint_E functions return E(k,φ) as denoted in [ISO:31, Item No. 11-14.17].

26.x.16.2 (complete) elliptic integral of the second kind

  double       ellint_E2( double k );
  float        ellint_E2( float k );
  long double  ellint_E2( long double k );

Effects:
These functions compute the complete elliptic integral of the second kind (as defined by [Abramowitz & Stegun, §17.2.9, etc.) of their respective arguments k and phi. A domain error occurs (a) if k is less than or equal to zero, or (b) if k is greater than or equal to one. A range error occurs if the magnitude of k is too close to one.

Returns:
The ellint_E2 functions return E(k) as denoted in [ISO:31, Item No. 11-14.17].

26.x.17.1 (incomplete) elliptic integral of the third kind

  double       ellint_P( double k, double n, double phi );
  float        ellint_P( float k, float n, float phi );
  long double  ellint_P( long double k, long double n, long double phi );

Effects:
These functions compute the incomplete elliptic integral of the third kind (as defined by [Abramowitz & Stegun, §17.7.1, etc.) of their respective arguments k, n, and phi. A domain error occurs (a) if k is less than or equal to zero, or (b) if k is greater than or equal to one, or (c) if phi is negative. A range error occurs if the magnitude of k is too close to one.

Returns:
The ellint_P functions return Π(k,n,φ) as denoted in [ISO:31, Item No. 11-14.18].

26.x.17.2 (complete) elliptic integral of the third kind

  double       ellint_P2( double k, double n );
  float        ellint_P2( float k, float n );
  long double  ellint_P2( long double k, long double n );

Effects:
These functions compute the complete elliptic integral of the third kind (as defined by [Abramowitz & Stegun, §17.7.2, etc.) of their respective arguments k, n, and phi. A domain error occurs (a) if k is less than or equal to zero, or (b) if k is greater than or equal to one. A range error occurs if the magnitude of k is too close to one.

Returns:
The ellint_P2 functions return Π(k,n,π/2) as denoted in [ISO:31, Item No. 11-14.18].

26.x.19 beta function

  double       beta( double x, double y );
  float        beta( float x, float y );
  long double  beta( long double x, long double y );

Effects:
These functions compute the beta function (as defined by [Abramowitz & Stegun, §6.2 and §6.1]) of their respective arguments x and y. A domain error occurs (a) if either x or y is a negative integer, or (b) if either x or y is zero. A range error occurs if the magnitude of x or the magnitude of y is too large or too small.

Returns:
The beta functions return Β(x,y) as denoted in [ISO:31, Item No. 11-14.20].

26.x.20 exponential integral

  double       ei( double x );
  float        ei( float x );
  long double  ei( long double x );

Effects:
These functions compute the exponential integral (as defined by [Abramowitz & Stegun, §5.1.1, etc.]) of their respective arguments x. A domain error occurs if x is less than or equal to zero. A range error (due to overflow) occurs (a) if x is to small, or (b) if x x is too large. A range error (due to underflow) occurs if x is too close to the single root of the ei function.

Returns:
The ei functions return Ei(x) as denoted in [ISO:31, Item No. 11-14.21].

26.x.22 Riemann zeta function

  double       zeta( double x );
  float        zeta( float x );
  long double  zeta( long double x );

Effects:
These functions compute the Riemann zeta function (as defined by [Abramowitz & Stegun, §23.2.1, etc.]) of their respective arguments x. A range error (due to overflow) occurs if x is too close to one.

Returns:
The zeta functions return ζ(x) as denoted in [ISO:31, Item No. 11-14.23].

V. Relationship to Earlier Proposal

We noted in section II that there is a small potential overlap between this proposal (which is based on [ISO:31]) and an earlier proposal (which is based on [C99]). However, no competition is intended between the two proposals. Indeed, we have coordinated efforts to ensure that no incompatibilities result.

In particular, this section describes three functions (erf, erfc, and tgamma) originating in the C99 library and previously proposed for C++ standardization as part of [Plauger]. They are mentioned in the present proposal for two reasons, however: (1) these functions are traditionally characterized as mathematical special functions, and (2) they form part of [ISO:31], on which our proposal is based.

We emphasize that these functions are here described solely in the interest of completeness so that the full spectrum of envisioned mathematical special functions may be viewed together. For purposes of the C++ standardization effort, these functions should be considered part of the [Plauger] proposal. The language in the remainder of this section should therefore be considered solely for informative purposes.

To be inserted into Table 80 (Clause 26)

erf
erfc
tgamma

26.x.1 Synopsis

  // (26.x.18) gamma function:
  double       tgamma( double );
  float        tgamma( float );
  long double  tgamma( long double );

  // (26.x.21.1) error function:
  double       erf( double );
  float        erf( float );
  long double  erf( long double );

  // (26.x.21.2) complementary error function:
  double       erfc( double );
  float        erfc( float );
  long double  erfc( long double );

26.x.18 gamma function

  double       tgamma( double x );
  float        tgamma( float x );
  long double  tgamma( long double x );

Effects:
These functions compute the gamma function (as defined by [Abramowitz & Stegun, §6.1.1, etc.]) of their respective arguments x. A domain error occurs (a) if x is a negative integer, or (b) if x is zero. A range error occurs if the magnitude of x is too large or too small.

Returns:
The tgamma functions return Γ(x) as denoted in [ISO:31, Item No. 11-14.19].

Note:
This family of overloaded functions corresponds to the functions specified in [C99, §7.12.8.4].

26.x.21.1 error function

  double       erf( double x );
  float        erf( float x );
  long double  erf( long double x );

Effects:
These functions compute the error function (as defined by [Abramowitz & Stegun, §7.1.1, etc.]) of their respective arguments x.

Returns:
The erf functions return erf(x) as denoted in [ISO:31, Item No. 11-14.22].

Note:
This family of overloaded functions corresponds to the functions specified in [C99, §7.12.8.1].

26.x.21.2 complementary error function

  double       erfc( double x );
  float        erfc( float x );
  long double  erfc( long double x );

Effects:
These functions compute the complementary error function (as defined by [Abramowitz & Stegun, §7.1.2, etc.]) of their respective arguments x.

Returns:
The erf functions return erfc(x) as denoted in [ISO:31, Item No. 11-14.22].

Note:
This family of overloaded functions corresponds to the functions specified in [C99, §7.12.8.2].

VI. Acknowledgements

It is a pleasure to acknowledge the significant contributions provided by a number of colleagues at Fermilab: James Amundson, Mark Fischler, Jeffrey Kallenbach, Leo Michelotti, and Marc Paterno. Their active participation has materially improved this proposal, and their inspiration and support are deeply appreciated.

We extend special thanks to our Fermilab colleague John Marraffino for his extensive involvement with all aspects of this proposal's development.

We have also received valuable input and feedback from Matt Austern, Beman Dawes, Gabriel Dos Reis, and Fred J. Tydeman. We are grateful to them and to our outside reviewers for their careful consideration of earlier drafts of this document: Their thoughtful comments and suggestions inspired several refinements and enhancements to this proposal.

Finally, we appreciate the support of the Fermi National Accelerator Laboratory's Computing Division, sponsors of our participation in the C++ standards effort. Thank you, one and all.

VII. References

[Abramowitz & Stegun]
Abramowitz, Milton, and Irene A. Stegun (eds.): Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables, volume 55 of National Bureau of Standards Applied Mathematics Series. U. S. Government Printing Office, Washington, DC: 1964. Reprinted with corrections, Dover Publications: 1972. http://members.fortunecity.com/aands.
[Austern]
Austern, Matt: Notes on Standard Library Extensions. WG21/N1314 (same as J16/01-0028): 17 May 2001.
[Dawes]
Dawes, Beman: Library Technical Report Proposals and Issues List (Revision 4). WG21/N1397 (same as J16/02-0055): 2002.
[GNU C]
GNU C Library Steering Committee The GNU C Library. http://www.gnu.org/manual/glibc-2.0.6/html_node/libc_toc.html.
[GSL]
Galassi, Mark, et al.: The GNU Scientific Library. 2002. http://www.gnu.org/software/gsl/gsl.html.
[IMSL]
Visual Numerics, Inc.: IMSL C Numerical Library Version 5.0. Unpublished. http://www.vni.com/products/imsl/docs/cmath.pdf.
[ISO:31]
International Standards Organization: Quantities and units, Third edition. International Standard ISO 31-11:1992. ISBN 92-67-10185-4.
[ISO:9899]
International Standards Organization: Programming Languages -- C, Second edition. International Standard ISO/IEC 9899:1999.
[ISO:10967-1]
International Standards Organization: Information technology -- Language independent arithmetic -- Part 1: Integer and floating point arithmetic. International Standard ISO/IEC 10967-1:1994.
[ISO:10967-2]
International Standards Organization: Information technology -- Language independent arithmetic -- Part 2: Elementary numerical functions. Draft International Standard ISO/IEC FDIS 10967-2:2000.
[ISO:10967-3]
International Standards Organization: Information technology -- Language independent arithmetic -- Part 3: Complex integer and floating point arithmetic and complex elementary numerical functions. Draft International Standard ISO/IEC CD 10967-3.1:2002.
[ISO:14882]
International Standards Organization: Programming Languages -- C++. International Standard ISO/IEC 14882:1998.
[Josuttis]
Josuttis, Nicolai: A New Work Item Proposal: Technical Report for Library Issues. WG21/N1283 (same as J16/00-0060): 26 October 2000.
[NAG]
The Numerical Algorithms Group: The NAG C Library Manual, Mark 7. The Numerical Algorithms Group, Ltd., Oxford, UK: 2002. http://www.nag.com/numeric/CL/manual/html/CLlibrarymanual.asp.
[Plauger]
Plauger, P. J.: Proposed C99 Library Additions to C++ (Revised). WG21/N1372 (same as J16/02-0030): 2002. http://anubis.dkuug.dk/jtc1/sc22/wg21/docs/papers/2002/n1372.txt.
[SGI]
Silicon Graphics, Inc.: Introduction to mathematical library functions. 2002. http://techpubs.sgi.com/library/tpl/cgi-bin/getdoc.cgi?coll=0650&db=man&fname=/usr/share/catman/p_man/cat3/standard/math.z&srch=math.
[SLATEC]
Fong, Kirby W., et al.: Guide to the SLATEC Common Mathematical Library. July 1993. http://www.netlib.org/slatec/guide.