Most of the functions for functions of a single real argument have been based on truncated Chebyshev expansions. This method of approximation was adopted as a compromise between the conflicting requirements of efficiency and ease of implementation on many different machine ranges. For details of the reasons behind this choice and the production and testing procedures followed in constructing this chapter see Schonfelder (1976).

Basically, if the function to be approximated is , then for  an approximation of the form

is used ( denotes, according to the usual convention, a summation in which the first term is halved), where  is some suitable auxiliary function which extracts any singularities, asymptotes and, if possible, zeros of the function in the range in question and  is a mapping of the general range  to the specific range [] required by the Chebyshev polynomials, . For a detailed description of the properties of the Chebyshev polynomials see Clenshaw (1962) and Fox and Parker (1968).

The essential property of these polynomials for the purposes of function approximation is that  oscillates between  and it takes its extreme values  times in the interval []. Therefore, provided the coefficients  decrease in magnitude sufficiently rapidly the error made by truncating the Chebyshev expansion after  terms is approximately given by

That is, the error oscillates between  and takes its extreme value  times in the interval in question. Now this is just the condition that the approximation be a mini-max representation, one which minimizes the maximum error. By suitable choice of the interval, [], the auxiliary function, , and the mapping of the independent variable, , it is almost always possible to obtain a Chebyshev expansion with rapid convergence and hence truncations that provide near mini-max polynomial approximations to the required function. The difference between the true mini-max polynomial and the truncated Chebyshev expansion is seldom sufficiently great enough to be of significance.

The evaluation of the Chebyshev expansions follows one of two methods. The first and most efficient, and hence the most commonly used, works with the equivalent simple polynomial. The second method, which is used on the few occasions when the first method proves to be unstable, is based directly on the truncated Chebyshev-series, and uses backward recursion to evaluate the sum. For the first method, a suitably truncated Chebyshev expansion (truncation is chosen so that the error is less than the machine precision) is converted to the equivalent simple polynomial. That is, we evaluate the set of coefficients  such that

The polynomial can then be evaluated by the efficient Horner's method of nested multiplications,

This method of evaluation results in efficient functions but for some expansions there is considerable loss of accuracy due to cancellation effects. In these cases the second method is used. It is well known that if

then

and this is always stable. This method is most efficiently implemented by using three variables cyclically and explicitly constructing the recursion.

That is,

The auxiliary functions used are normally functions compounded of simple polynomial (usually linear) factors extracting zeros, and the primary compiler-provided functions, sin, cos, ln, exp, sqrt, which extract singularities and/or asymptotes or in some cases basic oscillatory behaviour, leaving a smooth well-behaved function to be approximated by the Chebyshev expansion which can therefore be rapidly convergent.

The mappings of [] to [] used range from simple linear mappings to the case when  is infinite, and considerable improvement in convergence can be obtained by use of a bilinear form of mapping. Another common form of mapping is used when the function is even; that is, it involves only even powers in its expansion. In this case an approximation over the whole interval [] can be provided using a mapping . This embodies the evenness property but the expansion in  involves all powers and hence removes the necessity of working with an expansion with half its coefficients zero.

For many of the functions an analysis of the error in principle is given, namely, if  and  are the absolute errors in function and argument and  and  are the corresponding relative errors, then

If we ignore errors that arise in the argument of the function by propagation of data errors, etc., and consider only those errors that result from the fact that a real number is being represented in the computer in floating-point form with finite precision, then  is bounded and this bound is independent of the magnitude of . For example, on an 11-digit machine

(This of course implies that the absolute error  is also bounded but the bound is now dependent on .) However, because of this the last two relations above are probably of more interest. If possible the relative error propagation is discussed; that is, the behaviour of the error amplification factor  is described, but in some cases, such as near zeros of the function which cannot be extracted explicitly, absolute error in the result is the quantity of significance and here the factor  is described. In general, testing of the functions has shown that their error behaviour follows fairly well these theoretical error behaviours. In regions where the error amplification factors are less than or of the order of one, the errors are slightly larger than the above predictions. The errors are here limited largely by the finite precision of arithmetic in the machine, but  is normally no more than a few times greater than the bound on . In regions where the amplification factors are large, of order ten or greater, the theoretical analysis gives a good measure of the accuracy obtainable.

It should be noted that the definitions and notations used for the functions in this chapter are all taken from Abramowitz and Stegun (1972). You are strongly recommended to consult this book for details before using the functions in this chapter.

Four functions provided here are symmetrised variants of the classical (Legendre) elliptic integrals. These alternative definitions have been suggested by Carlson (1965), Carlson (1977a) and Carlson (1977b) and he also developed the basic algorithms used in this chapter.

The standard Carlson integral of the first kind is represented by

where  and at most one may be equal to zero.

The normalization factor, , is chosen so as to make

If any two of the variables are equal,  degenerates into the second function

where the argument restrictions are now  and .

This function is related to the logarithm or inverse hyperbolic functions if , and to the inverse circular functions if .

The Carlson integral of the second kind is defined by

with ,  and , but only one of  or  may be zero.

The function is a degenerate special case of the Carlson integral of the third kind

with  and  with at most one equality holding. Thus . The normalization of both these functions is chosen so that

The algorithms used for all these functions are based on duplication theorems. These allow a recursion system to be established which constructs a new set of arguments from the old using a combination of arithmetic and geometric means. The value of the function at the original arguments can then be simply related to the value at the new arguments. These recursive reductions are used until the arguments differ from the mean by an amount small enough for a Taylor series about the mean to give sufficient accuracy when retaining terms of order less than six. Each step of the recurrences reduces the difference from the mean by a factor of four, and as the truncation error is of order six, the truncation error goes like , where  is the number of iterations.

The above forms can be related to the more traditional canonical forms (see Section 17.2 in Abramowitz and Stegun (1972)), as follows.

If we write , where , we have

the classical elliptic integral of the first kind:

the classical elliptic integral of the second kind:

the classical elliptic integral of the third kind:

Also the classical complete elliptic integral of the first kind:

the complete elliptic integral of the second kind:

For convenience, Chapter s contains functions to evaluate classical and Carlson elliptic integrals.

The functions for Bessel and Airy functions of a real argument are based on Chebyshev expansions, as described in Section [Functions of a Single Real Argument]. The functions for functions of a complex argument, however, use different methods. These functions relate all functions to the modified Bessel functions  and  computed in the right-half complex plane, including their analytic continuations.  and  are computed by different methods according to the values of  and . The methods include power series, asymptotic expansions and Wronskian evaluations. The relations between functions are based on well known formulae (see Abramowitz and Stegun (1972)).

IMPORTANT ADVICE: users who encounter elliptic integrals in the course of their work are strongly recommended to look at transforming their analysis directly to one of the Carlson forms, rather than to the traditional canonical Legendre forms. In general, the extra symmetry of the Carlson forms is likely to simplify the analysis, and these symmetric forms are much more stable to calculate.

The function s21bac (nag_elliptic_integral_rc) for  is largely included as an auxiliary to the other functions for elliptic integrals. This integral essentially calculates elementary functions, e.g.,

In general this method of calculating these elementary functions is not recommended as there are usually much more efficient specific functions available in the Library. However, s21bac (nag_elliptic_integral_rc) may be used, for example, to compute  when  is close to 1, without the loss of significant figures that occurs when  and  are computed separately.

For computing the Bessel functions , ,  and  where  is real and , special functions are provided, which are much faster than the more general functions that allow a complex argument and arbitrary real . Similarly, special functions are provided for computing the Airy functions and their derivatives , , ,  for a real argument which are much faster than the functions for complex arguments.