path: root/src/hyper/pages/nags.ht

\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}