\begin{page}{manpageXXs}{NAG On-line Documentation: s} \beginscroll \begin{verbatim} S(3NAG) Foundation Library (12/10/92) S(3NAG) S -- Approximations of Special Functions Introduction -- S Chapter S Approximations of Special Functions 1. Scope of the Chapter This chapter is concerned with the provision of some commonly occurring physical and mathematical functions. 2. Background to the Problems The majority of the routines in this chapter approximate real- valued functions of a single real argument, and the techniques involved are described in Section 2.1. In addition the chapter contains routines for elliptic integrals (see Section 2.2), Bessel and Airy functions of a complex argument (see Section 2.3) , and exponential of a complex argument. 2.1. Functions of a Single Real Argument Most of the routines 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 [7]. Basically if the function to be approximated is f(x), then for xis in [a,b] an approximation of the form --' f(x)=g(x) > C T (t) -- r r r=0 --' is used, ( > denotes, according to the usual convention, a -- summation in which the first term is halved), where g(x) is some suitable auxiliary function which extracts any singularities, asymptotes and, if possible, zeros of the function in the range in question and t=t(x) is a mapping of the general range [a,b] to the specific range [-1,+1] required by the Chebyshev polynomials, T (t). For a detailed description of the properties of the r Chebyshev polynomials see Clenshaw [5] and Fox and Parker [6]. The essential property of these polynomials for the purposes of function approximation is that T (t) oscillates between +-1 and n it takes its extreme values n+1 times in the interval [-1,+1]. Therefore, provided the coefficients C decrease in magnitude r sufficiently rapidly the error made by truncating the Chebyshev expansion after n terms is approximately given by E(t)~=C T (t) n n That is the error oscillates between +-C and takes its extreme n value n+1 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, [a,b], the auxiliary function, g(x), and the mapping of the independent variable, t(x), 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 to be of significance. The evaluation of the Chebyshev expansions follows one of two methods. The first and most efficient, and hence 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 b such that r n-1 n-1 -- r --' y(t)= > b t = > C T (t). -- r -- r r r=0 r=0 The polynomial can then be evaluated by the efficient Horner's method of nested multiplications, y(t)=(b +t(b +t(b +...t(b +tb )))...). 0 1 2 n-2 n-1 This method of evaluation results in efficient routines 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 b =C n-1 n-1 b =2tb +C n-2 n-1 n-2 b =2tb -b +C , j=n-3,n-4,...,0 j j+1 j+2 j then --' 1 > C T (t)= -(b -b ) -- r r 2 0 2 r=0 and this is always stable. This method is most efficiently implemented by using three variables cyclically and explicitly constructing the recursion. That is, (alpha) = C n-1 (beta) = 2t(alpha)+C n-2 (gamma) = 2t(beta)-(alpha)+C n-3 (alpha) = 2t(gamma)-(beta)+C n-4 (beta) = 2t(alpha)-(gamma)+C n-5 ... ... (alpha) = 2t(gamma)-(beta)+C (say) 2 (beta) = 2t(alpha)-(gamma)+C 1 1 y(t) = t(beta)-(alpha)+ -C 2 0 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 [a,b] to [-1,+1] used, range from simple linear mappings to the case when b 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 (x)2 [-a,a] can be provided using a mapping t=2(-) -1. This embodies (a) the evenness property but the expansion in t involves all powers and hence removes the necessity of working with an expansion with half its coefficients zero. For many of the routines an analysis of the error in principle is given, viz, if E and (nabla) are the absolute errors in function and argument and (epsilon) and (delta) are the corresponding relative errors, then E~=|f'(x)|(nabla) E~=|xf'(x)|(delta) | xf'(x)| (epsilon)~=| ------|(delta) | f(x) | 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 (delta) is bounded and this bound is independent of the magnitude of x; e.g. on an 11-digit machine -11 |(delta)|<=10 . (This of course implies that the absolute error (nabla)=x(delta) is also bounded but the bound is now dependent on x). 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 |xf'(x)/f(x)| 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 |xf'(x)| 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 (epsilon) is normally no more than a few times greater than the bound on (delta). In regions where the amplification factors are large, order of 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 [1]. Users are strongly recommended to consult this book for details before using the routines in this chapter. 2.2. Approximations to Elliptic Integrals The functions provided here are symmetrised variants of the classic elliptic integrals. These alternative definitions have been suggested by Carlson (see [2], [3] and [4]) and he also developed the basic algorithms used in this chapter. The standard integral of the first kind is represented by infty 1 / dt R (x,y,z)= - | ----------------- F 2 / 0 \/(t+x)(t+y)(t+z) where x,y,z>=0 and at most one may be equal to zero. 1 The normalisation factor, -, is chosen so as to make 2 R (x,x,x)=1/\/x. F If any two of the variables are equal, R degenerates into the F second function infty 1 / dt R (x,y)=R (x,y,y)= - | ---------- C F 2 / 0 \/t+x(t+y) where the argument restrictions are now x>=0 and y/=0. This function is related to the logarithm or inverse hyperbolic functions if 00, x>=0 and y>=0 but only one of x or y may be zero. The function is a degenerate special case of the integral of the third kind infty 3 / dt R (x,y,z,(rho))= - | ---------------------------- J 2 / ( ) 0 (\/(t+x)(t+y)(t+z))(t+(rho)) with (rho)/=0, x,y,z>=0 with at most one equality holding. Thus R (x,y,z)=R (x,y,z,z). The normalisation of both these functions D J is chosen so that R (x,x,x)=R (x,x,x,x)=1/(x\/x) D J 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, -n the truncation error goes like (4096) , where n is the number of iterations. The above forms can be related to the more traditional canonical forms (see Abramowitz and Stegun [1], 17.2). 2 2 2 If we write q=cos (phi),r=1-m.sin (phi),s=1+n.sin (phi), where 1 0<(phi)<= -(pi), we have: the elliptic integral of the first 2 kind: sin(phi) -1/2 1-/2 / 2 2 F((phi)|m)= | (1-t ) (1-mt ) dt=sin(phi).R (q,r,1); / F 0 the elliptic integral of the second kind: sin(phi) -1/2 -1/2 / 2 2 E((phi)|m) = | (1-t ) (1-mt ) dt / 0 1 3 = sin(phi).R (q,r,1)- -m.sin (phi).R (q,r,1) F 3 D the elliptic integral of the third kind: sin(phi) -1/2 -1/2 / 2 2 2 -1 (Pi)(n;(phi)|m) = | (1-t ) (1-mt ) (1+nt ) dt / 0 1 3 =sin(phi).R (q,r,1)- -n.sin (phi).R (q,r,1,s) F 3 J Also the complete elliptic integral of the first kind: (pi)/2 / 2 -1/2 K(m)= | (1-m.sin (theta)) d(theta)=R (0,1-m,1); / F 0 the complete elliptic integral of the second kind: (pi)/2 / 2 1/2 1 E(m)= | (1-m.sin (theta)) d(theta)=R (0,1-m,1)- -mR (0,1-m,1). / F 3 D 0 2.3. Bessel and Airy Functions of a Complex Argument The routines for Bessel and Airy functions of a real argument are based on Chebyshev expansions, as described in Section 2.1. The routines for functions of a complex argument, however, use different methods. These routines relate all functions to the modified Bessel functions I (z) and K (z) computed in the (nu) (nu) right-half complex plane, including their analytic continuations. I and K are computed by different methods according to (nu) (nu) the values of zand (nu). The methods include power series, asymptotic expansions and Wronskian evaluations. The relations between functions are based on well known formulae (see Abramowitz and Stegun [1]). 2.4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. [2] Carlson B C (1977) Special Functions of Applied Mathematics. Academic Press. [3] Carlson B C (1965) On Computing Elliptic Integrals and Functions. J Math Phys. 44 36--51. [4] Carlson B C (1977) Elliptic Integrals of the First Kind. SIAM J Math Anal. 8 231--242. [5] Clenshaw C W (1962) Mathematical Tables. Chebyshev Series for Mathematical Functions. HMSO. [6] Fox L and Parker I B (1968) Chebyshev Polynomials in Numerical Analysis. Oxford University Press. [7] Schonfelder J L (1976 ) The Production of Special Function Routines for a Multi-Machine Library. Software Practice and Experience. 6(1) 3. Recommendations on Choice and Use of Routines 3.1. Elliptic Integrals 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 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 routine S21BAF for R is largely included as an auxiliary to C the other routines for elliptic integrals. This integral essentially calculates elementary functions, e.g. (( 1+x)2 ) lnx=(x-1).R (( ---) ,x),x>0; C(( 2 ) ) 2 arcsinx=x.R (1-x ,1),|x|<=1; C 2 arcsinhx=x.R (1+x ,1), etc C In general this method of calculating these elementary functions is not recommended as there are usually much more efficient specific routines available in the Library. However, S21BAF may be used, for example, to compute lnx/(x-1) when x is close to 1, without the loss of significant figures that occurs when lnx and x-1 are computed separately. 3.2. Bessel and Airy Functions For computing the Bessel functions J (x), Y (x), I (x) (nu) (nu) (nu) and K (x) where x is real and (nu)=0 or 1, special routines (nu) are provided, which are much faster than the more general routines that allow a complex argument and arbitrary real (nu)>=0 functions and their derivatives Ai(x), Bi(x), Ai'(x), Bi'(x) for a real argument which are much faster than the routines for complex arguments. 3.3. Index Airy function, Ai, real argument S17AGF Airy function, Ai', real argument S17AJF Airy function, Ai or Ai', complex argument, optionally S17DGF scaled Airy function, Bi, real argument S17AHF Airy function, Bi', real argument S17AKF Airy function, Bi or Bi', complex argument, optionally S17DHF scaled Bessel function, J , real argument S17AEF 0 Bessel function, J , real argument S17AFF 1 Bessel function, J , complex argument, optionally S17DEF (nu) scaled Bessel function, Y , real argument S17ACF 0 Bessel function, Y , real argument S17ADF 1 Bessel function, Y , complex argument, optionally S17DCF (nu) scaled Complement of the Error function S15ADF Cosine Integral S13ACF Elliptic integral, symmetrised, degenerate of 1st kind, S21BAF R C Elliptic integral, symmetrised, of 1st kind, R S21BBF F Elliptic integral, symmetrised, of 2nd kind, R S21BCF D Elliptic integral, symmetrised, of 3rd kind, R S21BDF J Erf, real argument S15AEF Erfc, real argument S15ADF Error function S15AEF Exponential, complex S01EAF Exponential Integral S13AAF Fresnel Integral, C S20ADF Fresnel Integral, S S20ACF Gamma function S14AAF Gamma function, incomplete S14BAF Generalized Factorial function S14AAF (1) (2) Hankel function H or H , complex argument, S17DLF (nu) (nu) optionally scaled Incomplete Gamma function S14BAF Jacobian elliptic functions, sn, cn, dn S21CAF Kelvin function, bei x S19ABF Kelvin function, ber x S19AAF Kelvin function, kei x S19ADF Kelvin function, ker x S19ACF Logarithm of Gamma function S14ABF Modified Bessel function, I , real argument S18AEF 0 Modified Bessel function, I , real argument S18AFF 1 Modified Bessel function, I , complex argument, S18DEF (nu) optionally scaled Modified Bessel function, K , real argument S18ACF 0 Modified Bessel function, K , real argument S18ADF 1 Modified Bessel function, K , complex argument, S18DCF (nu) optionally scaled Sine integral S13ADF S -- Approximations of Special Functions Contents -- S Chapter S Approximations of Special Functions z S01EAF Complex exponential, e S13AAF Exponential integral E (x) 1 S13ACF Cosine integral Ci(x) S13ADF Sine integral Si(x) S14AAF Gamma function S14ABF Log Gamma function S14BAF Incomplete gamma functions P(a,x) and Q(a,x) S15ADF Complement of error function erfc x S15AEF Error function erf x S17ACF Bessel function Y (x) 0 S17ADF Bessel function Y (x) 1 S17AEF Bessel function J (x) 0 S17AFF Bessel function J (x) 1 S17AGF Airy function Ai(x) S17AHF Airy function Bi(x) S17AJF Airy function Ai'(x) S17AKF Airy function Bi'(x) S17DCF Bessel functions Y (z), real a>=0, complex z, (nu)+a (nu)=0,1,2,... S17DEF Bessel functions J (z), real a>=0, complex z, (nu)+a (nu)=0,1,2,... S17DGF Airy functions Ai(z) and Ai'(z), complex z S17DHF Airy functions Bi(z) and Bi'(z), complex z (j) S17DLF Hankel functions H (z), j=1,2, real a>=0, complex z, (nu)+a (nu)=0,1,2,... S18ACF Modified Bessel function K (x) 0 S18ADF Modified Bessel function K (x) 1 S18AEF Modified Bessel function I (x) 0 S18AFF Modified Bessel function I (x) 1 S18DCF Modified Bessel functions K (z), real a>=0, complex (nu)+a z, (nu)=0,1,2,... S18DEF Modified Bessel functions I (z), real a>=0, complex (nu)+a z, (nu)=0,1,2,... S19AAF Kelvin function ber x S19ABF Kelvin function bei x S19ACF Kelvin function ker x S19ADF Kelvin function kei x S20ACF Fresnel integral S(x) S20ADF Fresnel integral C(x) S21BAF Degenerate symmetrised elliptic integral of 1st kind R (x,y) C S21BBF Symmetrised elliptic integral of 1st kind R (x,y,z) F S21BCF Symmetrised elliptic integral of 2nd kind R (x,y,z) D S21BDF Symmetrised elliptic integral of 3rd kind R (x,y,z,r) J \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs01eaf}{NAG On-line Documentation: s01eaf} \beginscroll \begin{verbatim} S01EAF(3NAG) Foundation Library (12/10/92) S01EAF(3NAG) S01 -- Approximations of Special Functions S01EAF S01EAF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose z S01EAF evaluates the exponential function e , for complex z. 2. Specification COMPLEX(KIND(1.0D0)) FUNCTION S01EAF (Z, IFAIL) INTEGER IFAIL COMPLEX(KIND(1.0D0)) Z 3. Description z This routine evaluates the exponential function e , taking care to avoid machine overflow, and giving a warning if the result cannot be computed to more than half precision. The function is z x evaluated as e =e (cosy+isiny), where x and y are the real and imaginary parts respectively of z. Since cosy and siny are less than or equal to 1 in magnitude, it x x x is possible that e may overflow although e cosy or e siny does x+ln|cosy| not. In this case the alternative formula sign(cosy)e is used for the real part of the result, and x+ln|siny| sign(siny)e for the imaginary part. If either part of the result still overflows, a warning is returned through parameter IFAIL. If Im z is too large, precision may be lost in the evaluation of siny and cosy. Again, a warning is returned through IFAIL. 4. References None. 5. Parameters 1: Z -- COMPLEX(KIND(1.0D0)) Input On entry: the argument z of the function. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 The real part of the result overflows, and is set to the largest safe number with the correct sign. The imaginary part of the result is meaningful. IFAIL= 2 The imaginary part of the result overflows, and is set to the largest safe number with the correct sign. The real part of the result is meaningful. IFAIL= 3 Both real and imaginary parts of the result overflow, and are set to the largest safe number with the correct signs. IFAIL= 4 The computed result is accurate to less than half precision, due to the size of Im z. IFAIL= 5 The computed result has no precision, due to the size of Im z, and is set to zero. 7. Accuracy Accuracy is limited in general only by the accuracy of the Fortran intrinsic functions in the computation of siny, cosy and x e , where x=Re z, y=Im z. As y gets larger, precision will probably be lost due to argument reduction in the evaluation of the sine and cosine functions, until the warning error IFAIL = 4 occurs when y gets larger than \/1/(epsilon), where (epsilon) is the machine precision. Note that on some machines, the intrinsic functions SIN and COS will not operate on arguments larger than about \/1/(epsilon), and so IFAIL can never return as 4. In the comparatively rare event that the result is computed by x+ln|cosy| x+ln|siny| the formulae sign(cosy)e and sign(siny)e , a further small loss of accuracy may be expected due to rounding errors in the logarithmic function. 8. Further Comments None. 9. Example The example program reads values of the argument z from a file, evaluates the function at each value of z and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs13aaf}{NAG On-line Documentation: s13aaf} \beginscroll \begin{verbatim} S13AAF(3NAG) Foundation Library (12/10/92) S13AAF(3NAG) S13 -- Approximations of Special Functions S13AAF S13AAF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S13AAF returns the value of the exponential integral E (x), via 1 the routine name. 2. Specification DOUBLE PRECISION FUNCTION S13AAF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description The routine calculates an approximate value for infty / e-u E (x)= | - du , x>0. 1 / u x For 0 'a T (t)-lnx , where t= -x-1. 1 -- r r 2 r For x>4, -x -x e e --' E (x)= ---y(t)= --- > a T (t), 1 x x -- r r r 11.25-x where t=-1.0+14.5/(x+3.25)= -------. 3.25+x In both cases, -1<=t<=+1. To guard against producing underflows, if x>x the result is set hi directly to zero. For the value x see the Users' Note for your hi implementation. 4. References [1] Abramowitz M and Stegun I A (1965) Handbook of Mathematical Functions. Dover Publications. Ch. 26. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. Constraint: X > 0. 0. 2: IFAIL -- INTEGER Input/Output Before entry, IFAIL must be assigned a value. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. Unless the routine detects an error (see Section 6), IFAIL contains 0 on exit. 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 The routine has been called with an argument less than or equal to zero for which the function is not defined. The result returned is zero. 7. Accuracy If (delta) and (epsilon) are the relative errors in argument and result respectively, then in principle, | -x | | e | |(epsilon)|~=| -----*(delta)| | E (x) | | 1 | so the relative error in the argument is amplified in the result -x by at least a factor e /E (x). The equality should hold if 1 (delta) is greater than the machine precision ((delta) due to data errors etc) but if (delta) is simply a result of round-off in the machine representation, it is possible that an extra figure may be lost in internal calculation and round-off. The behaviour of this amplification factor is shown in Figure 1. Figure 1 Please see figure in printed Reference Manual It should be noted that, for small x, the amplification factor tends to zero and eventually the error in the result will be limited by machine precision. For large x, (epsilon)~x(delta)=(Delta), the absolute error in the argument. 8. Further Comments None. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs13acf}{NAG On-line Documentation: s13acf} \beginscroll \begin{verbatim} S13ACF(3NAG) Foundation Library (12/10/92) S13ACF(3NAG) S13 -- Approximations of Special Functions S13ACF S13ACF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S13ACF returns the value of the cosine integral x / cosu-1 Ci(x)=(gamma)+lnx+ | ------du, x>0 / u 0 via the routine name, where (gamma) denotes Euler's constant. 2. Specification DOUBLE PRECISION FUNCTION S13ACF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description The routine calculates an approximate value for Ci(x). For 0 a T (t) , t=2( --) -1. -- r r ( 16) r=0 For 16 f T (t) and g(x)= > g T (t), t=2( --) -1. -- r r -- r r ( x ) r=0 r=0 For x>=x , Ci(x)=0 to within the accuracy possible (see Section hi 7). 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. Constraint: X > 0. 0. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 The routine has been called with an argument less than or equal to zero for which the function is not defined. The result returned is zero. 7. Accuracy If E and (epsilon) are the absolute and relative errors in the result and (delta) is the relative error in the argument then in principle these are related by | (delta)cosx| |E|~=|(delta)cosx| and |(epsilon)|~=| -----------|. | Ci(x) | That is accuracy will be limited by machine precision near the origin and near the zeros of cosx, but near the zeros of Ci(x) only absolute accuracy can be maintained. The behaviour of this amplification is shown in Figure 1. Figure 1 Please see figure in printed Reference Manual sinx For large values of x, Ci(x)~ ---- therefore x (epsilon)~(delta)xcotx and since (delta) is limited by the finite precision of the machine it becomes impossible to return results which have any relative accuracy. That is, when x>=1/(delta) we have that |Ci(x)|<=1/x~E and hence is not significantly different from zero. Hence x is chosen such that for values of x>=x , Ci(x) in hi hi principle would have values less than the machine precision and so is essentially zero. 8. Further Comments For details of the time taken by the routine see the Users' Note for your implementation. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs13adf}{NAG On-line Documentation: s13adf} \beginscroll \begin{verbatim} S13ADF(3NAG) Foundation Library (12/10/92) S13ADF(3NAG) S13 -- Approximations of Special Functions S13ADF S13ADF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S13ADF returns the value of the sine integral x / sinu Si(x)= | ----du, / u 0 via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S13ADF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description The routine calculates an approximate value for Si(x). For |x|<=16.0 it is based on the Chebyshev expansion --' ( x )2 Si(x)=x > a T (t) , t=2( --) -1. -- r r ( 16) r=0 For 16<|x| f T (t) and g(x)= > g T (t), t=2( --) -1. -- r r -- r r ( x ) r=0 r=0 1 For |x|>=x , Si(x)= -(pi)signx to within machine precision. hi 2 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings There are no failure exits from this routine. The parameter IFAIL has been included for consistency with other routines in this chapter. 7. Accuracy If (delta) and (epsilon) are the relative errors in the argument and result, respectively, then in principle | (delta)sinx| |(epsilon)|~=| -----------|. | Si(x) | The equality may hold if (delta) is greater than the machine precision ((delta) due to data errors etc) but if (delta) is simply due to round-off in the machine representation, then since the factor relating (delta) to (epsilon) is always less than one, the accuracy will be limited by machine precision. 8. Further Comments For details of the time taken by the routine see the Users' Note for your implementation. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs14aaf}{NAG On-line Documentation: s14aaf} \beginscroll \begin{verbatim} S14AAF(3NAG) Foundation Library (12/10/92) S14AAF(3NAG) S14 -- Approximations of Special Functions S14AAF S14AAF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S14AAF returns the value of the Gamma function (Gamma)(x), via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S14AAF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the Gamma function (Gamma)(x). The routine is based on the Chebyshev expansion: --' (Gamma)(1+u)= > a T (t) , where 0<=u<1 , t=2u-1, -- r r r=0 and uses the property (Gamma)(1+x)=x(Gamma)(x). If x=N+1+u where N is integral and 0<=u<1 then it follows that: for N>0 (Gamma)(x)=(x-1)(x-2)...(x-N)(Gamma)(1+u), for N=0 (Gamma)(x)=(Gamma)(1+u), (Gamma)(1+u) for N<0 (Gamma)(x)= ---------------------. x(x+1)(x+2)...(x-N-1) There are four possible failures for this routine: (i) if x is too large, there is a danger of overflow since (Gamma)(x) could become too large to be represented in the machine; (ii) if x is too large and negative, there is a danger of underflow; (iii) if x is equal to a negative integer, (Gamma)(x) would overflow since it has poles at such points; (iv) if x is too near zero, there is again the danger of 1 overflow on some machines. For small x, (Gamma)(x)~= -, and x on some machines there exists a range of non-zero but small values of x for which 1/x is larger than the greatest representable value. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. Constraint: X must not be a negative integer. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 The argument is too large. On soft failure the routine returns the approximate value of (Gamma)(x) at the nearest valid argument. IFAIL= 2 The argument is too large and negative. On soft failure the routine returns zero. IFAIL= 3 The argument is too close to zero. On soft failure the routine returns the approximate value of (Gamma)(x) at the nearest valid argument. IFAIL= 4 The argument is a negative integer, at which value (Gamma)(x) is infinite. On soft failure the routine returns a large positive value. 7. Accuracy Let (delta) and (epsilon) be the relative errors in the argument and the result respectively. If (delta) is somewhat larger than the machine precision (i.e., is due to data errors etc), then (epsilon) and (delta) are approximately related by: (epsilon)~=|x(Psi)(x)|(delta) (provided (epsilon) is also greater than the representation (Gamma)'(x) error). Here (Psi)(x) is the digamma function -----------. (Gamma)(x) Figure 1 shows the behaviour of the error amplification factor |x(Psi)(x)|. Figure 1 Please see figure in printed Reference Manual If (delta) is of the same order as machine precision, then rounding errors could make (epsilon) slightly larger than the above relation predicts. There is clearly a severe, but unavoidable, loss of accuracy for arguments close to the poles of (Gamma)(x) at negative integers. However relative accuracy is preserved near the pole at x=0 right up to the point of failure arising from the danger of overflow. Also accuracy will necessarily be lost as x becomes large since in this region (epsilon)~=(delta)xlnx. However since (Gamma)(x) increases rapidly with x, the routine must fail due to the danger of overflow before this loss of accuracy is too great. (e.g. for x=20, the amplification factor ~= 60.) 8. Further Comments For details of the time taken by the routine see the Users' Note for your implementation. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs14abf}{NAG On-line Documentation: s14abf} \beginscroll \begin{verbatim} S14ABF(3NAG) Foundation Library (12/10/92) S14ABF(3NAG) S14 -- Approximations of Special Functions S14ABF S14ABF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S14ABF returns a value for the logarithm of the Gamma function, ln(Gamma)(x), via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S14ABF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to ln(Gamma)(x). It is based on two Chebyshev expansions. For 0 a T (t), t=2u-1. -- r r r=0 Once (Gamma)(x) has been calculated, the required result is produced by taking the logarithm. For 15.0 b T (t), t=2( --) -1. -- r r ( x ) r=0 For x x there is a danger of setting overflow so the routine vbig must fail. For x<=0.0 the function is not defined and the routine fails. Note: x is calculated so that if xx , big big 1 1 ln(Gamma)(x)=(x- -)lnx-x+ -ln2(pi) 2 2 to within machine accuracy. x is calculated so that vbig ln(Gamma)(x ) is close to the value returned by X02ALF(*). vbig 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. Constraint: X > 0. 0. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 X <= 0.0, the function is undefined. On soft failure, the routine returns zero. IFAIL= 2 X is too large, the function would overflow. On soft failure, the routine returns the value of the function at the largest permissible argument. 7. Accuracy Let (delta) and (epsilon) be the relative errors in the argument and result respectively, and E be the absolute error in the result. If (delta) is somewhat larger than the relative machine precision, then | x*(Psi)(x) | E~=|x*(Psi)(x)|(delta) and (epsilon)~=| ------------|(delta) | ln(Gamma)(x)| (Gamma)'(x) where (Psi)(x) is the digamma function -----------. Figure 1 and (Gamma)(x) Figure 2 show the behaviour of these error amplification factors. Figure 1 Please see figure in printed Reference Manual Figure 2 Please see figure in printed Reference Manual These show that relative error can be controlled, since except near x=1 or 2 relative error is attenuated by the function or at least is not greatly amplified. ( 1 ) For large x, (epsilon)~=(1+ ---)(delta) and for small x, ( lnx) 1 (epsilon)~= ---(delta). lnx The function ln(Gamma)(x) has zeros at x=1 and 2 and hence relative accuracy is not maintainable near those points. However absolute accuracy can still be provided near those zeros as is shown above. If however, (delta) is of the order of the machine precision, then rounding errors in the routine's internal arithmetic may result in errors which are slightly larger than those predicted by the equalities. It should be noted that even in areas where strong attenuation of errors is predicted the relative precision is bounded by the effective machine precision. 8. Further Comments For details of the time taken by the routine see the Users' Note for your implementation. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs14baf}{NAG On-line Documentation: s14baf} \beginscroll \begin{verbatim} S14BAF(3NAG) Foundation Library (12/10/92) S14BAF(3NAG) S14 -- Approximations of Special Functions S14BAF S14BAF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S14BAF computes values for the incomplete gamma functions P(a,x) and Q(a,x). 2. Specification SUBROUTINE S14BAF (A, X, TOL, P, Q, IFAIL) INTEGER IFAIL DOUBLE PRECISION A, X, TOL, P, Q 3. Description This subroutine evaluates the incomplete gamma functions in the normalised form x 1 / a-1 -t P(a,x)= ---------- |t e dt, (Gamma)(a) / 0 infty 1 / a-1 -t Q(a,x)= ---------- | t e dt, (Gamma)(a) / x with x>=0 and a>0, to a user-specified accuracy. With this normalisation, P(a,x)+Q(a,x)=1. Several methods are used to evaluate the functions depending on the arguments a and x, the methods including Taylor expansion for P(a,x), Legendre's continued fraction for Q(a,x), and power series for Q(a,x). When both a and x are large, and a~=x, the uniform asymptotic expansion of Temme [3] is employed for greater efficiency - specifically, this expansion is used when a>=20 and 0.7a<=x<=1.4a. Once either of P or Q is computed, the other is obtained by subtraction from 1. In order to avoid loss of relative precision in this subtraction, the smaller of P and Q is computed first. This routine is derived from subroutine GAM in Gautschi [2]. 4. References [1] Gautschi W (1979) A Computational Procedure for Incomplete Gamma Functions. ACM Trans. Math. Softw. 5 466--481. [2] Gautschi W (1979) Algorithm 542: Incomplete Gamma Functions. ACM Trans. Math. Softw. 5 482--489. [3] Temme N M (1987) On the Computation of the Incomplete Gamma Functions for Large Values of the Parameters. Algorithms for Approximation. (ed J C Mason and M G Cox) Oxford University Press. 5. Parameters 1: A -- DOUBLE PRECISION Input On entry: the argument a of the functions. Constraint: A > 0.0. 2: X -- DOUBLE PRECISION Input On entry: the argument x of the functions. Constraint: X >= 0.0. 3: TOL -- DOUBLE PRECISION Input On entry: the relative accuracy required by the user in the results. If S14BAF is entered with TOL greater than 1.0 or less than machine precision, then the value of machine precision is used instead. 4: P -- DOUBLE PRECISION Output 5: Q -- DOUBLE PRECISION Output On exit: the values of the functions P(a,x) and Q(a,x) respectively. 6: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 On entry A <= 0.0. IFAIL= 2 On entry X < 0.0. IFAIL= 3 Convergence of the Taylor series or Legendre continued fraction fails within 600 iterations. This error is extremely unlikely to occur; if it does, contact NAG. 7. Accuracy There are rare occasions when the relative accuracy attained is somewhat less than that specified by parameter TOL. However, the error should never exceed more than one or two decimal places. Note also that there is a limit of 18 decimal places on the achievable accuracy, because constants in the routine are given to this precision. 8. Further Comments The time taken for a call of S14BAF depends on the precision requested through TOL, and also varies slightly with the input arguments a and x. 9. Example The example program reads values of the argument a and x from a file, evaluates the function and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs15adf}{NAG On-line Documentation: s15adf} \beginscroll \begin{verbatim} S15ADF(3NAG) Foundation Library (12/10/92) S15ADF(3NAG) S15 -- Approximations of Special Functions S15ADF S15ADF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S15ADF returns the value of the complementary error function, erfcx, via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S15ADF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description The routine calculates an approximate value for the complement of the error function 2 infty -u 2 / erfc x= ------ | e du=1-erf x. / \/(pi) x For x>=0, it is based on the Chebyshev expansion 2 -x erfc x=e y(x), --' where y(x)= > a T (t) and t=(x-3.75)/(x+3.75), -1<=t<=+1. -- r r r=0 For x>=x , where there is a danger of setting underflow, the hi result is returned as zero. 2 -x For x<0, erfc x=2-e y(|x|). For x a T (t), where t= -x -1. -- r r 2 r=0 For 2<|x| 'b T (t)}, where t= ---. { -- r r } x+3 { |x|\/(pi) r=0 } For |x|>=x , hi erf x = signx. x is the value above which erf x=+-1 within machine precision. hi Its value is given in the Users' Note for your implementation. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings There are no failure exits from this routine. The parameter IFAIL has been included for consistency with other routines in this chapter. 7. Accuracy On a machine with approximately 11 significant figures the routine agrees with available tables to 10 figures and consistency checking with S15ADF of the relation erf x+erfc x=1.0 shows errors in at worst the 11th figure. 8. Further Comments None. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs17acf}{NAG On-line Documentation: s17acf} \beginscroll \begin{verbatim} S17ACF(3NAG) Foundation Library (12/10/92) S17ACF(3NAG) S17 -- Approximations of Special Functions S17ACF S17ACF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S17ACF returns the value of the Bessel Function Y (x), via the 0 routine name. 2. Specification DOUBLE PRECISION FUNCTION S17ACF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the Bessel Function of the second kind Y (x). 0 Note: Y (x) is undefined for x<=0 and the routine will fail for 0 such arguments. The routine is based on four Chebyshev expansions: For 0 a T (t)+ > b T (t), with t=2( -) -1, 0 (pi) -- r r -- r r ( 8) r=0 r=0 For x>8, / 2 { ( (pi)) ( (pi))} Y (x)= / -----{P (x)sin(x- ----)+Q (x)cos(x- ----)} 0 \/ (pi)x{ 0 ( 4 ) 0 ( 4 )} --' where P (x)= > c T (t), 0 -- r r r=0 2 8 --' ( 8) and Q (x)= - > d T (t), with t=2( -) -1. 0 x -- r r ( x) r=0 2 ( ( x) ) For x near zero, Y (x)~= ----(ln( -)+(gamma)), where (gamma) 0 (pi)( ( 2) ) denotes Euler's constant. This approximation is used when x is sufficiently small for the result to be correct to machine precision. For very large x, it becomes impossible to provide results with any reasonable accuracy (see Section 7), hence the routine fails. Such arguments contain insufficient information to determine the / 2 phase of oscillation of Y (x); only the amplitude, / -, can be 0 \/ x determined and this is returned on soft failure. The range for which this occurs is roughly related to the machine precision: the routine will fail if x>~1/machine precision (see the Users' Note for your implementation for details). 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. [2] Clenshaw C W (1962) Mathematical Tables. Chebyshev Series for Mathematical Functions. HMSO. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. Constraint: X > 0. 0. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 X is too large. On soft failure the routine returns the amplitude of the Y oscillation, \/2/(pi)x. 0 IFAIL= 2 X <= 0.0, Y is undefined. On soft failure the routine 0 returns zero. 7. Accuracy Let (delta) be the relative error in the argument and E be the absolute error in the result. (Since Y (x) oscillates about zero, 0 absolute error and not relative error is significant, except for very small x.) If (delta) is somewhat larger than the machine representation error (e.g. if (delta) is due to data errors etc), then E and (delta) are approximately related by E~=|xY (x)|(delta) 1 (provided E is also within machine bounds). Figure 1 displays the behaviour of the amplification factor |xY (x)|. 1 Figure 1 Please see figure in printed Reference Manual However, if (delta) is of the same order as the machine representation errors, then rounding errors could make E slightly larger than the above relation predicts. For very small x, the errors are essentially independent of (delta) and the routine should provide relative accuracy bounded by the machine precision. For very large x, the above relation ceases to apply. In this / 2 ( (pi)) / 2 region, Y (x)~= / -----sin(x- ----). The amplitude / ----- 0 \/ (pi)x ( 4 ) / (pi)x can be calculated with reasonable accuracy for all x, but ( (pi)) (pi) sin(x- ----) cannot. If x- ---- is written as 2N(pi)+(theta) ( 4 ) 4 ( (pi)) where N is an integer and 0<=(theta)<2(pi), then sin(x- ----) is ( 4 ) -1 determined by (theta) only. If x>~(delta) , (theta) cannot be determined with any accuracy at all. Thus if x is greater than, or of the order of the inverse of machine precision, it is impossible to calculate the phase of Y (x) and the routine must 0 fail. 8. Further Comments For details of the time taken by the routine see the Users' Note for your implementation. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs17adf}{NAG On-line Documentation: s17adf} \beginscroll \begin{verbatim} S17ADF(3NAG) Foundation Library (12/10/92) S17ADF(3NAG) S17 -- Approximations of Special Functions S17ADF S17ADF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S17ADF returns the value of the Bessel Function Y (x), via the 1 routine name. 2. Specification DOUBLE PRECISION FUNCTION S17ADF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the Bessel Function of the second kind Y (x). 1 Note: Y (x) is undefined for x<=0 and the routine will fail for 1 such arguments. The routine is based on four Chebyshev expansions: For 0 a T (t)- ----x+ - > b T (t) , with t 1 (pi) 8 -- r r (pi) 8 -- r r r=0 r=0 2 ( x) =2( -) -1; ( 8) For x>8, / 2 { ( (pi)) ( (pi))} Y (x)= / -----{P (x)sin(x-3 ----)+Q (x)cos(x-3 ----)} 1 \/ (pi)x{ 1 ( 4 ) 1 ( 4 )} --' where P (x)= > c T (t), 1 -- r r r=0 2 8 --' ( 8) and Q (x)= - > d T (t), with t=2( -) -1. 1 x -- r r ( x) r=0 2 For x near zero, Y (x)~=- -----. This approximation is used when 1 (pi)x x is sufficiently small for the result to be correct to machine precision. For extremely small x, there is a danger of overflow 2 in calculating - ----- and for such arguments the routine will (pi)x fail. For very large x, it becomes impossible to provide results with any reasonable accuracy (see Section 7), hence the routine fails. Such arguments contain insufficient information to determine the / 2 phase of oscillation of Y (x), only the amplitude, / -----, 1 \/ (pi)x can be determined and this is returned on soft failure. The range for which this occurs is roughly related to machine precision; the routine will fail if x>~1/machine precision (see the Users' Note for your implementation for details). 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. [2] Clenshaw C W (1962) Mathematical Tables. Chebyshev Series for Mathematical Functions. HMSO. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. Constraint: X > 0. 0. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 X is too large. On soft failure the routine returns the / 2 amplitude of the Y oscillation, / -----. 1 \/ (pi)x IFAIL= 2 X <= 0.0, Y is undefined. On soft failure the routine 1 returns zero. IFAIL= 3 X is too close to zero, there is a danger of overflow. On soft failure, the routine returns the value of Y (x) at the 1 smallest valid argument. 7. Accuracy Let (delta) be the relative error in the argument and E be the absolute error in the result. (Since Y (x) oscillates about zero, 1 absolute error and not relative error is significant, except for very small x.) If (delta) is somewhat larger than the machine precision (e.g. if (delta) is due to data errors etc), then E and (delta) are approximately related by: E~=|xY (x)-Y (x)|(delta) 0 1 (provided E is also within machine bounds). Figure 1 displays the behaviour of the amplification factor |xY (x)-Y (x)|. 0 1 Figure 1 Please see figure in printed Reference Manual However, if (delta) is of the same order as machine precision, then rounding errors could make E slightly larger than the above relation predicts. For very small x, absolute error becomes large, but the relative error in the result is of the same order as (delta). For very large x, the above relation ceases to apply. In this 2(pi) ( 3(pi)) 2(pi) region, Y (x)~= -----sin(x- -----). The amplitude ----- can be 1 x ( 4 ) x ( 3(pi)) calculated with reasonable accuracy for all x, but sin(x- -----) ( 4 ) 3(pi) cannot. If x- ----- is written as 2N(pi)+(theta) where N is an 4 ( 3(pi)) integer and 0<=(theta)<2(pi), then sin(x- -----) is determined by ( 4 ) -1 (theta) only. If x>(delta) , (theta) cannot be determined with any accuracy at all. Thus if x is greater than, or of the order of, the inverse of the machine precision, it is impossible to calculate the phase of Y (x) and the routine must fail. 1 8. Further Comments For details of the time taken by the routine see the Users' Note for your implementation. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs17aef}{NAG On-line Documentation: s17aef} \beginscroll \begin{verbatim} S17AEF(3NAG) Foundation Library (12/10/92) S17AEF(3NAG) S17 -- Approximations of Special Functions S17AEF S17AEF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S17AEF returns the value of the Bessel Function J (x), via the 0 routine name. 2. Specification DOUBLE PRECISION FUNCTION S17AEF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the Bessel Function of the first kind J (x). 0 Note: J (-x)=J (x), so the approximation need only consider x>=0. 0 0 The routine is based on three Chebyshev expansions: For 0 a T (t), with t=2( -) -1. 0 -- r r ( 8) r=0 For x>8, / 2 { ( (pi)) ( (pi))} J (x)= / -----{P (x)cos(x- ----)-Q (x)sin(x- ----)} 0 \/ (pi)x{ 0 ( 4 ) 0 ( 4 )} --' where P (x) = > b T (t), 0 -- r r r=0 2 8 --' ( 8) and Q (x) = - > c T (t), with t=2( -) -1. 0 x -- r r ( x) r=0 For x near zero, J (x)~=1. This approximation is used when x is 0 sufficiently small for the result to be correct to machine precision. For very large x, it becomes impossible to provide results with any reasonable accuracy (see Section 7), hence the routine fails. Such arguments contain insufficient information to determine the / 2 phase of oscillation of J (x); only the amplitude, / -------, 0 \/ (pi)|x| can be determined and this is returned on soft failure. The range for which this occurs is roughly related to the machine precision; the routine will fail if |x|>~1/machine precision (see the Users' Note for your implementation). 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. [2] Clenshaw C W (1962) Mathematical Tables. Chebyshev Series for Mathematical Functions. HMSO. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 X is too large. On soft failure the routine returns the / 2 amplitude of the J oscillation, / -------. 0 \/ (pi)|x| 7. Accuracy Let (delta) be the relative error in the argument and E be the absolute error in the result. (Since J (x) oscillates about zero, 0 absolute error and not relative error is significant.) If (delta) is somewhat larger than the machine precision (e.g. if (delta) is due to data errors etc), then E and (delta) are approximately related by: E~=|xJ (x)|(delta) 1 (provided E is also within machine bounds). Figure 1 displays the behaviour of the amplification factor |xJ (x)|. 1 Figure 1 Please see figure in printed Reference Manual However, if (delta) is of the same order as machine precision, then rounding errors could make E slightly larger than the above relation predicts. For very large x, the above relation ceases to apply. In this / 2 ( (pi)) region, J (x)~= / -------cos(x- ----). The amplitude 0 \/ (pi)|x| ( 4 ) / 2 / ------- can be calculated with reasonable accuracy for all x \/ (pi)|x| ( (pi)) (pi) but cos(x- ----) cannot. If x- ---- is written as ( 4 ) 4 2N(pi)+(theta) where N is an integer and 0 <= (theta) < 2(pi), ( (pi)) -1 then cos(x- ----) is determined by (theta) only. If x>~(delta) , ( 4 ) (theta) cannot be determined with any accuracy at all. Thus if x is greater than, or of the order of, the inverse of the machine precision, it is impossible to calculate the phase of J (x) and 0 the routine must fail. 8. Further Comments For details of the time taken by the routine see the Users' Note for your implementation. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs17aff}{NAG On-line Documentation: s17aff} \beginscroll \begin{verbatim} S17AFF(3NAG) Foundation Library (12/10/92) S17AFF(3NAG) S17 -- Approximations of Special Functions S17AFF S17AFF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S17AFF returns the value of the Bessel Function J (x), via the 1 routine name. 2. Specification DOUBLE PRECISION FUNCTION S17AFF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the Bessel Function of the first kind J (x). 1 Note: J (-x)=-J (x), so the approximation need only consider x>=0 1 1 The routine is based on three Chebyshev expansions: For 0 a T (t) , with t=2( -) -1. 1 8 -- r r ( 8) r=0 For x>8, / 2 { ( 3(pi)) ( 3(pi))} J (x)= / ----x{P (x)cos(x- -----)-Q (x)sin(x- -----)} 1 \/ (pi) { 1 ( 4 ) 1 ( 4 )} --' whereP (x)= > b T (t), 1 -- r r r=0 2 8 --' ( 8) and Q (x)= - > c T (t), with t=2( -) -1. 1 x -- r r ( x) r=0 x For x near zero, J (x)~= -. This approximation is used when x is 1 2 sufficiently small for the result to be correct to machine precision. For very large x, it becomes impossible to provide results with any reasonable accuracy (see Section 7), hence the routine fails. Such arguments contain insufficient information to determine the / 2 phase of oscillation of J (x); only the amplitude, / -------, 1 \/ (pi)|x| can be determined and this is returned on soft failure. The range for which this occurs is roughly related to the machine precision; the routine will fail if |x|>~1/machine precision (see the Users' Note for your implementation for details). 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. [2] Clenshaw C W (1962) Mathematical Tables. Chebyshev Series for Mathematical Functions. HMSO. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 X is too large. On soft failure the routine returns the / 2 amplitude of the J oscillation, / -------. 1 \/ (pi)|x| 7. Accuracy Let (delta) be the relative error in the argument and E be the absolute error in the result. (Since J (x) oscillates about zero, 1 absolute error and not relative error is significant.) If (delta) is somewhat larger than machine precision (e.g. if (delta) is due to data errors etc), then E and (delta) are approximately related by: E~=|xJ (x)-J (x)|(delta) 0 1 (provided E is also within machine bounds). Figure 1 displays the behaviour of the amplification factor |xJ (x)-J (x)|. 0 1 Figure 1 Please see figure in printed Reference Manual However, if (delta) is of the same order as machine precision, then rounding errors could make E slightly larger than the above relation predicts. For very large x, the above relation ceases to apply. In this / 2 ( 3(pi)) region, J (x)~= / -------cos(x- -----). The amplitude 1 \/ (pi)|x| ( 4 ) / 2 / ------- can be calculated with reasonable accuracy for all x \/ (pi)|x| ( 3(pi)) 3(pi) but cos(x- -----) cannot. If x- ----- is written as ( 4 ) 4 2N(pi)+(theta) where N is an integer and 0<=(theta)<2(pi), then ( 3(pi)) -1 cos(x- -----) is determined by (theta) only. If x>~(delta) , ( 4 ) (theta) cannot be determined with any accuracy at all. Thus if x is greater than, or of the order of, machine precision, it is impossible to calculate the phase of J (x) and the routine must 1 fail. 8. Further Comments For details of the time taken by the routine see the Users' Note for your implementation. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs17agf}{NAG On-line Documentation: s17agf} \beginscroll \begin{verbatim} S17AGF(3NAG) Foundation Library (12/10/92) S17AGF(3NAG) S17 -- Approximations of Special Functions S17AGF S17AGF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S17AGF returns a value for the Airy function, Ai(x), via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S17AGF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the Airy function, Ai(x). It is based on a number of Chebyshev expansions: For x<-5, a(t)sinz-b(t)cosz Ai(x)= ----------------- 1/4 (-x) (pi) 2 / 3 where z= ----+ -\/ -x , and a(t) and b(t) are expansions in the 4 3 3 ( 5) variable t=-2( -) -1. ( x) For -5<=x<=0, Ai(x)=f(t)-xg(t), 3 ( x) where f and g are expansions in t=-2( -) -1. ( 5) For 0=9, -z e v(t) Ai(x)= -------, 1/4 x 2 / 3 ( 18) where z= -\/ x and v is an expansion in t=2( --)-1. 3 ( z ) For |x|< the machine precision, the result is set directly to Ai(0). This both saves time and guards against underflow in intermediate calculations. For large negative arguments, it becomes impossible to calculate the phase of the oscillatory function with any precision and so 2/3 ( 3 ) the routine must fail. This occurs if x<-( ----------) , where ( 2(epsilon)) (epsilon) is the machine precision. For large positive arguments, where Ai decays in an essentially exponential manner, there is a danger of underflow so the routine must fail. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 X is too large and positive. On soft failure, the routine returns zero. IFAIL= 2 X is too large and negative. On soft failure, the routine returns zero. 7. Accuracy For negative arguments the function is oscillatory and hence absolute error is the appropriate measure. In the positive region the function is essentially exponential-like and here relative error is appropriate. The absolute error, E, and the relative error, (epsilon), are related in principle to the relative error in the argument, (delta), by | xAi'(x)| E~=|xAi'(x)|(delta), (epsilon)~=| -------|(delta). | Ai(x) | In practice, approximate equality is the best that can be expected. When (delta), (epsilon) or E is of the order of the machine precision, the errors in the result will be somewhat larger. For small x, errors are strongly damped by the function and hence will be bounded by the machine precision. For moderate negative x, the error behaviour is oscillatory but the amplitude of the error grows like 5/4 ( E ) |x| amplitude ( -------)~ ------. ( (delta)) \/(pi) 2 / 3 However the phase error will be growing roughly like -\/ |x| 3 and hence all accuracy will be lost for large negative arguments due to the impossibility of calculating sin and cos to any 2 / 3 1 accuracy if -\/ |x| > -------. 3 (delta) For large positive arguments, the relative error amplification is considerable: (epsilon) / 3 ---------~\/ x . (delta) This means a loss of roughly two decimal places accuracy for arguments in the region of 20. However very large arguments are not possible due to the danger of setting underflow and so the errors are limited in practice. 8. Further Comments None. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs17ahf}{NAG On-line Documentation: s17ahf} \beginscroll \begin{verbatim} S17AHF(3NAG) Foundation Library (12/10/92) S17AHF(3NAG) S17 -- Approximations of Special Functions S17AHF S17AHF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S17AHF returns a value of the Airy function, Bi(x), via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S17AHF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the Airy function Bi(x). It is based on a number of Chebyshev expansions. For x<-5, a(t)cosz+b(t)sinz Bi(x)= -----------------, 1/4 (-x) (pi) 2 / 3 where z= ----+ -\/ -x and a(t) and b(t) are expansions in the 4 3 3 ( 5) variable t=-2( -) -1. ( x) For -5<=x<=0, Bi(x)=\/3(f(t)+xg(t)), 3 ( x) where f and g are expansions in t=-2( -) -1. ( 5) For 0=9, z e u(t) Bi(x)= ------, 1/4 x 2 / 3 ( 18) where z= -\/ x and u is an expansion in t=2( --)-1. 3 ( z ) For |x|< the machine precision, the result is set directly to Bi(0). This both saves time and avoids possible intermediate underflows. For large negative arguments, it becomes impossible to calculate the phase of the oscillating function with any accuracy so the 2/3 ( 3 ) routine must fail. This occurs if x<-( ----------) , where ( 2(epsilon)) (epsilon) is the machine precision. For large positive arguments, there is a danger of causing overflow since Bi grows in an essentially exponential manner, so the routine must fail. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 X is too large and positive. On soft failure, the routine returns zero. IFAIL= 2 X is too large and negative. On soft failure, the routine returns zero. 7. Accuracy For negative arguments the function is oscillatory and hence absolute error is the appropriate measure. In the positive region the function is essentially exponential-like and here relative error is appropriate. The absolute error, E, and the relative error, (epsilon), are related in principle to the relative error in the argument, (delta), by | xBi'(x)| E~=|xBi'(x)|(delta), (epsilon)~=| -------|(delta). | Bi(x) | In practice, approximate equality is the best that can be expected. When (delta), (epsilon) or E is of the order of the machine precision, the errors in the result will be somewhat larger. For small x, errors are strongly damped and hence will be bounded essentially by the machine precision. For moderate to large negative x, the error behaviour is clearly oscillatory but the amplitude of the error grows like amplitude 5/4 ( E ) |x| ( -------)~ ------. ( (delta)) \/(pi) 2 / 3 However the phase error will be growing roughly as -\/ |x| and 3 hence all accuracy will be lost for large negative arguments. This is due to the impossibility of calculating sin and cos to 2 / 3 1 any accuracy if -\/ |x| > -------. 3 (delta) For large positive arguments, the relative error amplification is considerable: (epsilon) / 3 ---------~\/ x . (delta) This means a loss of roughly two decimal places accuracy for arguments in the region of 20. However very large arguments are not possible due to the danger of causing overflow and errors are therefore limited in practice. 8. Further Comments For details of the time taken by the routine see the Users' Note for your implementation. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs17ajf}{NAG On-line Documentation: s17ajf} \beginscroll \begin{verbatim} S17AJF(3NAG) Foundation Library (12/10/92) S17AJF(3NAG) S17 -- Approximations of Special Functions S17AJF S17AJF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S17AJF returns a value of the derivative of the Airy function Ai(x), via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S17AJF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the derivative of the Airy function Ai(x). It is based on a number of Chebyshev expansions. For x<-5, 4 [ b(t) ] Ai'(x)=\/-x[a(t)cosz+ ------sinz], [ (zeta) ] (pi) 2 / 3 where z= ----+(zeta), (zeta)= -\/ -x and a(t) and b(t) are 4 3 3 ( 5) expansions in variable t=-2( -) -1. ( x) For -5<=x<=0, 2 Ai'(x)=x f(t)-g(t), 3 ( x) where f and g are expansions in t=-2( -) -1. ( 5) For 0=9, 4 -z Ai'(x)=\/-xe u(t), 2 / 3 ( 18) where z= -\/ x and u(t) is an expansion in t=2( --)-1. 3 ( z ) For |x|< the square of the machine precision, the result is set directly to Ai'(0). This both saves time and avoids possible intermediate underflows. For large negative arguments, it becomes impossible to calculate a result for the oscillating function with any accuracy and so 4/7 ( ) ( \/(pi) ) the routine must fail. This occurs for x<-( ---------) , where ( (epsilon)) (epsilon) is the machine precision. For large positive arguments, where Ai' decays in an essentially exponential manner, there is a danger of underflow so the routine must fail. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 X is too large and positive. On soft failure, the routine returns zero. IFAIL= 2 X is too large and negative. On soft failure, the routine returns zero. 7. Accuracy For negative arguments the function is oscillatory and hence absolute error is the appropriate measure. In the positive region the function is essentially exponential in character and here relative error is needed. The absolute error, E, and the relative error, (epsilon), are related in principle to the relative error in the argument, (delta), by | 2 | 2 | x Ai(x)| E~=|x Ai(x)|(delta) (epsilon)~=| -------|(delta). | Ai'(x) | In practice, approximate equality is the best that can be expected. When (delta), (epsilon) or E is of the order of the machine precision, the errors in the result will be somewhat larger. For small x, positive or negative, errors are strongly attenuated by the function and hence will be roughly bounded by the machine precision. For moderate to large negative x, the error, like the function, is oscillatory; however the amplitude of the error grows like 7/4 |x| ------. \/(pi) Therefore it becomes impossible to calculate the function with 7/4 \/(pi) any accuracy if |x| > -------. (delta) For large positive x, the relative error amplification is considerable: (epsilon) / 3 ---------~=\/ x . (delta) However, very large arguments are not possible due to the danger of underflow. Thus in practice error amplification is limited. 8. Further Comments None. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs17akf}{NAG On-line Documentation: s17akf} \beginscroll \begin{verbatim} S17AKF(3NAG) Foundation Library (12/10/92) S17AKF(3NAG) S17 -- Approximations of Special Functions S17AKF S17AKF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S17AKF returns a value for the derivative of the Airy function Bi(x), via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S17AKF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine calculates an approximate value for the derivative of the Airy function Bi(x). It is based on a number of Chebyshev expansions. For x<-5, 4 [ b(t) ] Bi'(x)=\/-x[-a(t)sinz+ ------cosz], [ (zeta) ] (pi) 2 / 3 where z= ----+(zeta),(zeta)= -\/ -x and a(t) and b(t) are 4 3 3 ( 5) expansions in the variable t=-2( -) -1. ( x) For -5<=x<=0, 2 Bi'(x)=\/3(x f(t)+g(t)), 3 ( x) where f and g are expansions in t=-2( -) -1. ( 5) For 0=9, 4 z Bi'(x)=\/xe v(t), 2 / 3 ( 18) where z= -\/ x and v(t) is an expansion in t=2( --)-1. 3 ( z ) For |x|< the square of the machine precision, the result is set directly to Bi'(0). This saves time and avoids possible underflows in calculation. For large negative arguments, it becomes impossible to calculate a result for the oscillating function with any accuracy so the 4/7 ( ) ( \/(pi) ) routine must fail. This occurs for x<-( ---------) , where ( (epsilon)) (epsilon) is the machine precision. For large positive arguments, where Bi' grows in an essentially exponential manner, there is a danger of overflow so the routine must fail. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 X is too large and positive. On soft failure the routine returns zero. IFAIL= 2 X is too large and negative. On soft failure the routine returns zero. 7. Accuracy For negative arguments the function is oscillatory and hence absolute error is appropriate. In the positive region the function has essentially exponential behaviour and hence relative error is needed. The absolute error, E, and the relative error (epsilon), are related in principle to the relative error in the argument (delta), by | 2 | 2 | x Bi(x)| E~=|x Bi(x)|(delta) (epsilon)~=| -------|(delta). | Bi'(x) | In practice, approximate equality is the best that can be expected. When (delta), (epsilon) or E is of the order of the machine precision, the errors in the result will be somewhat larger. For small x, positive or negative, errors are strongly attenuated by the function and hence will effectively be bounded by the machine precision. For moderate to large negative x, the error is, like the function, oscillatory. However, the amplitude of the absolute 7/4 |x| error grows like ------. Therefore it becomes impossible to \/(pi) 7/4 \/(pi) calculate the function with any accuracy if |x| > -------. (delta) For large positive x, the relative error amplification is (epsilon) / 3 considerable: ---------~\/ x . However, very large arguments are (delta) not possible due to the danger of overflow. Thus in practice the actual amplification that occurs is limited. 8. Further Comments None. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs17dcf}{NAG On-line Documentation: s17dcf} \beginscroll \begin{verbatim} S17DCF(3NAG) Foundation Library (12/10/92) S17DCF(3NAG) S17 -- Approximations of Special Functions S17DCF S17DCF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S17DCF returns a sequence of values for the Bessel functions Y (z) for complex z, non-negative (nu) and n=0,1,...,N-1, (nu)+n with an option for exponential scaling. 2. Specification SUBROUTINE S17DCF (FNU, Z, N, SCALE, CY, NZ, CWRK, IFAIL) INTEGER N, NZ, IFAIL DOUBLE PRECISION FNU COMPLEX(KIND(1.0D0)) Z, CY(N), CWRK(N) CHARACTER*1 SCALE 3. Description This subroutine evaluates a sequence of values for the Bessel function Y (z), where z is complex, -(pi) < arg z <= (pi), and (nu) (nu) is the real, non-negative order. The N-member sequence is generated for orders (nu), (nu)+1,...,(nu)+N-1. Optionally, the -|Im z| sequence is scaled by the factor e . Note: although the routine may not be called with (nu) less than zero, for negative orders the formula Y (z)=Y (z)cos((pi)(nu))+J (z)sin((pi)(nu)) may be used -(nu) (nu) (nu) (for the Bessel function J (z), see S17DEF). (nu) The routine is derived from the routine CBESY in Amos [2]. It is (1) (2) H (z)-H (z) (nu) (nu) (1) based on the relation Y (z)= -----------------, where H (z) (nu) 2i (nu) (2) and H (z) are the Hankel functions of the first and second (nu) kinds respectively (see S17DLF). When N is greater than 1, extra values of Y (z) are computed (nu) using recurrence relations. For very large |z| or ((nu)+N-1), argument reduction will cause total loss of accuracy, and so no computation is performed. For slightly smaller |z| or ((nu)+N-1), the computation is performed but results are accurate to less than half of machine precision. If |z| is very small, near the machine underflow threshold, or ((nu)+N-1) is too large, there is a risk of overflow and so no computation is performed. In all the above cases, a warning is given by the routine. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. [2] Amos D E (1986) Algorithm 644: A Portable Package for Bessel Functions of a Complex Argument and Nonnegative Order. ACM Trans. Math. Softw. 12 265--273. 5. Parameters 1: FNU -- DOUBLE PRECISION Input On entry: the order, (nu), of the first member of the sequence of functions. Constraint: FNU >= 0.0. 2: Z -- COMPLEX(KIND(1.0D0)) Input On entry: the argument, z, of the functions. Constraint: Z /= (0.0, 0.0). 3: N -- INTEGER Input On entry: the number, N, of members required in the sequence Y (z), Y (z),...,Y (z). Constraint: N >= 1. (nu) (nu)+1 (nu)+N-1 4: SCALE -- CHARACTER*1 Input On entry: the scaling option. If SCALE = 'U', the results are returned unscaled. If SCALE = 'S', the results are returned scaled by the -|Imz| factor e . Constraint: SCALE = 'U' or 'S'. 5: CY(N) -- COMPLEX(KIND(1.0D)) array Output On exit: the N required function values: CY(i) contains Y (z), for i=1,2,...,N. (nu)+i-1 6: NZ -- INTEGER Output On exit: the number of components of CY that are set to zero due to underflow. The positions of such components in the array CY are arbitrary. 7: CWRK(N) -- COMPLEX(KIND(1.0D)) array Workspace 8: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry FNU < 0.0, or Z = (0.0, 0.0), or N < 1, or SCALE /= 'U' or 'S'. IFAIL= 2 No computation has been performed due to the likelihood of overflow, because ABS(Z) is less than a machine-dependent threshold value (given in the Users' Note for your implementation). IFAIL= 3 No computation has been performed due to the likelihood of overflow, because FNU + N - 1 is too large - how large depends on Z as well as the overflow threshold of the machine. IFAIL= 4 The computation has been performed, but the errors due to argument reduction in elementary functions make it likely that the results returned by S17DCF are accurate to less than half of machine precision. This error exit may occur if either ABS(Z) or FNU + N - 1 is greater than a machine- dependent threshold value (given in the Users' Note for your implementation). IFAIL= 5 No computation has been performed because the errors due to argument reduction in elementary functions mean that all precision in results returned by S17DCF would be lost. This error exit may occur if either ABS(Z) or FNU + N - 1 is greater than a machine-dependent threshold value (given in the Users' Note for your implementation). IFAIL= 6 No results are returned because the algorithm termination condition has not been met. This may occur because the parameters supplied to S17DCF would have caused overflow or underflow. 7. Accuracy All constants in subroutine S17DCF are given to approximately 18 digits of precision. Calling the number of digits of precision in the floating-point arithmetic being used t, then clearly the maximum number of correct digits in the results obtained is limited by p=min(t,18). Because of errors in argument reduction when computing elementary functions inside S17DCF, the actual number of correct digits is limited, in general, by p-s, where s~max(1,|log |z||,|log (nu)|) represents the number of digits 10 10 lost due to the argument reduction. Thus the larger the values of |z| and (nu), the less the precision in the result. If S17DCF is called with N>1, then computation of function values via recurrence may lead to some further small loss of accuracy. If function values which should nominally be identical are computed by calls to S17DCF with different base values of (nu) and different N, the computed values may not agree exactly. Empirical tests with modest values of (nu) and z have shown that the discrepancy is limited to the least significant 3-4 digits of precision. 8. Further Comments The time taken by the routine for a call of S17DCF is approximately proportional to the value of N, plus a constant. In general it is much cheaper to call S17DCF with N greater than 1, rather than to make N separate calls to S17DCF. Paradoxically, for some values of z and (nu), it is cheaper to call S17DCF with a larger value of N than is required, and then discard the extra function values returned. However, it is not possible to state the precise circumstances in which this is likely to occur. It is due to the fact that the base value used to start recurrence may be calculated in different regions for different N, and the costs in each region may differ greatly. Note that if the function required is Y (x) or Y (x), i.e., (nu) 0 1 = 0.0 or 1.0, where x is real and positive, and only a single unscaled function value is required, then it may be much cheaper to call S17ACF or S17ADF respectively. 9. Example The example program prints a caption and then proceeds to read sets of data from the input data stream. The first datum is a value for the order FNU, the second is a complex value for the argument, Z, and the third is a value for the parameter SCALE. The program calls the routine with N = 2 to evaluate the function for orders FNU and FNU + 1, and it prints the results. The process is repeated until the end of the input data stream is encountered. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs17def}{NAG On-line Documentation: s17def} \beginscroll \begin{verbatim} S17DEF(3NAG) Foundation Library (12/10/92) S17DEF(3NAG) S17 -- Approximations of Special Functions S17DEF S17DEF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S17DEF returns a sequence of values for the Bessel functions J (z) for complex z, non-negative (nu) and n=0,1,...,N-1, (nu)+n with an option for exponential scaling. 2. Specification SUBROUTINE S17DEF (FNU, Z, N, SCALE, CY, NZ, IFAIL) INTEGER N, NZ, IFAIL DOUBLE PRECISION FNU COMPLEX(KIND(1.0D0)) Z, CY(N) CHARACTER*1 SCALE 3. Description This subroutine evaluates a sequence of values for the Bessel function J (z), where z is complex, -(pi) < arg z <= (pi), and (nu) (nu) is the real, non-negative order. The N-member sequence is generated for orders (nu), (nu)+1,...,(nu)+N-1. Optionally, the -|Im z| sequence is scaled by the factor e . Note: although the routine may not be called with (nu) less than zero, for negative orders the formula J (z)=J (z)cos((pi)(nu))-Y (z)sin((pi)(nu)) may be used -(nu) (nu) (nu) (for the Bessel function Y (z), see S17DCF). (nu) The routine is derived from the routine CBESJ in Amos [2]. It is (nu)(pi)i/2 based on the relations J (z)=e I (-iz), Im z>=0.0 (nu) (nu) -(nu)(pi)i/2 and J (z)=e I (iz), Im z<0.0. (nu) (nu) The Bessel function I (z) is computed using a variety of (nu) techniques depending on the region under consideration. When N is greater than 1, extra values of J (z) are computed (nu) using recurrence relations. For very large |z| or ((nu)+N-1), argument reduction will cause total loss of accuracy, and so no computation is performed. For slightly smaller |z| or ((nu)+N-1), the computation is performed but results are accurate to less than half of machine precision. If Im z is large, there is a risk of overflow and so no computation is performed. In all the above cases, a warning is given by the routine. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. [2] Amos D E (1986) Algorithm 644: A Portable Package for Bessel Functions of a Complex Argument and Nonnegative Order. ACM Trans. Math. Softw. 12 265--273. 5. Parameters 1: FNU -- DOUBLE PRECISION Input On entry: the order, (nu), of the first member of the sequence of functions. Constraint: FNU >= 0.0. 2: Z -- COMPLEX(KIND(1.0D0)) Input On entry: the argument z of the functions. 3: N -- INTEGER Input On entry: the number, N, of members required in the sequence J (z),J (z),...,J (z). Constraint: N >= 1. (nu) (nu)+1 (nu)+N-1 4: SCALE -- CHARACTER*1 Input On entry: the scaling option. If SCALE = 'U', the results are returned unscaled. If SCALE = 'S', the results are returned scaled by the -|Imz| factor e . Constraint: SCALE = 'U' or 'S'. 5: CY(N) -- COMPLEX(KIND(1.0D)) array Output On exit: the N required function values: CY(i) contains J (z), for i=1,2,...,N. (nu)+i-1 6: NZ -- INTEGER Output On exit: the number of components of CY that are set to zero due to underflow. If NZ > 0, then elements CY(N-NZ+1), CY(N- NZ+2),...,CY(N) are set to zero. 7: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry FNU < 0.0, or N < 1, or SCALE /= 'U' or 'S'. IFAIL= 2 No computation has been performed due to the likelihood of overflow, because Im Z is larger than a machine-dependent threshold value (given in the Users' Note for your implementation). This error exit can only occur when SCALE = 'U'. IFAIL= 3 The computation has been performed, but the errors due to argument reduction in elementary functions make it likely that the results returned by S17DEF are accurate to less than half of machine precision. This error exit may occur if either ABS(Z) or FNU + N - 1 is greater than a machine- dependent threshold value (given in the Users' Note for your implementation). IFAIL= 4 No computation has been performed because the errors due to argument reduction in elementary functions mean that all precision in results returned by S17DEF would be lost. This error exit may occur when either ABS(Z) or FNU + N - 1 is greater than a machine-dependent threshold value (given in the Users' Note for your implementation). IFAIL= 5 No results are returned because the algorithm termination condition has not been met. This may occur because the parameters supplied to S17DEF would have caused overflow or underflow. 7. Accuracy All constants in subroutine S17DEF are given to approximately 18 digits of precision. Calling the number of digits of precision in the floating-point arithmetic being used t, then clearly the maximum number of correct digits in the results obtained is limited by p=min(t,18). Because of errors in argument reduction when computing elementary functions inside S17DEF, the actual number of correct digits is limited, in general, by p-s, where s~max(1,|log |(|z|),|log (nu)|) represents the number of digits 10 10 lost due to the argument reduction. Thus the larger the values of |z| and (nu), the less the precision in the result. If S17DEF is called with N>1, then computation of function values via recurrence may lead to some further small loss of accuracy. If function values which should nominally be identical are computed by calls to S17DEF with different base values of (nu) and different N, the computed values may not agree exactly. Empirical tests with modest values of (nu) and z have shown that the discrepancy is limited to the least significant 3-4 digits of precision. 8. Further Comments The time taken by the routine for a call of S17DEF is approximately proportional to the value of N, plus a constant. In general it is much cheaper to call S17DEF with N greater than 1, rather than to make N separate calls to S17DEF. Paradoxically, for some values of z and (nu), it is cheaper to call S17DEF with a larger value of N than is required, and then discard the extra function values returned. However, it is not possible to state the precise circumstances in which this is likely to occur. It is due to the fact that the base value used to start recurrence may be calculated in different regions for different N, and the costs in each region may differ greatly. Note that if the function required is J (x) or J (x), i.e., (nu) 0 1 = 0.0 or 1.0, where x is real and positive, and only a single unscaled function value is required, then it may be much cheaper to call S17AEF or S17AFF respectively. 9. Example The example program prints a caption and then proceeds to read sets of data from the input data stream. The first datum is a value for the order FNU, the second is a complex value for the argument, Z, and the third is a value for the parameter SCALE. The program calls the routine with N = 2 to evaluate the function for orders FNU and FNU + 1, and it prints the results. The process is repeated until the end of the input data stream is encountered. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs17dgf}{NAG On-line Documentation: s17dgf} \beginscroll \begin{verbatim} S17DGF(3NAG) Foundation Library (12/10/92) S17DGF(3NAG) S17 -- Approximations of Special Functions S17DGF S17DGF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S17DGF returns the value of the Airy function Ai(z) or its derivative Ai'(z) for complex z, with an option for exponential scaling. 2. Specification SUBROUTINE S17DGF (DERIV, Z, SCALE, AI, NZ, IFAIL) INTEGER NZ, IFAIL COMPLEX(KIND(1.0D0)) Z, AI CHARACTER*1 DERIV, SCALE 3. Description This subroutine returns a value for the Airy function Ai(z) or its derivative Ai'(z), where z is complex, -(pi) < arg z <= (pi). 2z\/z/3 Optionally, the value is scaled by the factor e . The routine is derived from the routine CAIRY in Amos [2]. It is \/zK (w) -zK (w) 1/3 2/3 based on the relations Ai(z)= ----------, and Ai'(z)= ---------, (pi)\/3 (pi)/3 where K is the modified Bessel function and w=2z\/z/3. (nu) For very large |z|, argument reduction will cause total loss of accuracy, and so no computation is performed. For slightly smaller |z|, the computation is performed but results are accurate to less than half of machine precision. If Re w is too large, and the unscaled function is required, there is a risk of overflow and so no computation is performed. In all the above cases, a warning is given by the routine. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. [2] Amos D E (1986) Algorithm 644: A Portable Package for Bessel Functions of a Complex Argument and Nonnegative Order. ACM Trans. Math. Softw. 12 265--273. 5. Parameters 1: DERIV -- CHARACTER*1 Input On entry: specifies whether the function or its derivative is required. If DERIV = 'F', Ai(z) is returned. If DERIV = 'D', Ai'(z) is returned. Constraint: DERIV = 'F' or 'D'. 2: Z -- COMPLEX(KIND(1.0D0)) Input On entry: the argument z of the function. 3: SCALE -- CHARACTER*1 Input On entry: the scaling option. If SCALE = 'U', the result is returned unscaled. If SCALE = 'S', the result is returned scaled by the factor 2z\/z/3 e . Constraint: SCALE = 'U' or 'S'. 4: AI -- COMPLEX(KIND(1.0D0)) Output On exit: the required function or derivative value. 5: NZ -- INTEGER Output On exit: NZ indicates whether or not AI is set to zero due to underflow. This can only occur when SCALE = 'U'. If NZ = 0, AI is not set to zero. If NZ = 1, AI is set to zero. 6: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry DERIV /= 'F' or 'D'. or SCALE /= 'U' or 'S'. IFAIL= 2 No computation has been performed due to the likelihood of overflow, because Re w is too large, where w=2Z\/Z/3 -- how large depends on Z and the overflow threshold of the machine. This error exit can only occur when SCALE = 'U'. IFAIL= 3 The computation has been performed, but the errors due to argument reduction in elementary functions make it likely that the result returned by S17DGF is accurate to less than half of machine precision. This error exit may occur if ABS (Z) is greater than a machine-dependent threshold value (given in the Users' Note for your implementation). IFAIL= 4 No computation has been performed because the errors due to argument reduction in elementary functions mean that all precision in the result returned by S17DGF would be lost. This error exit may occur if ABS(Z) is greater than a machine-dependent threshold value (given in the Users' Note for your implementation). IFAIL= 5 No result is returned because the algorithm termination condition has not been met. This may occur because the parameters supplied to S17DGF would have caused overflow or underflow. 7. Accuracy All constants in subroutine S17DGF are given to approximately 18 digits of precision. Calling the number of digits of precision in the floating-point arithmetic being used t, then clearly the maximum number of correct digits in the results obtained is limited by p=min(t,18). Because of errors in argument reduction when computing elementary functions inside S17DGF, the actual number of correct digits is limited, in general, by p-s, where s~max(1,|log |z||) represents the number of digits lost due to 10 the argument reduction. Thus the larger the value of |z|, the less the precision in the result. Empirical tests with modest values of z, checking relations between Airy functions Ai(z), Ai'(z), Bi(z) and Bi'(z), have shown errors limited to the least significant 3-4 digits of precision. 8. Further Comments Note that if the function is required to operate on a real argument only, then it may be much cheaper to call S17AGF or S17AJF. 9. Example The example program prints a caption and then proceeds to read sets of data from the input data stream. The first datum is a value for the parameter DERIV, the second is a complex value for the argument, Z, and the third is a value for the parameter SCALE. The program calls the routine and prints the results. The process is repeated until the end of the input data stream is encountered. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs17dhf}{NAG On-line Documentation: s17dhf} \beginscroll \begin{verbatim} S17DHF(3NAG) Foundation Library (12/10/92) S17DHF(3NAG) S17 -- Approximations of Special Functions S17DHF S17DHF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S17DHF returns the value of the Airy function Bi(z) or its derivative Bi'(z) for complex z, with an option for exponential scaling. 2. Specification SUBROUTINE S17DHF (DERIV, Z, SCALE, BI, IFAIL) INTEGER IFAIL COMPLEX(KIND(1.0D0)) Z, BI CHARACTER*1 DERIV, SCALE 3. Description This subroutine returns a value for the Airy function Bi(z) or its derivative Bi'(z), where z is complex, -(pi) < argz <= (pi). |Re (2z\/z/3)| Optionally, the value is scaled by the factor e . The routine is derived from the routine CBIRY in Amos [2]. It is \/z based on the relations Bi(z)= ---(I (w)+I (w)), and -1/3 1/3 \/3 z Bi'(z)= ---(I (w)+I (w)), where I is the modified Bessel -2/3 2/3 (nu) \/3 function and w=2z\/z/3. For very large |z|, argument reduction will cause total loss of accuracy, and so no computation is performed. For slightly smaller |z|, the computation is performed but results are accurate to less than half of machine precision. If Re z is too large, and the unscaled function is required, there is a risk of overflow and so no computation is performed. In all the above cases, a warning is given by the routine. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. [2] Hammersley J M and Handscomb D C (1967) Monte-Carlo Methods. Methuen. 5. Parameters 1: DERIV -- CHARACTER*1 Input On entry: specifies whether the function or its derivative is required. If DERIV = 'F', Bi(z) is returned. If DERIV = 'D', Bi'(z) is returned. Constraint: DERIV = 'F' or 'D'. 2: Z -- COMPLEX(KIND(1.0D0)) Input On entry: the argument z of the function. 3: SCALE -- CHARACTER*1 Input On entry: the scaling option. If SCALE = 'U', the result is returned unscaled. If SCALE = 'S', the result is returned scaled by the |Re(2z\/z/3)| factor e . Constraint: SCALE = 'U' or 'S'. 4: BI -- COMPLEX(KIND(1.0D0)) Output On exit: the required function or derivative value. 5: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry DERIV /= 'F' or 'D'. or SCALE /= 'U' or 'S'. IFAIL= 2 No computation has been performed due to the likelihood of overflow, because real(Z) is too large - how large depends on the overflow threshold of the machine. This error exit can only occur when SCALE = 'U'. IFAIL= 3 The computation has been performed, but the errors due to argument reduction in elementary functions make it likely that the result returned by S17DHF is accurate to less than half of machine precision. This error exit may occur if ABS (Z) is greater than a machine-dependent threshold value (given in the Users' Note for your implementation). IFAIL= 4 No computation has been performed because the errors due to argument reduction in elementary functions mean that all precision in the result returned by S17DHF would be lost. This error exit may occur if ABS(Z) is greater than a machine-dependent threshold value (given in the Users' Note for your implementation). IFAIL= 5 No result is returned because the algorithm termination condition has not been met. This may occur because the parameters supplied to S17DHF would have caused overflow or underflow. 7. Accuracy All constants in subroutine S17DHF are given to approximately 18 digits of precision. Calling the number of digits of precision in the floating-point arithmetic being used t, then clearly the maximum number of correct digits in the results obtained is limited by p=min(t,18). Because of errors in argument reduction when computing elementary functions inside S17DHF, the actual number of correct digits is limited, in general, by p-s, where s~max(1,|log |z||) represents the number of digits lost due to 10 the argument reduction. Thus the larger the value of |z|, the less the precision in the result. Empirical tests with modest values of z, checking relations between Airy functions Ai(z), Ai'(z), Bi(z) and Bi'(z), have shown errors limited to the least significant 3-4 digits of precision. 8. Further Comments Note that if the function is required to operate on a real argument only, then it may be much cheaper to call S17AHF or S17AKF. 9. Example The example program prints a caption and then proceeds to read sets of data from the input data stream. The first datum is a value for the parameter DERIV, the second is a complex value for the argument, Z, and the third is a value for the parameter SCALE. The program calls the routine and prints the results. The process is repeated until the end of the input data stream is encountered. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs17dlf}{NAG On-line Documentation: s17dlf} \beginscroll \begin{verbatim} S17DLF(3NAG) Foundation Library (12/10/92) S17DLF(3NAG) S17 -- Approximations of Special Functions S17DLF S17DLF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S17DLF returns a sequence of values for the Hankel functions (1) (2) H (z) or H (z) for complex z, non-negative (nu) and (nu)+n (nu)+n n=0,1,...,N-1, with an option for exponential scaling. 2. Specification SUBROUTINE S17DLF (M, FNU, Z, N, SCALE, CY, NZ, IFAIL) INTEGER M, N, NZ, IFAIL DOUBLE PRECISION FNU COMPLEX(KIND(1.0D0)) Z, CY(N) CHARACTER*1 SCALE 3. Description This subroutine evaluates a sequence of values for the Hankel (1) (2) function H (z) or H (z), where z is complex, -(pi) < argz (nu) (nu) <= (pi), and (nu) is the real, non-negative order. The N-member sequence is generated for orders (nu), (nu)+1,...,(nu)+N-1. -iz Optionally, the sequence is scaled by the factor e if the (1) iz function is H (z) or by the factor e if the function is (nu) (2) H (z). (nu) Note: although the routine may not be called with (nu) less than zero, for negative orders the formulae (1) (nu)(pi)i (1) (2) -(nu)(pi)i (2) H (z)=e H (z), and H (z)=e H (z) -(nu) (nu) -(nu) (nu) may be used. The routine is derived from the routine CBESH in Amos [2]. It is based on the relation (m) 1 -p(nu) -p H (z)= -e K (ze ), (nu) p (nu) (pi) (pi) where p=i ---- if m=1 and p=-i ---- if m=2, and the Bessel 2 2 function K (z) is computed in the right half-plane only. (nu) Continuation of K (z) to the left half-plane is computed in (nu) terms of the Bessel function I (z). These functions are (nu) evaluated using a variety of different techniques, depending on the region under consideration. (m) When N is greater than 1, extra values of H (z) are computed (nu) using recurrence relations. For very large |z| or ((nu)+N-1), argument reduction will cause total loss of accuracy, and so no computation is performed. For slightly smaller |z| or ((nu)+N-1), the computation is performed but results are accurate to less than half of machine precision. If |z| is very small, near the machine underflow threshold, or ((nu)+N-1) is too large, there is a risk of overflow and so no computation is performed. In all the above cases, a warning is given by the routine. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. [2] Amos D E (1986) Algorithm 644: A Portable Package for Bessel Functions of a Complex Argument and Nonnegative Order. ACM Trans. Math. Softw. 12 265--273. 5. Parameters 1: M -- INTEGER Input On entry: the kind of functions required. (1) If M = 1, the functions are H (z). (nu) (2) If M = 2, the functions are H (z). (nu) Constraint: M = 1 or 2. 2: FNU -- DOUBLE PRECISION Input On entry: the order, (nu), of the first member of the sequence of functions. Constraint: FNU >= 0.0. 3: Z -- COMPLEX(KIND(1.0D0)) Input On entry: the argument z of the functions. Constraint: Z /= (0.0, 0.0). 4: N -- INTEGER Input On entry: the number, N, of members required in the sequence (M) (M) (M) H , H ,..., H . Constraint: N >= 1. (nu) (nu)+1 (nu)+N-1 5: SCALE -- CHARACTER*1 Input On entry: the scaling option. If SCALE = 'U', the results are returned unscaled. If SCALE = 'S', the results are returned scaled by the -iz iz factor e when M = 1, or by the factor e when M = 2. Constraint: SCALE = 'U' or 'S'. 6: CY(N) -- COMPLEX(KIND(1.0D)) array Output On exit: the N required function values: CY(i) contains (M) H , for i=1,2,...,N. (nu)+i-1 7: NZ -- INTEGER Output On exit: the number of components of CY that are set to zero due to underflow. If NZ > 0, then if Imz>0.0 and M = 1, or Imz<0.0 and M = 2, elements CY(1), CY(2),...,CY(NZ) are set to zero. In the complementary half-planes, NZ simply states the number of underflows, and not which elements they are. 8: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry M /= 1 and M /= 2, or FNU < 0.0, or Z = (0.0, 0.0), or N < 1, or SCALE /= 'U' or 'S'. IFAIL= 2 No computation has been performed due to the likelihood of overflow, because ABS(Z) is less than a machine-dependent threshold value (given in the Users' Note for your implementation). IFAIL= 3 No computation has been performed due to the likelihood of overflow, because FNU + N - 1 is too large - how large depends on Z and the overflow threshold of the machine. IFAIL= 4 The computation has been performed, but the errors due to argument reduction in elementary functions make it likely that the results returned by S17DLF are accurate to less than half of machine precision. This error exit may occur if either ABS(Z) or FNU + N - 1 is greater than a machine- dependent threshold value (given in the Users' Note for your implementation). IFAIL= 5 No computation has been performed because the errors due to argument reduction in elementary functions mean that all precision in results returned by S17DLF would be lost. This error exit may occur when either of ABS(Z) or FNU + N - 1 is greater than a machine-dependent threshold value (given in the Users' Note for your implementation). IFAIL= 6 No results are returned because the algorithm termination condition has not been met. This may occur because the parameters supplied to S17DLF would have caused overflow or underflow. 7. Accuracy All constants in subroutine S17DLF are given to approximately 18 digits of precision. Calling the number of digits of precision in the floating-point arithmetic being used t, then clearly the maximum number of correct digits in the results obtained is limited by p=min(t,18). Because of errors in argument reduction when computing elementary functions inside S17DLF, the actual number of correct digits is limited, in general, by p-s, where s~max(1,|log |z||,|log (nu)|) represents the number of digits 10 10 lost due to the argument reduction. Thus the larger the values of |z| and (nu), the less the precision in the result. If S17DLF is called with N>1, then computation of function values via recurrence may lead to some further small loss of accuracy. If function values which should nominally be identical are computed by calls to S17DLF with different base values of (nu) and different N, the computed values may not agree exactly. Empirical tests with modest values of (nu) and z have shown that the discrepancy is limited to the least significant 3-4 digits of precision. 8. Further Comments The time taken by the routine for a call of S17DLF is approximately proportional to the value of N, plus a constant. In general it is much cheaper to call S17DLF with N greater than 1, rather than to make N separate calls to S17DLF. Paradoxically, for some values of z and (nu), it is cheaper to call S17DLF with a larger value of N than is required, and then discard the extra function values returned. However, it is not possible to state the precise circumstances in which this is likely to occur. It is due to the fact that the base value used to start recurrence may be calculated in different regions for different N, and the costs in each region may differ greatly. 9. Example The example program prints a caption and then proceeds to read sets of data from the input data stream. The first datum is a value for the kind of function, M, the second is a value for the order FNU, the third is a complex value for the argument, Z, and the fourth is a value for the parameter SCALE. The program calls the routine with N = 2 to evaluate the function for orders FNU and FNU + 1, and it prints the results. The process is repeated until the end of the input data stream is encountered. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs18acf}{NAG On-line Documentation: s18acf} \beginscroll \begin{verbatim} S18ACF(3NAG) Foundation Library (12/10/92) S18ACF(3NAG) S18 -- Approximations of Special Functions S18ACF S18ACF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S18ACF returns the value of the modified Bessel Function K (x), 0 via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S18ACF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the modified Bessel Function of the second kind K (x). 0 Note: K (x) is undefined for x<=0 and the routine will fail for 0 such arguments. The routine is based on five Chebyshev expansions: For 0 a T (t)+ > b T (t) , where t=2x -1; 0 -- r r -- r r r=0 r=0 For 1 c T (t) , where t=2x-3; 0 -- r r r=0 For 2 d T (t) , where t=x-3; 0 -- r r r=0 For x>4, -x e --' 9-x K (x)= --- > e T (t) , where t= ---. 0 -- r r 1+x \/x r=0 ( x) For x near zero, K (x)~=-(gamma)-ln( -), where (gamma) denotes 0 ( 2) Euler's constant. This approximation is used when x is sufficiently small for the result to be correct to machine precision. For large x, where there is a danger of underflow due to the smallness of K , the result is set exactly to zero. 0 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. Constraint: X > 0. 0. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 X <= 0.0, K is undefined. On soft failure the routine 0 returns zero. 7. Accuracy Let (delta) and (epsilon) be the relative errors in the argument and result respectively. If (delta) is somewhat larger than the machine precision (i.e., if (delta) is due to data errors etc), then (epsilon) and (delta) are approximately related by: | xK (x)| | 1 | (epsilon)~=| ------|(delta). | K (x) | | 0 | Figure 1 shows the behaviour of the error amplification factor | xK (x)| | 1 | | ------|. | K (x) | | 0 | Figure 1 Please see figure in printed Reference Manual However, if (delta) is of the same order as machine precision, then rounding errors could make (epsilon) slightly larger than the above relation predicts. | 1 | For small x, the amplification factor is approximately | ---|, | lnx| which implies strong attenuation of the error, but in general (epsilon) can never be less than the machine precision. For large x, (epsilon)~=x(delta) and we have strong amplification of the relative error. Eventually K , which is asymptotically 0 -x e given by ---, becomes so small that it cannot be calculated \/x without underflow and hence the routine will return zero. Note that for large x the errors will be dominated by those of the Fortran intrinsic function EXP. 8. Further Comments For details of the time taken by the routine see the appropriate the Users' Note for your implementation. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs18adf}{NAG On-line Documentation: s18adf} \beginscroll \begin{verbatim} S18ADF(3NAG) Foundation Library (12/10/92) S18ADF(3NAG) S18 -- Approximations of Special Functions S18ADF S18ADF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S18ADF returns the value of the modified Bessel Function K (x), 1 via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S18ADF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the modified Bessel Function of the second kind K (x). 1 Note: K (x) is undefined for x<=0 and the routine will fail for 1 such arguments. The routine is based on five Chebyshev expansions: For 0 a T (t)-x > b T (t) , where t=2x -1; 1 x -- r r -- r r r=0 r=0 For 1 c T (t) , where t=2x-3; 1 -- r r r=0 For 2 d T (t) , where t=x-3; 1 -- r r r=0 For x>4, -x e --' 9-x K (x)= --- > e T (t) , where t= ---. 1 -- r r 1+x \/x r=0 1 For x near zero, K (x)~= -. This approximation is used when x is 1 x sufficiently small for the result to be correct to machine precision. For very small x on some machines, it is impossible to 1 calculate - without overflow and the routine must fail. x For large x, where there is a danger of underflow due to the smallness of K , the result is set exactly to zero. 1 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. Constraint: X > 0. 0. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 X <= 0.0, K is undefined. On soft failure the routine 1 returns zero. IFAIL= 2 X is too small, there is a danger of overflow. On soft failure the routine returns approximately the largest representable value. 7. Accuracy Let (delta) and (epsilon) be the relative errors in the argument and result respectively. If (delta) is somewhat larger than the machine precision (i.e., if (delta) is due to data errors etc), then (epsilon) and (delta) are approximately related by: | xK (x)-K (x)| | 0 1 | (epsilon)~=| ------------|(delta). | K (x) | | 1 | Figure 1 shows the behaviour of the error amplification factor | xK (x)-K (x)| | 0 1 | | ------------|. | K (x) | | 1 | Figure 1 Please see figure in printed Reference Manual However if (delta) is of the same order as the machine precision, then rounding errors could make (epsilon) slightly larger than the above relation predicts. For small x, (epsilon)~=(delta) and there is no amplification of errors. For large x, (epsilon)~=x(delta) and we have strong amplification of the relative error. Eventually K , which is asymptotically 1 -x e given by ---, becomes so small that it cannot be calculated \/x without underflow and hence the routine will return zero. Note that for large x the errors will be dominated by those of the Fortran intrinsic function EXP. 8. Further Comments For details of the time taken by the routine see the Users' Note for your implementation. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs18aef}{NAG On-line Documentation: s18aef} \beginscroll \begin{verbatim} S18AEF(3NAG) Foundation Library (12/10/92) S18AEF(3NAG) S18 -- Approximations of Special Functions S18AEF S18AEF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S18AEF returns the value of the modified Bessel Function I (x), 0 via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S18AEF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the modified Bessel Function of the first kind I (x). 0 Note: I (-x)=I (x), so the approximation need only consider x>=0. 0 0 The routine is based on three Chebyshev expansions: For 0 a T (t) , where t=2( -)-1; 0 -- r r ( 4) r=0 For 4 b T (t) , where t= ---; 0 -- r r 4 r=0 For x>12, x e --' ( 12) I (x)= --- > c T (t) , where t=2( --)-1. 0 -- r r ( x ) \/x r=0 For small x, I (x)~=1. This approximation is used when x is 0 sufficiently small for the result to be correct to machine precision. For large x, the routine must fail because of the danger of x overflow in calculating e . 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 X is too large. On soft failure the routine returns the approximate value of I (x) at the nearest valid argument. 0 7. Accuracy Let (delta) and (epsilon) be the relative errors in the argument and result respectively. If (delta) is somewhat larger than the machine precision (i.e., if (delta) is due to data errors etc), then (epsilon) and (delta) are approximately related by: | xI (x)| | 1 | (epsilon)~=| ------|(delta). | I (x) | | 0 | Figure 1 shows the behaviour of the error amplification factor | xI (x)| | 1 | | ------|. | I (x) | | 0 | Figure 1 Please see figure in printed Reference Manual However if (delta) is of the same order as machine precision, then rounding errors could make (epsilon) slightly larger than the above relation predicts. 2 x For small x the amplification factor is approximately --, which 2 implies strong attenuation of the error, but in general (epsilon) can never be less than the machine precision. For large x, (epsilon)~=x(delta) and we have strong amplification of errors. However the routine must fail for quite moderate values of x, because I (x) would overflow; hence in practice the 0 loss of accuracy for large x is not excessive. Note that for large x the errors will be dominated by those of the Fortran intrinsic function EXP. 8. Further Comments For details of the time taken by the routine see the Users' Note for your implementation. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs18aff}{NAG On-line Documentation: s18aff} \beginscroll \begin{verbatim} S18AFF(3NAG) Foundation Library (12/10/92) S18AFF(3NAG) S18 -- Approximations of Special Functions S18AFF S18AFF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S18AFF returns a value for the modified Bessel Function I (x), 1 via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S18AFF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the modified Bessel Function of the first kind I (x). 1 Note: I (-x)=-I (x), so the approximation need only consider x>=0 1 1 The routine is based on three Chebyshev expansions: For 0 a T (t) , where t=2( -) -1; 1 -- r r ( 4) r=0 For 4 b T (t) , where t= ---; 1 -- r r 4 r=0 For x>12, x e --' ( 12) I (x)= --- > c T (t) , where t=2( --)-1. 1 -- r r ( x ) \/x r=0 For small x, I (x)~=x. This approximation is used when x is 1 sufficiently small for the result to be correct to machine precision. For large x, the routine must fail because I (x) cannot be 1 represented without overflow. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 X is too large. On soft failure the routine returns the approximate value of I (x) at the nearest valid argument. 1 7. Accuracy Let (delta) and (epsilon) be the relative errors in the argument and result respectively. If (delta) is somewhat larger than the machine precision (i.e., if (delta) is due to data errors etc), then (epsilon) and (delta) are approximately related by: | xI (x)-I (x)| | 0 1 | (epsilon)~=| ------------|(delta). | I (x) | | 1 | Figure 1 shows the behaviour of the error amplification factor | xI (x)-I (x)| | 0 1 | | ------------|, | I (x) | | 1 | Figure 1 Please see figure in printed Reference Manual However if (delta) is of the same order as machine precision, then rounding errors could make (epsilon) slightly larger than the above relation predicts. For small x, (epsilon)~=(delta) and there is no amplification of errors. For large x, (epsilon)~=x(delta) and we have strong amplification of errors. However the routine must fail for quite moderate values of x because I (x) would overflow; hence in practice the 1 loss of accuracy for large x is not excessive. Note that for large x, the errors will be dominated by those of the Fortran intrinsic function EXP. 8. Further Comments For details of the time taken by the routine see the Users' Note for your implementation. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs18dcf}{NAG On-line Documentation: s18dcf} \beginscroll \begin{verbatim} S18DCF(3NAG) Foundation Library (12/10/92) S18DCF(3NAG) S18 -- Approximations of Special Functions S18DCF S18DCF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S18DCF returns a sequence of values for the modified Bessel functions K (z) for complex z, non-negative (nu) and (nu)+n n=0,1,...,N-1, with an option for exponential scaling. 2. Specification SUBROUTINE S18DCF (FNU, Z, N, SCALE, CY, NZ, IFAIL) INTEGER N, NZ, IFAIL DOUBLE PRECISION FNU COMPLEX(KIND(1.0D0)) Z, CY(N) CHARACTER*1 SCALE 3. Description This subroutine evaluates a sequence of values for the modified Bessel function K (z), where z is complex, -(pi) < arg z <= (nu) (pi), and (nu) is the real, non-negative order. The N-member sequence is generated for orders (nu), (nu)+1,...,(nu)+N-1. z Optionally, the sequence is scaled by the factor e . The routine is derived from the routine CBESK in Amos [2]. Note: although the routine may not be called with (nu) less than zero, for negative orders the formula K (z)=K (z) may be -(nu) (nu) used. When N is greater than 1, extra values of K (z) are computed (nu) using recurrence relations. For very large |z| or ((nu)+N-1), argument reduction will cause total loss of accuracy, and so no computation is performed. For slightly smaller |z| or ((nu)+N-1), the computation is performed but results are accurate to less than half of machine precision. If |z| is very small, near the machine underflow threshold, or ((nu)+N-1) is too large, there is a risk of overflow and so no computation is performed. In all the above cases, a warning is given by the routine. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. [2] Amos D E (1986) Algorithm 644: A Portable Package for Bessel Functions of a Complex Argument and Nonnegative Order. ACM Trans. Math. Softw. 12 265--273. 5. Parameters 1: FNU -- DOUBLE PRECISION Input On entry: the order, (nu), of the first member of the sequence of functions. Constraint: FNU >= 0.0. 2: Z -- COMPLEX(KIND(1.0D0)) Input On entry: the argument z of the functions. Constraint: Z /= (0.0, 0.0). 3: N -- INTEGER Input On entry: the number, N, of members required in the sequence K (z), K (z),..., K (z). Constraint: N >= 1. (nu) (nu)+1 (nu)+N-1 4: SCALE -- CHARACTER*1 Input On entry: the scaling option. If SCALE = 'U', the results are returned unscaled. If SCALE = 'S', the results are returned scaled by the z factor e . Constraint: SCALE = 'U' or 'S'. 5: CY(N) -- COMPLEX(KIND(1.0D)) array Output On exit: the N required function values: CY(i) contains K (z), for i=1,2,...,N. (nu)+i-1 6: NZ -- INTEGER Output On exit: the number of components of CY that are set to zero due to underflow. If NZ > 0 and Rez>=0.0, elements CY(1),CY (2),...,CY(NZ) are set to zero. If Rez<0.0, NZ simply states the number of underflows, and not which elements they are. 7: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry FNU < 0.0, or Z = (0.0, 0.0), or N < 1, or SCALE /= 'U' or 'S'. IFAIL= 2 No computation has been performed due to the likelihood of overflow, because ABS(Z) is less than a machine-dependent threshold value (given in the Users' Note for your implementation). IFAIL= 3 No computation has been performed due to the likelihood of overflow, because FNU + N - 1 is too large - how large depends on Z and the overflow threshold of the machine. IFAIL= 4 The computation has been performed, but the errors due to argument reduction in elementary functions make it likely that the results returned by S18DCF are accurate to less than half of machine precision. This error exit may occur if either ABS(Z) or FNU + N - 1 is greater than a machine- dependent threshold value (given in the Users' Note for your implementation). IFAIL= 5 No computation has been performed because the errors due to argument reduction in elementary functions mean that all precision in results returned by S18DCF would be lost. This error exit may occur when either ABS(Z) or FNU + N - 1 is greater than a machine-dependent threshold value (given in the Users' Note for your implementation). IFAIL= 6 No results are returned because the algorithm termination condition has not been met. This may occur because the parameters supplied to S18DCF would have caused overflow or underflow. 7. Accuracy All constants in subroutine S18DCF are given to approximately 18 digits of precision. Calling the number of digits of precision in the floating-point arithmetic being used t, then clearly the maximum number of correct digits in the results obtained is limited by p=min(t,18). Because of errors in argument reduction when computing elementary functions inside S18DCF, the actual number of correct digits is limited, in general, by p-s, where s~max(1,|log |z||,|log (nu)|) represents the number of digits 10 10 lost due to the argument reduction. Thus the larger the values of |z| and (nu), the less the precision in the result. If S18DCF is called with N>1, then computation of function values via recurrence may lead to some further small loss of accuracy. If function values which should nominally be identical are computed by calls to S18DCF with different base values of (nu) and different N, the computed values may not agree exactly. Empirical tests with modest values of (nu) and z have shown that the discrepancy is limited to the least significant 3-4 digits of precision. 8. Further Comments The time taken by the routine for a call of S18DCF is approximately proportional to the value of N, plus a constant. In general it is much cheaper to call S18DCF with N greater than 1, rather than to make N separate calls to S18DCF. Paradoxically, for some values of z and (nu), it is cheaper to call S18DCF with a larger value of N than is required, and then discard the extra function values returned. However, it is not possible to state the precise circumstances in which this is likely to occur. It is due to the fact that the base value used to start recurrence may be calculated in different regions for different N, and the costs in each region may differ greatly. Note that if the function required is K (x) or K (x), i.e., 0 1 (nu)=0.0 or 1.0, where x is real and positive, and only a single function value is required, then it may be much cheaper to call S18ACF, S18ADF, S18CCF(*) or S18CDF(*), depending on whether a scaled result is required or not. 9. Example The example program prints a caption and then proceeds to read sets of data from the input data stream. The first datum is a value for the order FNU, the second is a complex value for the argument, Z, and the third is a value for the parameter SCALE. The program calls the routine with N = 2 to evaluate the function for orders FNU and FNU + 1, and it prints the results. The process is repeated until the end of the input data stream is encountered. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs18def}{NAG On-line Documentation: s18def} \beginscroll \begin{verbatim} S18DEF(3NAG) Foundation Library (12/10/92) S18DEF(3NAG) S18 -- Approximations of Special Functions S18DEF S18DEF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S18DEF returns a sequence of values for the modified Bessel functions I (z) for complex z, non-negative (nu) and (nu)+n n=0,1,...,N-1, with an option for exponential scaling. 2. Specification SUBROUTINE S18DEF (FNU, Z, N, SCALE, CY, NZ, IFAIL) INTEGER N, NZ, IFAIL DOUBLE PRECISION FNU COMPLEX(KIND(1.0D0)) Z, CY(N) CHARACTER*1 SCALE 3. Description This subroutine evaluates a sequence of values for the modified Bessel function I (z), where z is complex, -(pi) < argz <= (nu) (pi), and (nu) is the real, non-negative order. The N-member sequence is generated for orders (nu), (nu)+1,...,(nu)+N-1. -|Rez| Optionally, the sequence is scaled by the factor e . The routine is derived from the routine CBESI in Amos [2]. Note: although the routine may not be called with (nu) less than zero, for negative orders the formula 2 I (z)=I (z)+ ----sin((pi)(nu))K (z) may be used (for -(nu) (nu) (pi) (nu) the Bessel function K (z), see S18DCF). (nu) When N is greater than 1, extra values of I (z) are computed (nu) using recurrence relations. For very large |z| or ((nu)+N-1), argument reduction will cause total loss of accuracy, and so no computation is performed. For slightly smaller |z| or ((nu)+N-1), the computation is performed but results are accurate to less than half of machine precision. If Re(z) is too large and the unscaled function is required, there is a risk of overflow and so no computation is performed. In all the above cases, a warning is given by the routine. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. [2] Amos D E (1986) Algorithm 644: A Portable Package for Bessel Functions of a Complex Argument and Nonnegative Order. ACM Trans. Math. Softw. 12 265--273. 5. Parameters 1: FNU -- DOUBLE PRECISION Input On entry: the order, (nu), of the first member of the sequence of functions. Constraint: FNU >= 0.0. 2: Z -- COMPLEX(KIND(1.0D0)) Input On entry: the argument z of the functions. 3: N -- INTEGER Input On entry: the number, N, of members required in the sequence I (z),I (z),...,I (z). Constraint: N >= 1. (nu) (nu)+1 (nu)+N-1 4: SCALE -- CHARACTER*1 Input On entry: the scaling option. If SCALE = 'U', the results are returned unscaled. If SCALE = 'S', the results are returned scaled by the -|Rez| factor e . Constraint: SCALE = 'U' or 'S'. 5: CY(N) -- COMPLEX(KIND(1.0D)) array Output On exit: the N required function values: CY(i) contains I (z), for i=1,2,...,N. (nu)+i-1 6: NZ -- INTEGER Output On exit: the number of components of CY that are set to zero due to underflow. If NZ > 0, then elements CY(N-NZ+1),CY(N-NZ+2),...,CY(N) are set to zero. 7: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry FNU < 0.0, or N < 1, or SCALE /= 'U' or 'S'. IFAIL= 2 No computation has been performed due to the likelihood of overflow, because real(Z) is greater than a machine- dependent threshold value (given in the Users' Note for your implementation). This error exit can only occur when SCALE = 'U'. IFAIL= 3 The computation has been performed, but the errors due to argument reduction in elementary functions make it likely that the results returned by S18DEF are accurate to less than half of machine precision. This error exit may occur when either ABS(Z) or FNU + N - 1 is greater than a machine- dependent threshold value (given in the Users' Note for your implementation). IFAIL= 4 No computation has been performed because the errors due to argument reduction in elementary functions mean that all precision in results returned by S18DEF would be lost. This error exit may occur when either ABS(Z) or FNU + N - 1 is greater than a machine-dependent threshold value (given in the Users' Note for your implementation). IFAIL= 5 No results are returned because the algorithm termination condition has not been met. This may occur because the parameters supplied to S18DEF would have caused overflow or underflow. 7. Accuracy All constants in subroutine S18DEF are given to approximately 18 digits of precision. Calling the number of digits of precision in the floating-point arithmetic being used t, then clearly the maximum number of correct digits in the results obtained is limited by p=min(t,18). Because of errors in argument reduction when computing elementary functions inside S18DEF, the actual number of correct digits is limited, in general, by p-s, where s~max(1,|log |z||,|log (nu)|) represents the number of digits 10 10 lost due to the argument reduction. Thus the larger the values of |z| and (nu), the less the precision in the result. If S18DEF is called with N>1, then computation of function values via recurrence may lead to some further small loss of accuracy. If function values which should nominally be identical are computed by calls to S18DEF with different base values of (nu) and different N, the computed values may not agree exactly. Empirical tests with modest values of (nu) and z have shown that the discrepancy is limited to the least significant 3-4 digits of precision. 8. Further Comments The time taken by the routine for a call of S18DEF is approximately proportional to the value of N, plus a constant. In general it is much cheaper to call S18DEF with N greater than 1, rather than to make N separate calls to S18DEF. Paradoxically, for some values of z and (nu), it is cheaper to call S18DEF with a larger value of N than is required, and then discard the extra function values returned. However, it is not possible to state the precise circumstances in which this is likely to occur. It is due to the fact that the base value used to start recurrence may be calculated in different regions for different N, and the costs in each region may differ greatly. Note that if the function required is I (x) or I (x), i.e., 0 1 (nu)=0.0 or 1.0, where x is real and positive, and only a single function value is required, then it may be much cheaper to call S18AEF, S18AFF, S18CEF(*) or S18CFF(*), depending on whether a scaled result is required or not. 9. Example The example program prints a caption and then proceeds to read sets of data from the input data stream. The first datum is a value for the order FNU, the second is a complex value for the argument, Z, and the third is a value for the parameter SCALE. The program calls the routine with N = 2 to evaluate the function for orders FNU and FNU + 1, and it prints the results. The process is repeated until the end of the input data stream is encountered. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs19aaf}{NAG On-line Documentation: s19aaf} \beginscroll \begin{verbatim} S19AAF(3NAG) Foundation Library (12/10/92) S19AAF(3NAG) S19 -- Approximations of Special Functions S19AAF S19AAF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S19AAF returns a value for the Kelvin function ber x via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S19AAF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the Kelvin function berx. Note: ber(-x)=berx, so the approximation need only consider x>=0.0. The routine is based on several Chebyshev expansions: For 0<=x<=5, --' ( x)4 berx= > a T (t) with t=2( -) -1; -- r r ( 5) r=0 For x>5, x/\/2 e [( 1 ) 1 ] berx= --------[(1+ -a(t))cos(alpha)+ -b(t)sin(alpha)] [( x ) x ] \/2(pi)x -x/\/2 e [( 1 ) 1 ] + --------[(1+ -c(t))sin(beta)+ -d(t)cos(beta)] [( x ) x ] \/2(pi)x x (pi) x (pi) where (alpha)= ---- ----, (beta)= ---+ ----, 8 8 \/2 /2 and a(t), b(t), c(t), and d(t) are expansions in the variable 10 t= ---1. x When x is sufficiently close to zero, the result is set directly to ber 0=1.0. For large x, there is a danger of the result being totally inaccurate, as the error amplification factor grows in an essentially exponential manner; therefore the routine must fail. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 On entry ABS(X) is too large for an accurate result to be returned. On soft failure, the routine returns zero. 7. Accuracy Since the function is oscillatory, the absolute error rather than the relative error is important. Let E be the absolute error in the result and (delta) be the relative error in the argument. If (delta) is somewhat larger than the machine precision, then we have: | x | E~=| ---(ber x+bei x)|(delta) | 1 1 | | \/2 | (provided E is within machine bounds). For small x the error amplification is insignificant and thus the absolute error is effectively bounded by the machine precision. For medium and large x, the error behaviour is oscillatory and / x x/\/2 its amplitude grows like / -----e . Therefore it is not \/ 2(pi) possible to calculate the function with any accuracy when x/\/2 /2(pi) \/xe > -------. Note that this value of x is much smaller than (delta) the minimum value of x for which the function overflows. 8. Further Comments For details of the time taken by the routine see the Users' Note for your implementation. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs19abf}{NAG On-line Documentation: s19abf} \beginscroll \begin{verbatim} S19ABF(3NAG) Foundation Library (12/10/92) S19ABF(3NAG) S19 -- Approximations of Special Functions S19ABF S19ABF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S19ABF returns a value for the Kelvin function bei x via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S19ABF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the Kelvin function beix. Note: bei(-x)=beix, so the approximation need only consider x>=0.0. The routine is based on several Chebyshev expansions: For 0<=x<=5, 2 x --' ( x)4 bei x= -- > a T (t), with t=2( -) -1; 4 -- r r ( 5) r=0 For x>5, x/\/2 e [( 1 ) 1 ] bei x= --------[(1+ -a(t))sin(alpha)- -b(t)cos(alpha)] [( x ) x ] \/2(pi)x x/\/2 e [( 1 ) 1 ] + --------[(1+ -c(t))cos(beta)- -d(t)sin(beta)] [( x ) x ] \/2(pi)x x (pi) x (pi) where(alpha)= ---- ----,(beta)= ---+ ----, 8 8 \/2 /2 and a(t), b(t), c(t), and d(t) are expansions in the variable 10 t= ---1. x When x is sufficiently close to zero, the result is computed as 2 x bei x= --. If this result would underflow, the result returned is 4 beix=0.0. For large x, there is a danger of the result being totally inaccurate, as the error amplification factor grows in an essentially exponential manner; therefore the routine must fail. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 On entry ABS(X) is too large for an accurate result to be returned. On soft failure, the routine returns zero. 7. Accuracy Since the function is oscillatory, the absolute error rather than the relative error is important. Let E be the absolute error in the function, and (delta) be the relative error in the argument. If (delta) is somewhat larger than the machine precision, then we have: | x | E~=| ---(-ber x+bei x)|(delta) | 1 1 | | \/2 | (provided E is within machine bounds). For small x the error amplification is insignificant and thus the absolute error is effectively bounded by the machine precision. For medium and large x, the error behaviour is oscillatory and / x x\/2 its amplitude grows like / -----e . Therefore it is \/ 2(pi) impossible to calculate the functions with any accuracy when x/\/2 /2(pi) \/xe > -------. Note that this value of x is much smaller than (delta) the minimum value of x for which the function overflows. 8. Further Comments For details of the time taken by the routine see the Users' Note for your implementation. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs19acf}{NAG On-line Documentation: s19acf} \beginscroll \begin{verbatim} S19ACF(3NAG) Foundation Library (12/10/92) S19ACF(3NAG) S19 -- Approximations of Special Functions S19ACF S19ACF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S19ACF returns a value for the Kelvin function ker x, via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S19ACF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the Kelvin function ker x. Note: for x<0 the function is undefined and at x=0 it is infinite so we need only consider x>0. The routine is based on several Chebyshev expansions: For 03, -x/\/2 / (pi) [( 1 ) 1 ] ker x= / ----e [(1+ -c(t))cos(beta)- -d(t)sin(beta)] \/ 2x [( x ) x ] x (pi) where (beta)= ---+ ----, and c(t) and d(t) are expansions in the 8 \/2 6 variable t= --1. x When x is sufficiently close to zero, the result is computed as 2 ( x) ( 3 2) x ker x=-(gamma)-log( -)+((pi)- -x ) -- ( 2) ( 8 ) 16 and when x is even closer to zero, simply as ( x) ker x=-(gamma)-log( -). ( 2) / (pi) -x/\/2 For large x, ker x is asymptotically given by / ----e and \/ 2x this becomes so small that it cannot be computed without underflow and the routine fails. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. Constraint: X > 0. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 On entry X is too large, the result underflows. On soft failure, the routine returns zero. IFAIL= 2 On entry X <= 0, the function is undefined. On soft failure the routine returns zero. 7. Accuracy Let E be the absolute error in the result, (epsilon) be the relative error in the result and (delta) be the relative error in the argument. If (delta) is somewhat larger than the machine precision, then we have: | x | E~=| ---(ker x+kei x)|(delta), | 1 1 | | \/2 | | ker x+kei x| | x 1 1 | (epsilon)~=| --- ------------|(delta). | ker x | | \/2 | For very small x, the relative error amplification factor is 1 approximately given by ------, which implies a strong |logx| attenuation of relative error. However, (epsilon) in general cannot be less than the machine precision. For small x, errors are damped by the function and hence are limited by the machine precision. For medium and large x, the error behaviour, like the function itself, is oscillatory, and hence only the absolute accuracy for the function can be maintained. For this range of x, the / (pi)x -x/\/2 amplitude of the absolute error decays like / -----e \/ 2 which implies a strong attenuation of error. Eventually, ker x, / (pi) -x/\/2 which asymptotically behaves like / ----e , becomes so \/ 2x small that it cannot be calculated without causing underflow, and the routine returns zero. Note that for large x the errors are dominated by those of the Fortran intrinsic function EXP. 8. Further Comments Underflow may occur for a few values of x close to the zeros of ker x, below the limit which causes a failure with IFAIL = 1. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs19adf}{NAG On-line Documentation: s19adf} \beginscroll \begin{verbatim} S19ADF(3NAG) Foundation Library (12/10/92) S19ADF(3NAG) S19 -- Approximations of Special Functions S19ADF S19ADF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S19ADF returns a value for the Kelvin function keix via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S19ADF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the Kelvin function keix. Note: for x<0 the function is undefined, so we need only consider x>=0. The routine is based on several Chebyshev expansions: For 0<=x<=1, 2 (pi) x keix=- ----f(t)+ --[-g(t)logx+v(t)] 4 4 4 where f(t), g(t) and v(t) are expansions in the variable t=2x -1; For 13, / (pi) -x/\/2[( 1) 1 ] keix= / ----e [(1+ -)c(t)sin(beta)+ -d(t)cos(beta)] \/ 2x [( x) x ] x (pi) where (beta)= ---+ ----, and c(t) and d(t) are expansions in the 8 \/2 6 variable t= --1. x For x<0, the function is undefined, and hence the routine fails and returns zero. When x is sufficiently close to zero, the result is computed as 2 (pi) ( ( x)) x keix=- ----+(1-(gamma)-log( -)) -- 4 ( ( 2)) 4 and when x is even closer to zero simply as (pi) keix=- ----. 4 / (pi) -x/\/2 For large x, keix is asymptotically given by / ----e and \/ 2x this becomes so small that it cannot be computed without underflow and the routine fails. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. Constraint: X >= 0. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 On entry X is too large, the result underflows. On soft failure, the routine returns zero. IFAIL= 2 On entry X < 0, the function is undefined. On soft failure the routine returns zero. 7. Accuracy Let E be the absolute error in the result, and (delta) be the relative error in the argument. If (delta) is somewhat larger than the machine representation error, then we have: | x | E~=| ---(-ker x+kei x)|(delta). | 1 1 | | \/2 | For small x, errors are attenuated by the function and hence are limited by the machine precision. For medium and large x, the error behaviour, like the function itself, is oscillatory and hence only absolute accuracy of the function can be maintained. For this range of x, the amplitude of / (pi)x -x/\/2 the absolute error decays like / -----e , which implies a \/ 2 strong attenuation of error. Eventually, keix, which is / (pi) -x/\/2 asymptotically given by / ----e , becomes so small that it \/ 2x cannot be calculated without causing underflow and therefore the routine returns zero. Note that for large x, the errors are dominated by those of the Fortran intrinsic function EXP. 8. Further Comments Underflow may occur for a few values of x close to the zeros of keix, below the limit which causes a failure with IFAIL = 1. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs20acf}{NAG On-line Documentation: s20acf} \beginscroll \begin{verbatim} S20ACF(3NAG) Foundation Library (12/10/92) S20ACF(3NAG) S20 -- Approximations of Special Functions S20ACF S20ACF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S20ACF returns a value for the Fresnel Integral S(x), via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S20ACF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the Fresnel Integral x / ( (pi) 2) S(x)= |sin( ----t )dt. / ( 2 ) 0 Note: S(x)=-S(-x), so the approximation need only consider x>=0.0 The routine is based on three Chebyshev expansions: For 0 a T (t) , with t=2( -) -1; -- r r ( 3) r=0 For x>3, 1 f(x) ( (pi) 2) g(x) ( (pi) 2) S(x)= -- ----cos( ----x )- ----sin( ----x ), 2 x ( 2 ) 3 ( 2 ) x --' where f(x) = > b T (t), -- r r r=0 --' ( 3)4 and g(x) = > c T (t), with t=2( -) -1. -- r r ( x) r=0 (pi) 3 For small x, S(x)~= ----x . This approximation is used when x is 6 sufficiently small for the result to be correct to machine precision. For very small x, this approximation would underflow; the result is then set exactly to zero. 1 1 For large x, f(x)~= ---- and g(x)~= -----. Therefore for (pi) 2 (pi) 1 1 moderately large x, when ------- is negligible compared with -, 2 3 2 (pi) x the second term in the approximation for x>3 may be dropped. For 1 1 very large x, when ----- becomes negligible, S(x)~= -. However (pi)x 2 there will be considerable difficulties in calculating ( (pi) 2) cos( ----x ) accurately before this final limiting value can be ( 2 ) ( (pi) 2) used. Since cos( ----x ) is periodic, its value is essentially ( 2 ) 2 2 determined by the fractional part of x . If x =N+(theta) where N ( (pi) 2) is an integer and 0<=(theta)<1, then cos( ----x ) depends on ( 2 ) (theta) and on N modulo 4. By exploiting this fact, it is possible to retain significance in the calculation of ( (pi) 2) cos( ----x ) either all the way to the very large x limit, or at ( 2 ) x least until the integer part of - is equal to the maximum 2 integer allowed on the machine. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings There are no failure exits from this routine. The parameter IFAIL has been included for consistency with other routines in this chapter. 7. Accuracy Let (delta) and (epsilon) be the relative errors in the argument and result respectively. If (delta) is somewhat larger than the machine precision (i.e., if (delta) is due to data errors etc), then (epsilon) and (delta) are approximately related by: | ( (pi) 2)| | xsin( ----x )| | ( 2 )| (epsilon)~=| -------------|(delta). | S(x) | Figure 1 shows the behaviour of the error amplification factor | ( (pi) 2)| | xsin( ----x )| | ( 2 )| | -------------|. | S(x) | Figure 1 Please see figure in printed Reference Manual However if (delta) is of the same order as the machine precision, then rounding errors could make (epsilon) slightly larger than the above relation predicts. For small x, (epsilon)~=3(delta) and hence there is only moderate amplification of relative error. Of course for very small x where the correct result would underflow and exact zero is returned, relative error-control is lost. For moderately large values of x, | ( (pi) 2)| |(epsilon)|~=|2xsin( ----x )||(delta)| | ( 2 )| and the result will be subject to increasingly large amplification of errors. However the above relation breaks down 1 for large values of x (i.e., when -- is of the order of the 2 x machine precision in this region the relative error in the result 2 is essentially bounded by -----). (pi)x Hence the effects of error amplification are limited and at worst the relative error loss should not exceed half the possible number of significant figures. 8. Further Comments None. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs20adf}{NAG On-line Documentation: s20adf} \beginscroll \begin{verbatim} S20ADF(3NAG) Foundation Library (12/10/92) S20ADF(3NAG) S20 -- Approximations of Special Functions S20ADF S20ADF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S20ADF returns a value for the Fresnel Integral C(x), via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S20ADF (X, IFAIL) INTEGER IFAIL DOUBLE PRECISION X 3. Description This routine evaluates an approximation to the Fresnel Integral x / ( (pi) 2) C(x)= |cos( ----t )dt. / ( 2 ) 0 Note: C(x)=-C(-x), so the approximation need only consider x>=0.0 The routine is based on three Chebyshev expansions: For 0 a T (t) , with t=2( -) -1; -- r r ( 3) r=0 For x>3, 1 f(x) ( (pi) 2) g(x) ( (pi) 2) C(x)= -+ ----sin( ----x )- ----cos( ----x ), 2 x ( 2 ) 3 ( 2 ) x --' where f(x) = > b T (t), -- r r r=0 --' ( 3)4 and g(x) = > c T (t), with t=2( -) -1. -- r r ( x) r=0 For small x, C(x)~=x. This approximation is used when x is sufficiently small for the result to be correct to machine precision. 1 1 For large x, f(x)~= ---- and g(x)~= -----. Therefore for (pi) 2 (pi) 1 1 moderately large x, when ------- is negligible compared with -, 2 3 2 (pi) x the second term in the approximation for x>3 may be dropped. For 1 1 very large x, when ----- becomes negligible, C(x)~= -. However (pi)x 2 there will be considerable difficulties in calculating ( (pi) 2) sin( ----x ) accurately before this final limiting value can be ( 2 ) ( (pi) 2) used. Since sin( ----x ) is periodic, its value is essentially ( 2 ) 2 2 determined by the fractional part of x . If x =N+(theta), where N ( (pi) 2) is an integer and 0<=(theta)<1, then sin( ----x ) depends on ( 2 ) (theta) and on N modulo 4. By exploiting this fact, it is possible to retain some significance in the calculation of ( (pi) 2) sin( ----x ) either all the way to the very large x limit, or at ( 2 ) x least until the integer part of - is equal to the maximum 2 integer allowed on the machine. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. 5. Parameters 1: X -- DOUBLE PRECISION Input On entry: the argument x of the function. 2: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings There are no failure exits from this routine. The parameter IFAIL has been included for consistency with other routines in this chapter. 7. Accuracy Let (delta) and (epsilon) be the relative errors in the argument and result respectively. If (delta) is somewhat larger than the machine precision (i.e if (delta) is due to data errors etc), then (epsilon) and (delta) are approximately related by: | ( (pi) 2)| | xcos( ----x )| | ( 2 )| (epsilon)~=| -------------|(delta). | C(x) | Figure 1 shows the behaviour of the error amplification factor | ( (pi) 2)| | xcos( ----x )| | ( 2 )| | -------------|. | C(x) | Figure 1 Please see figure in printed Reference Manual However if (delta) is of the same order as the machine precision, then rounding errors could make (epsilon) slightly larger than the above relation predicts. For small x, (epsilon)~=(delta) and there is no amplification of relative error. For moderately large values of x, | ( (pi) 2)| |(epsilon)|~=|2xcos( ----x )||(delta)| | ( 2 )| and the result will be subject to increasingly large amplification of errors. However the above relation breaks down 1 for large values of x (i.e., when -- is of the order of the 2 x machine precision); in this region the relative error in the 2 result is essentially bounded by -----). (pi)x Hence the effects of error amplification are limited and at worst the relative error loss should not exceed half the possible number of significant figures. 8. Further Comments None. 9. Example The example program reads values of the argument x from a file, evaluates the function at each value of x and prints the results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs21baf}{NAG On-line Documentation: s21baf} \beginscroll \begin{verbatim} S21BAF(3NAG) Foundation Library (12/10/92) S21BAF(3NAG) S21 -- Approximations of Special Functions S21BAF S21BAF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S21BAF returns a value of an elementary integral, which occurs as a degenerate case of an elliptic integral of the first kind, via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S21BAF (X, Y, IFAIL) INTEGER IFAIL DOUBLE PRECISION X, Y 3. Description This routine calculates an approximate value for the integral infty 1 / dt R (x,y)= - | ---------- C 2 / 0 \/t+x(t+y) where x>=0 and y/=0. This function, which is related to the logarithm or inverse hyperbolic functions for y= 0.0 and Y /= 0.0. 3: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry X < 0.0; the function is undefined. IFAIL= 2 On entryY = 0.0; the function is undefined. On soft failure the routine returns zero. 7. Accuracy In principle the routine is capable of producing full machine precision. However round-off errors in internal arithmetic will result in slight loss of accuracy. This loss should never be excessive as the algorithm does not involve any significant amplification of round-off error. It is reasonable to assume that the result is accurate to within a small multiple of the machine precision. 8. Further Comments Users should consult the Chapter Introduction which shows the relationship of this function to the classical definitions of the elliptic integrals. 9. Example This example program simply generates a small set of non-extreme arguments which are used with the routine to produce the table of low accuracy results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs21bbf}{NAG On-line Documentation: s21bbf} \beginscroll \begin{verbatim} S21BBF(3NAG) Foundation Library (12/10/92) S21BBF(3NAG) S21 -- Approximations of Special Functions S21BBF S21BBF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S21BBF returns a value of the symmetrised elliptic integral of the first kind, via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S21BBF (X, Y, Z, IFAIL) INTEGER IFAIL DOUBLE PRECISION X, Y, Z 3. Description This routine calculates an approximation to the integral infty 1 / dt R (x,y,z)= - | ----------------- F 2 / 0 \/(t+x)(t+y)(t+z) where x, y, z>=0 and at most one is zero. The basic algorithm, which is due to Carlson [2] and [3], is to reduce the arguments recursively towards their mean by the rule: x =min(x,y,z),z =max(x,y,z), 0 0 y = remaining third intermediate value argument. 0 (This ordering, which is possible because of the symmetry of the function, is done for technical reasons related to the avoidance of overflow and underflow.) (mu) = (x +y +3z )/3 n n n n X = (1-x )/(mu) n n n Y = (1-y )/(mu) n n n Z = (1-z )/(mu) n n n (lambda) = /x y + /y z + /z x n \/ n n / n n / n n x = (x +(lambda) )/4 n+1 n n y = (y +(lambda) )/4 n+1 n n z = (z +(lambda) )/4 n+1 n n (epsilon) =max(|X |,|Y |,|Z |) and the function may be n n n n approximated adequately by a 5th order power series: ( 2 ) ( E E 3E E E ) 1 ( 2 2 2 3 3) R (x,y,z)= --------(1- --+ --- -----+ --) F ( 10 24 44 14) /(mu) \/ n where E =X Y +Y Z +Z X ,E =X Y Z . 2 n n n n n n 3 n n n The truncation error involved in using this approximation is 6 bounded by (epsilon) /4(1-(epsilon) ) and the recursive process n n is stopped when this truncation error is negligible compared with the machine precision. Within the domain of definition, the function value is itself representable for all representable values of its arguments. However, for values of the arguments near the extremes the above algorithm must be modified so as to avoid causing underflows or overflows in intermediate steps. In extreme regions arguments are pre-scaled away from the extremes and compensating scaling of the result is done before returning to the calling program. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. [2] Carlson B C (1978) Computing Elliptic Integrals by Duplication. (Preprint) Department of Physics, Iowa State University. [3] Carlson B C (1988) A Table of Elliptic Integrals of the Third Kind. Math. Comput. 51 267--280. 5. Parameters 1: X -- DOUBLE PRECISION Input 2: Y -- DOUBLE PRECISION Input 3: Z -- DOUBLE PRECISION Input On entry: the arguments x, y and z of the function. Constraint: X, Y, Z >= 0.0 and only one of X, Y and Z may be zero. 4: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry one or more of X, Y and Z is negative; the function is undefined. IFAIL= 2 On entry two or more of X, Y and Z are zero; the function is undefined. On soft failure, the routine returns zero. 7. Accuracy In principle the routine is capable of producing full machine precision. However round-off errors in internal arithmetic will result in slight loss of accuracy. This loss should never be excessive as the algorithm does not involve any significant amplification of round-off error. It is reasonable to assume that the result is accurate to within a small multiple of the machine precision. 8. Further Comments Users should consult the Chapter Introduction which shows the relationship of this function to the classical definitions of the elliptic integrals. If two arguments are equal, the function reduces to the elementary integral R , computed by S21BAF. C 9. Example This example program simply generates a small set of non-extreme arguments which are used with the routine to produce the table of low accuracy results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs21bcf}{NAG On-line Documentation: s21bcf} \beginscroll \begin{verbatim} S21BCF(3NAG) Foundation Library (12/10/92) S21BCF(3NAG) S21 -- Approximations of Special Functions S21BCF S21BCF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S21BCF returns a value of the symmetrised elliptic integral of the second kind, via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S21BCF (X, Y, Z, IFAIL) INTEGER IFAIL DOUBLE PRECISION X, Y, Z 3. Description This routine calculates an approximate value for the integral infty 3 / dt R (x,y,z)= - | ------------------- D 2 / 0 / 3 \/ (t+x)(t+y)(t+z) where x, y>=0, at most one of x and y is zero, and z>0. The basic algorithm, which is due to Carlson [2] and [3], is to reduce the arguments recursively towards their mean by the rule: x = x 0 y = y 0 z = z 0 (mu) = (x +y +3z )/5 n n n n X = (1-x )/(mu) n n n Y = (1-y )/(mu) n n n Z = (1-z )/(mu) n n n (lambda) = /x y + /y z + /z x n\/ n n / n n / n n x = (x +(lambda) )/4 n+1 n n y = (y +(lambda) )/4 n+1 n n z = (z +(lambda) )/4 n+1 n n For n sufficiently large, ( 1)n (epsilon) =max(|X |,|Y |,|Z |)~( -) n n n n ( 4) and the function may be approximated adequately by a 5th order n-1 -m -- 4 power series R (x,y,z)=3 > -------------------+ D -- m=0 (z +(lambda) ) /z m n \/ m -n 4 [ 3 (2) 1 (3) 3 (2) 2 3 (4) 3 (2) (3) 3 (5)] ---------[1+ -S + -S + --(S ) + --S + --S S + --S ] [ 7 n 3 n 22 n 11 n 13 n n 13 n ] / 3 / (mu) \/ n where (m) ( m m m) S =(X +Y +3Z )/2m. n ( n n n) The truncation error in this expansion is bounded by 6 3(epsilon) n ------------------- and the recursive process is terminated when / 3 / (1-(epsilon) ) \/ n this quantity is negligible compared with the machine precision. The routine may fail either because it has been called with arguments outside the domain of definition, or with arguments so extreme that there is an unavoidable danger of setting underflow or overflow. -3/2 Note: R (x,x,x)=x , so there exists a region of extreme D arguments for which the function value is not representable. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. [2] Carlson B C (1978) Computing Elliptic Integrals by Duplication. (Preprint) Department of Physics, Iowa State University. [3] Carlson B C (1988) A Table of Elliptic Integrals of the Third Kind. Math. Comput. 51 267--280. 5. Parameters 1: X -- DOUBLE PRECISION Input 2: Y -- DOUBLE PRECISION Input 3: Z -- DOUBLE PRECISION Input On entry: the arguments x, y and z of the function. Constraint: X, Y >= 0.0, Z > 0.0 and only one of X and Y may be zero. 4: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry either X or Y is negative, or both X and Y are zero; the function is undefined. IFAIL= 2 On entry Z <= 0.0; the function is undefined. IFAIL= 3 On entry either Z is too close to zero or both X and Y are too close to zero: there is a danger of setting overflow. IFAIL= 4 On entry at least one of X, Y and Z is too large: there is a danger of setting underflow. On soft failure the routine returns zero. 7. Accuracy In principle the routine is capable of producing full machine precision. However round-off errors in internal arithmetic will result in slight loss of accuracy. This loss should never be excessive as the algorithm does not involve any significant amplification of round-off error. It is reasonable to assume that the result is accurate to within a small multiple of the machine precision. 8. Further Comments Users should consult the Chapter Introduction which shows the relationship of this function to the classical definitions of the elliptic integrals. 9. Example This example program simply generates a small set of non-extreme arguments which are used with the routine to produce the table of low accuracy results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXs21bdf}{NAG On-line Documentation: s21bdf} \beginscroll \begin{verbatim} S21BDF(3NAG) Foundation Library (12/10/92) S21BDF(3NAG) S21 -- Approximations of Special Functions S21BDF S21BDF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose S21BDF returns a value of the symmetrised elliptic integral of the third kind, via the routine name. 2. Specification DOUBLE PRECISION FUNCTION S21BDF (X, Y, Z, R, IFAIL) INTEGER IFAIL DOUBLE PRECISION X, Y, Z, R 3. Description This routine calculates an approximation to the integral infty 3 / dt R (x,y,z,(rho))= - | -------------------------- J 2 / 0 (t+(rho))\/(t+x)(t+y)(t+z) where x, y, z>=0, (rho)/=0 and at most one of x, y and z is zero. If p<0, the result computed is the Cauchy principal value of the integral. The basic algorithm, which is due to Carlson [2] and [3], is to reduce the arguments recursively towards their mean by the rule: x = x 0 y = y 0 z = z 0 (rho) = (rho) 0 (mu) = (x +y +z +2(rho) )/5 n n n n n X = (1-x )/(mu) n n n Y = (1-y )/(mu) n n n Z = (1-z )/(mu) n n n P = (1-(rho) )/(mu) n n n (lambda) = /x y + /y z + /z x n\/ n n / n n / n n x = (x +(lambda) )/4 n+1 n n y = (y +(lambda) )/4 n+1 n n z = (z +(lambda) )/4 n+1 n n (rho) = ((rho) +(lambda) )/4 n+1 n n 2 [ ] (alpha) =[(rho) ( /x + /y + /z )+ /x y z ] n [ n \/ n / n / n / n n n] 2 (beta) = (rho) ((rho) +(lambda) ) n n n n For n sufficiently large, 1 (epsilon) =max(|X |,|Y |,|Z |,|P |)~ -- n n n n n n 4 and the function may be approximated by a 5th order power series n-1 -- -m R (x,y,z,(rho))=3 > 4 R ((alpha) ,(beta) ) J -- C m m m=0 -n 4 [ 3 (2) 1 (3) 3 (2) 2 3 (4) 3 (2) (3) 3 (5)] + ---------[1+ -S + -S + --(S ) + --S + --S S + --S ] [ 7 n 3 n 22 n 11 n 13 n n 13 n ] / 3 / (mu) \/ n (m) m m m m where S = (X +Y +Z +2P )/2m. n n n n n The truncation error in this expansion is bounded by 6 / 3 3(epsilon) / / (1-(epsilon) ) and the recursion process is n \/ n terminated when this quantity is negligible compared with the machine precision. The routine may fail either because it has been called with arguments outside the domain of definition or with arguments so extreme that there is an unavoidable danger of setting underflow or overflow. 3 - - 2 Note: R (x,x,x,x)=x , so there exists a region of extreme J arguments for which the function value is not representable. 4. References [1] Abramowitz M and Stegun I A (1968) Handbook of Mathematical Functions. Dover Publications. [2] Carlson B C (1978) Computing Elliptic Integrals by Duplication. (Preprint) Department of Physics, Iowa State University. [3] Carlson B C (1988) A Table of Elliptic Integrals of the Third Kind. Math. Comput. 51 267--280. 5. Parameters 1: X -- DOUBLE PRECISION Input 2: Y -- DOUBLE PRECISION Input 3: Z -- DOUBLE PRECISION Input 4: R -- DOUBLE PRECISION Input On entry: the arguments x, y, z and (rho) of the function. Constraint: X, Y, Z >= 0.0, R /= 0.0 and at most one of X, Y and Z may be zero. 5: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry at least one of X, Y and Z is negative, or at least two of them are zero; the function is undefined. IFAIL= 2 On entry R = 0.0; the function is undefined. IFAIL= 3 On entry either R is too close to zero, or any two of X, Y and Z are too close to zero; there is a danger of setting overflow. IFAIL= 4 On entry at least one of X, Y, Z and R is too large; there is a danger of setting underflow. IFAIL= 5 An error has occurred in a call to S21BAF. Any such occurrence should be reported to NAG. On soft failure, the routine returns zero. 7. Accuracy In principle the routine is capable of producing full machine precision. However round-off errors in internal arithmetic will result in slight loss of accuracy. This loss should never be excessive as the algorithm does not involve any significant amplification of round-off error. It is reasonable to assume that the result is accurate to within a small multiple of the machine precision. 8. Further Comments Users should consult the Chapter Introduction which shows the relationship of this function to the classical definitions of the elliptic integrals. If the argument R is equal to any of the other arguments, the function reduces to the integral R , computed by S21BCF. D 9. Example This example program simply generates a small set of non-extreme arguments which are used with the routine to produce the table of low accuracy results. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page}