diff options
| author | dos-reis <gdr@axiomatics.org> | 2007-08-14 05:14:52 +0000 |
|---|---|---|
| committer | dos-reis <gdr@axiomatics.org> | 2007-08-14 05:14:52 +0000 |
| commit | ab8cc85adde879fb963c94d15675783f2cf4b183 (patch) | |
| tree | c202482327f474583b750b2c45dedfc4e4312b1d /src/hyper/pages/nags.ht | |
| download | open-axiom-ab8cc85adde879fb963c94d15675783f2cf4b183.tar.gz | |
Initial population.
Diffstat (limited to 'src/hyper/pages/nags.ht')
| -rw-r--r-- | src/hyper/pages/nags.ht | 7689 |
1 files changed, 7689 insertions, 0 deletions
diff --git a/src/hyper/pages/nags.ht b/src/hyper/pages/nags.ht new file mode 100644 index 00000000..c6ee4a2c --- /dev/null +++ b/src/hyper/pages/nags.ht @@ -0,0 +1,7689 @@ +\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 0<y<x, and to the inverse circular functions if + 0<=x<=y. + + The integrals of the second kind are defined by + + infty + 3 / dt + R (x,y,z)= - | ------------------- + D 2 / + 0 / 3 + \/ (t+x)(t+y)(t+z) + + with z>0, 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<x<=4, the approximation is based on the Chebyshev expansion + + -- 1 + E (x)=y(t)-lnx= > '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<x<=16 it is based on the Chebyshev expansion + + --' ( x )2 + Ci(x)=lnx+ > a T (t) , t=2( --) -1. + -- r r ( 16) + r=0 + + For 16<x<x where the value of x is given in the Users' Note + hi hi + for your implementation, + + f(x)sinx g(x)cosx + Ci(x)= --------- -------- + x 2 + x + + --' --' ( 16)2 + where f(x)= > 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|<x , where x is an implementation dependent number, + hi hi + + { (pi) f(x)cosx g(x)sinx} + Si(x)=sign(x){ ----- --------- --------} + { 2 x 2 } + { x } + + --' --' ( 16)2 + where f(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<x<=x , ln(Gamma)(x)=-lnx to within machine accuracy. + small + + For x <x<=15.0, the recursive relation + small + (Gamma)(1+x)=x(Gamma)(x) is used to reduce the calculation to one + involving (Gamma)(1+u), 0<=u<1 which is evaluated as: + + --' + (Gamma)(1+u)= > 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<x<=x , + big + + 1 1 + ln(Gamma)(x)=(x- -)lnx-x+ -ln2(pi)+y(x)/x + 2 2 + 2 + --' ( 15) + where y(x)= > b T (t), t=2( --) -1. + -- r r ( x ) + r=0 + + For x <x<=x the term y(x)/x is negligible and so its + big vbig + calculation is omitted. + + 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 x<x , (Gamma)(x)=1/x to + small small + within machine accuracy. x is calculated so that if x>x , + 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<x <0, the result is returned as 2.0 which is correct to + low + within machine precision. The values of x and x are given in + hi low + 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 + + If (delta) and (epsilon) are relative errors in the argument and + result, respectively, then in principle + + | 2 | + | -x | + | 2xe | + |(epsilon)|~=| ------------(delta)|. + | | + + | \/(pi)erfc x | + + That is, the relative error in the argument, x, is amplified by a + 2 + -x + 2xe + factor ------------ in the result. + + + \/(pi)erfc x + + The behaviour of this factor is shown in Figure 1. + + + Figure 1 + + Please see figure in printed Reference Manual + + 2x + It should be noted that near x=0 this factor behaves as ------ + + + \/(pi) + and hence the accuracy is largely determined by the machine + precision. Also for large negative x, where the factor is + 2 + -x + xe + ~ ------, accuracy is mainly limited by machine precision. + + + \/(pi) + 2 + However, for large positive x, the factor becomes ~2x and to an + extent relative accuracy is necessarily lost. The absolute + accuracy E is given by + + 2 + -x + 2xe + E~= ------(delta) + + + \/(pi) + + so absolute accuracy is guaranteed for all x. + + 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}{manpageXXs15aef}{NAG On-line Documentation: s15aef} +\beginscroll +\begin{verbatim} + + + + S15AEF(3NAG) Foundation Library (12/10/92) S15AEF(3NAG) + + + + S15 -- Approximations of Special Functions S15AEF + S15AEF -- 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 + + S15AEF returns the value of the error function erfx, via the + routine name. + + 2. Specification + + DOUBLE PRECISION FUNCTION S15AEF (X, IFAIL) + INTEGER IFAIL + DOUBLE PRECISION X + + 3. Description + + Evaluates the error function, + + 2 + x -t + 2 / + erf x = ------ |e dt. + / + \/(pi) 0 + + For |x|<=2, + --' 1 2 + erf x = x > a T (t), where t= -x -1. + -- r r 2 + r=0 + + For 2<|x|<x , + hi + { 2 } + { -x } + { e -- } x-7 + erf x = signx{1- --------- > '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<x<=8, + 2 + 2 --' --' ( x) + Y (x)= ----lnx > 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<x<=8, + + 2 x --' 2 x --' + Y (x)= ----lnx - > 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<x<=8, 2 + --' ( x) + J (x)= > 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<x<=8, 2 + x --' ( x) + J (x)= - > 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<x<4.5, + + -3x/2 + Ai(x)=e y(t), + + where y is an expansion in t=4x/9-1. + + For 4.5<=x<9, + + -5x/2 + Ai(x)=e u(t), + + where u is an expansion in t=4x/9-3. + + For x>=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<x<4.5, + + 11x/8 + Bi(x)=e y(t), + + where y is an expansion in t=4x/9-1. + + For 4.5<=x<=9, + + 5x/2 + Bi(x)=e v(t), + + where v is an expansion in t=4x/9-3. + + For x>=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<x<4.5, + + -11x/8 + Ai'(x)=e y(t), + + ( x) + where y(t) is an expansion in t=4( -)-1. + ( 9) + + For 4.5<=x<9, + + -5x/2 + Ai'(x)=e v(t), + + ( x) + where v(t) is an expansion in t=4( -)-3. + ( 9) + + For x>=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<x<4.5, + + 3x/2 + Bi'(x)=e y(t), + + where y(t) is an expansion in t=4x/9-1. + + For 4.5<=x<9, + + 21x/8 + Bi'(x)=e u(t), + + where u(t) is an expansion in t=4x/9-3. + + For x>=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<x<=1, + + --' --' 2 + K (x)=-lnx > a T (t)+ > b T (t) , where t=2x -1; + 0 -- r r -- r r + r=0 r=0 + + For 1<x<=2, + + -x --' + K (x)=e > c T (t) , where t=2x-3; + 0 -- r r + r=0 + + For 2<x<=4, + + -x --' + K (x)=e > 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<x<=1, + + 1 --' --' 2 + K (x)= -+xlnx > a T (t)-x > b T (t) , where t=2x -1; + 1 x -- r r -- r r + r=0 r=0 + + For 1<x<=2, + + -x --' + K (x)=e > c T (t) , where t=2x-3; + 1 -- r r + r=0 + + For 2<x<=4, + + -x --' + K (x)=e > 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<x<=4, + + x --' ( x) + I (x)=e > a T (t) , where t=2( -)-1; + 0 -- r r ( 4) + r=0 + + For 4<x<=12, + + x --' x-8 + I (x)=e > 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<x<=4, + + --' ( x)2 + I (x)=x > a T (t) , where t=2( -) -1; + 1 -- r r ( 4) + r=0 + + For 4<x<=12, + + x --' x-8 + I (x)=e > 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 0<x<=1, + + (pi) 2 + ker x=-f(t)logx+ ----x g(t)+y(t) + 16 + + 4 + where f(t), g(t) and y(t) are expansions in the variable t=2x -1; + + For 1<x<=3, + + ( 11 ) + ker x=exp(- --x)q(t) + ( 16 ) + + where q(t) is an expansion in the variable t=x-2; + + For x>3, + + + + -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 1<x<=3, + + ( 9 ) + keix=exp(- -x)u(t) + ( 8 ) + + where u(t) is an expansion in the variable t=x-2; + + For x>3, + + + + / (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<x<=3, + + 3 --' ( x)4 + S(x)=x > 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<x<=3, + + --' ( x)4 + C(x)=x > 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<x and to inverse circular functions if + x<y, arises as a degenerate form of the elliptic integral of the + first kind. If y<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 + system: + + x =x, + 0 + + y =y + 0 + + (mu) =(x +2y )/3, + n n n + + S =(y -x )/3(mu) + n n n n + + + + (lambda) =y +2 /x y + n n \/ n n + + x =(x +(lambda) )/4, + n+1 n n + + y =(y +(lambda) )/4. + n+1 n n + + The quantity |S | for n=0,1,2,3,... decreases with increasing n, + n + n + eventually |S |~1/4 . For small enough S the required function + n n + value can be approximated by the first few terms of the Taylor + series about the mean. That is + + ( 2 3 4 5) + ( 3S S 3S 9S ) + ( n n n n) + R (x,y)=(1+ ---+ --+ ---+ ---)/ /(mu) . + C ( 10 7 8 22 ) \/ n + + The truncation error involved in using this approximation is + 6 + bounded by 16|S | /(1-2|S |) and the recursive process is stopped + n n + when S is small enough for this truncation error to be + n + negligible compared to 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 + On entry: the arguments x and y of the function, + respectively. Constraint: X >= 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} |