Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox Chapter IntroductionS — Approximations of Special Functions

## Scope of the Chapter

This chapter is concerned with the provision of some commonly occurring physical and mathematical functions.

## Background to the Problems

The majority of the functions in this chapter approximate real-valued functions of a single real argument, and the techniques involved are described in Section [Functions of a Single Real Argument]. In addition the chapter contains functions for elliptic integrals (see Section [Approximations to Elliptic Integrals]), Bessel and Airy functions of a complex argument (see Section [Bessel and Airy Functions of a Complex Argument]), complementary error function of a complex argument, hypergeometric functions and various option pricing functions for use in financial applications.

### Functions of a Single Real Argument

Most of the functions provided for functions of a single real argument have been based on truncated Chebyshev expansions. This method of approximation was adopted as a compromise between the conflicting requirements of efficiency and ease of implementation on many different machine ranges. For details of the reasons behind this choice and the production and testing procedures followed in constructing this chapter see Schonfelder (1976).
Basically, if the function to be approximated is f(x)$f\left(x\right)$, then for x[a,b]$x\in \left[a,b\right]$ an approximation of the form
 f(x) = g(x) ∑′ CrTr(t) r = 0
$f(x)=g(x)∑′r=0CrTr(t)$
is used (${\sum }^{\prime }$ denotes, according to the usual convention, a summation in which the first term is halved), where g(x)$g\left(x\right)$ 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)$t=t\left(x\right)$ is a mapping of the general range [a,b]$\left[a,b\right]$ to the specific range [1, + 1$-1,+1$] required by the Chebyshev polynomials, Tr(t)${T}_{r}\left(t\right)$. For a detailed description of the properties of the Chebyshev polynomials see Clenshaw (1962) and Fox and Parker (1968).
The essential property of these polynomials for the purposes of function approximation is that Tn(t)${T}_{n}\left(t\right)$ oscillates between ± 1$±1$ and it takes its extreme values n + 1$n+1$ times in the interval [1, + 1$-1,+1$]. Therefore, provided the coefficients Cr${C}_{r}$ decrease in magnitude sufficiently rapidly the error made by truncating the Chebyshev expansion after n$n$ terms is approximately given by
 E(t) ≃ CnTn(t). $E(t)≃CnTn(t).$
That is, the error oscillates between ± Cn$±{C}_{n}$ and takes its extreme value n + 1$n+1$ times in the interval in question. Now this is just the condition that the approximation be a minimax representation, one which minimizes the maximum error. By suitable choice of the interval, [a,b$a,b$], the auxiliary function, g(x)$g\left(x\right)$, and the mapping of the independent variable, t(x)$t\left(x\right)$, it is almost always possible to obtain a Chebyshev expansion with rapid convergence and hence truncations that provide near minimax polynomial approximations to the required function. The difference between the true minimax polynomial and the truncated Chebyshev expansion is seldom sufficiently great enough to be of significance.
The evaluation of the Chebyshev expansions follows one of two methods. The first and most efficient, and hence the most commonly used, works with the equivalent simple polynomial. The second method, which is used on the few occasions when the first method proves to be unstable, is based directly on the truncated Chebyshev series, and uses backward recursion to evaluate the sum. For the first method, a suitably truncated Chebyshev expansion (truncation is chosen so that the error is less than the machine precision) is converted to the equivalent simple polynomial. That is, we evaluate the set of coefficients br${b}_{r}$ such that
 n − 1 n − 1 y(t) = ∑ brtr = ∑′ CrTr(t). r = 0 r = 0
$y(t)=∑r=0 n-1brtr=∑′r=0 n-1CrTr(t).$
The polynomial can then be evaluated by the efficient Horner's method of nested multiplications,
 y(t) = (b0 + t(b1 + t(b2 + … t(bn − 2 + tbn − 1))) … ). $y(t)=(b0+t(b1+t(b2+…t(bn- 2+tbn- 1)))…).$
This method of evaluation results in efficient functions but for some expansions there is considerable loss of accuracy due to cancellation effects. In these cases the second method is used. It is well known that if
 bn − 1 = Cn − 1 bn − 2 = 2tbn − 1 + Cn − 2 bj − 0 = 2tbj + 1 − bj + 2 + Cj,  j = n − 3,n − 4, … ,0
$bn-1=Cn-1 bn-2=2tbn-1+Cn-2 bj-0=2tbj+1-bj+2+Cj, j=n-3,n-4,…,0$
then
 ∑′ CrTr(t) = (1/2)(b0 − b2) r = 0
$∑′r=0 CrTr(t)=12(b0-b2)$
and this is always stable. This method is most efficiently implemented by using three variables cyclically and explicitly constructing the recursion.
That is,
 α = Cn − 1 β = 2tα + Cn − 2 γ = 2tβ − α + Cn − 3 α = 2tγ − β + Cn − 4 β = 2tα − γ + Cn − 5 ⋮ say ​α = 2tγ − β + C2 β = 2tα − γ + C1 y(t) = tβ − α + (1/2)C0
$α = Cn-1 β = 2tα+Cn-2 γ = 2tβ-α+Cn-3 α = 2tγ-β+Cn-4 β = 2tα-γ+Cn-5 ⋮ say ​α = 2tγ-β+C2 β = 2tα-γ+C1 y(t) = tβ-α+12C0$
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$a,b$] to [1, + 1$-1,+1$] used range from simple linear mappings to the case when b$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 [a,a$-a,a$] can be provided using a mapping t = 2 (x / a)21$t=2{\left(x/a\right)}^{2}-1$. This embodies the evenness property but the expansion in t$t$ involves all powers and hence removes the necessity of working with an expansion with half its coefficients zero.
For many of the functions an analysis of the error in principle is given, namely, if E$E$ and $\nabla$ are the absolute errors in function and argument and ε$\epsilon$ and δ$\delta$ are the corresponding relative errors, then
 E ≃ |f′(x)|∇ E ≃ |xf′(x)|δ ε ≃ |( x f′ (x) )/(f(x))|δ.
$E ≃ |f′(x)|∇ E ≃ |xf′(x)|δ ε ≃ | x f′ (x) 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$x$. For example, on an 11$11$-digit machine
 |δ| ≤ 10 − 11. $|δ|≤10-11.$
(This of course implies that the absolute error = xδ$\nabla =x\delta$ is also bounded but the bound is now dependent on x$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)|$|x{f}^{\prime }\left(x\right)/f\left(x\right)|$ 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)|$|x{f}^{\prime }\left(x\right)|$ 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, of order ten or greater, the theoretical analysis gives a good measure of the accuracy obtainable.
It should be noted that the definitions and notations used for the functions in this chapter are all taken from Abramowitz and Stegun (1972). You are strongly recommended to consult this book for details before using the functions in this chapter.

### Approximations to Elliptic Integrals

Four functions provided here are symmetrised variants of the classical (Legendre) elliptic integrals. These alternative definitions have been suggested by Carlson (1965), Carlson (1977b) and Carlson (1977a) and he also developed the basic algorithms used in this chapter.
The symmetrised elliptic integral of the first kind is represented by
 ∞ RF(x,y,z) = (1/2) ∫ (dt)/( sqrt( (t + x) (t + y) (t + z) ) ), 0
$RF (x,y,z) = 12 ∫0∞ dt (t+x) (t+y) (t+z) ,$
where x,y,z0$x,y,z\ge 0$ and at most one may be equal to zero.
The normalization factor, (1/2) $\frac{1}{2}$, is chosen so as to make
 RF(x,x,x) = 1 / sqrt(x). $RF(x,x,x)=1/x.$
If any two of the variables are equal, RF${R}_{F}$ degenerates into the second function
 ∞ RC(x,y) = RF(x,y,y) = (1/2) ∫ (dt)/( (t + y) . sqrt(t + x) ), 0
$RC (x,y) = RF (x,y,y) = 12 ∫0∞ dt (t+y) . t+x ,$
where the argument restrictions are now x0$x\ge 0$ and y0$y\ne 0$.
This function is related to the logarithm or inverse hyperbolic functions if 0 < y < x$0, and to the inverse circular functions if 0xy$0\le x\le y$.
The symmetrised elliptic integral of the second kind is defined by
 ∞ RD(x,y,z) = (3/2) ∫ (dt)/( sqrt( (t + x) (t + y) (t + z)3 ) ) 0
$RD (x,y,z) = 32 ∫0∞ dt (t+x) (t+y) (t+z)3$
with z > 0$z>0$, x0$x\ge 0$ and y0$y\ge 0$, but only one of x$x$ or y$y$ may be zero.
The function is a degenerate special case of the symmetrised elliptic integral of the third kind
 ∞ RJ(x,y,z,ρ) = (3/2) ∫ (dt)/( sqrt( (t + x) (t + y) (t + z) ) (t + ρ) ) 0
$RJ (x,y,z,ρ) = 32 ∫0∞ dt (t+x) (t+y) (t+z) (t+ρ)$
with ρ0$\rho \ne 0$ and x,y,z0$x,y,z\ge 0$ with at most one equality holding. Thus RD(x,y,z) = RJ(x,y,z,z)${R}_{D}\left(x,y,z\right)={R}_{J}\left(x,y,z,z\right)$. The normalization of both these functions is chosen so that
 RD(x,x,x) = RJ(x,x,x,x) = 1 / (x×sqrt(x)). $RD(x,x,x)=RJ(x,x,x,x)=1/(x⁢x).$
The algorithms used for all these functions are based on duplication theorems. These allow a recursion system to be established which constructs a new set of arguments from the old using a combination of arithmetic and geometric means. The value of the function at the original arguments can then be simply related to the value at the new arguments. These recursive reductions are used until the arguments differ from the mean by an amount small enough for a Taylor series about the mean to give sufficient accuracy when retaining terms of order less than six. Each step of the recurrences reduces the difference from the mean by a factor of four, and as the truncation error is of order six, the truncation error goes like (4096)n${\left(4096\right)}^{-n}$, where n$n$ is the number of iterations.
The above forms can be related to the more traditional canonical forms (see Section 17.2 of Abramowitz and Stegun (1972)), as follows.
If we write q = cos2φ, r = 1msin2φ, s = 1nsin2φ $q={\mathrm{cos}}^{2}\varphi ,r=1-m {\mathrm{sin}}^{2}\varphi ,s=1-n {\mathrm{sin}}^{2}\varphi$, where 0φ(1/2)π $0\le \varphi \le \frac{1}{2}\pi$, we have
the classical elliptic integral of the first kind:
 φ F(φ ∣ m) = ∫ (1 − m   sin2θ) − (1/2)dθ = sinφ   RF(q,r,1); 0
$F(ϕ∣m) = ∫0ϕ ( 1-m sin2⁡θ ) -12 dθ = sin⁡ϕ RF (q,r,1) ;$
the classical elliptic integral of the second kind:
E(φm)
 φ = ∫ (1 − m   sin2θ)(1/2)dθ 0
= sinφ RF (q,r,1) (1/3)m sin3 φ RD (q,r,1)
$E(ϕ∣m) = ∫0ϕ ( 1-m sin2⁡θ ) 12 dθ = sin⁡ϕ RF (q,r,1) -13m sin3 ϕ RD (q,r,1)$
the classical elliptic integral of the third kind:
Π(n;φm)
 φ = ∫ (1 − n   sin2θ) − 1(1 − m   sin2θ) − (1/2)dθ 0
= sinφ RF (q,r,1) + (1/3) n sin3 φ RJ (q,r,1,s) .
$Π(n; ϕ∣m) = ∫0ϕ ( 1-n sin2⁡θ ) -1 ( 1-m sin2⁡θ ) -12 dθ = sin⁡ϕ RF (q,r,1) + 13 n sin3 ϕ RJ (q,r,1,s) .$
Also the classical complete elliptic integral of the first kind:
 π/2 K(m) = ∫ (1 − m   sin2θ) − (1/2)dθ = RF(0,1 − m,1); 0
$K(m) = ∫ 0 π2 ( 1 - m sin2⁡θ ) -12 dθ = RF (0,1-m,1) ;$
the classical complete elliptic integral of the second kind:
 π/2 E(m) = ∫ (1 − m   sin2θ)(1/2)dθ = RF(0,1 − m,1) − (1/3)m   RD(0,1 − m,1). 0
$E(m) = ∫ 0 π2 ( 1-m sin2 θ ) 12 dθ = RF (0,1-m,1) - 13 m RD (0,1-m,1) .$
For convenience, Chapter S contains functions to evaluate classical and symmetrised elliptic integrals.

### Bessel and Airy Functions of a Complex Argument

The functions for Bessel and Airy functions of a real argument are based on Chebyshev expansions, as described in Section [Functions of a Single Real Argument]. The functions provided for functions of a complex argument, however, use different methods. These functions relate all functions to the modified Bessel functions Iν(z)${I}_{\nu }\left(z\right)$ and Kν(z)${K}_{\nu }\left(z\right)$ computed in the right-half complex plane, including their analytic continuations. Iν${I}_{\nu }$ and Kν${K}_{\nu }$ are computed by different methods according to the values of z$z$ and ν$\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 (1972)).

### Option Pricing Functions

The option pricing functions evaluate the closed form solutions or approximations to the equations that define mathematical models for the prices of selected financial option contracts. These solutions can be viewed as special functions determined by the underlying equations. The terminology associated with these functions arises from their setting in financial markets and is briefly outlined below. See Joshi (2003) for a comprehensive introduction to this subject. An option is a contract which gives the holder the right, but not the obligation, to buy (if it is a call) or sell (if it is a put) a particular asset, S$S$. A European option can be exercised only at the specified expiry time, T$T$, while an American option can be exercised at any time up to T$T$. For Asian options the average underlying price over a pre-set time period determines the payoff.
The asset is bought (if a call) or sold (if a put) at a pre-specified strike price X$X$. Thus, an option contract has a payoff to the holder of max {(STX),0}$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left\{\left({S}_{T}-X\right),0\right\}$ for a call or max {(XST),0}$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left\{\left(X-{S}_{T}\right),0\right\}$, for a put, which depends on whether the asset price at the time of exercise is above (call) or below (put) the strike, X$X$. If at any moment in time a contract is currently showing a theoretical profit then it is deemed ‘in-the-money’; otherwise it is deemed ‘out-of-the-money’.
The option contract itself therefore has a value and, in many cases, can be traded in markets. Mathematical models, such as the Black–Scholes model, give theoretical prices for particular option contracts using a number of assumptions about the behaviour of financial markets. Typically, the price, St${S}_{t}$, of the underlying asset at time t$t$, is modelled as the solution of a stochastic differential equation for the return, dSt / St$d{S}_{t}/{S}_{t}$, on the asset price over a time interval, dt$dt$,
 (dSt)/(St) = μ dt + σdWt , $dSt St = μ dt + σdWt ,$
where dWt$d{W}_{t}$ is a Brownian motion. The drift, μ$\mu$, defines the trend in the movements of S$S$, while the volatility, σ$\sigma$, measures the risk and may be taken to be the standard deviation of the returns on the asset price. In addition the model requires a riskless money market account or bond with value, Bt${B}_{t}$, at time t$t$ and risk-free rate, r$r$, such that
 dBt = rBt dt . $dBt = rBt dt .$
This leads to the determination of the Black–Scholes option price, P$P$, by a martingale method or via the derivation of the Black–Scholes partial differential equation,
 (∂P)/(∂t) + (∂P)/(∂S) rS + (1/2) (∂2P)/(∂S2) σ2 S2 − rP = 0 . $∂P ∂t + ∂P ∂S rS+ 12 ∂2P ∂S2 σ2 S2 - rP=0 .$
For this case a closed form solution exists which is evaluated by nag_specfun_opt_bsm_price (s30aa).
A number of different option types where the solution exists in closed form or as a closed form approximation are presented in this chapter. See Haug (2007) for an extensive listing of option pricing formulae.

### Hypergeometric Functions

The confluent hypergeometric function M(a,b,x)$M\left(a,b,x\right)$ (or 1F1 (a ; b ; x) ${}_{1}F_{1}\left(a;b;x\right)$) requires a number of techniques to approximate it over the whole parameter (a,b)$\left(a,b\right)$ space and for all argument (x)$\left(x\right)$ values. For x$x$ well within the unit circle |x|ρ < 1$|x|\le \rho <1$ (where ρ = 0.8$\rho =0.8$ say), and for relatively small parameter values, the function can be well approximated by Taylor expansions, continued functions or through the solution of the related ordinary differential equation by an explicit, adaptive integrator. For values of |x| > ρ$|x|>\rho$, one of several transformations can be performed (depending on the value of x$x$) to reformulate the problem in terms of a new argument x${x}^{\prime }$ such that |x|ρ$|{x}^{\prime }|\le \rho$. If one or more of the parameters is relatively large (e.g., |a| > 30$|a|>30$) then recurrence relations can be used in combination to reformulate the problem in terms of parameter values of small size (e.g., |a| < 1$|a|<1$).
Approximations to the hypergeometric functions can therefore require all of the above techniques in sequence: a transformation to get an argument well inside the unit circle, a combination of recurrence relations to reduce the parameter sizes, and the approximation of the resulting hypergeometric function by one of a set of approximation techniques.
All the techniques described above are based on those described in Pearson (2009).

## Recommendations on Choice and Use of Available Functions

### Vectorized Function Variants

Many functions in Chapter S which compute functions of a single real argument have variants which operate on vectors of arguments. For example, nag_specfun_bessel_i0_real (s18ae) computes the value of the I0${I}_{0}$ Bessel function for a single argument, and nag_specfun_bessel_i0_real_vector (s18as) computes the same function for multiple arguments. In general it should be more efficient to use vectorized functions where possible, though to some extent this will depend on the environment from which you call the functions. See Section [Functionality Index] for a complete list of vectorized functions.

### 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 to the traditional canonical Legendre forms. In general, the extra symmetry of the Carlson forms is likely to simplify the analysis, and these symmetric forms are much more stable to calculate.
The function nag_specfun_ellipint_symm_1_degen (s21ba) for RC${R}_{C}$ is largely included as an auxiliary to the other functions for elliptic integrals. This integral essentially calculates elementary functions, e.g.,
 lnx = (x − 1)   RC (((1 + x)/2)2,x) ,  x > 0; arcsinx = x   RC(1 − x2,1),|x| ≤ 1; arcsinhx = x   RC(1 + x2,1),etc.
$ln⁡x =(x-1) RC ( (1+x2) 2,x) , x>0; arcsin⁡x =x RC(1-x2,1),|x|≤1; arcsinh⁡x =x RC(1+x2,1),etc.$
In general this method of calculating these elementary functions is not recommended as there are usually much more efficient specific functions available in the Library. However, nag_specfun_ellipint_symm_1_degen (s21ba) may be used, for example, to compute lnx / (x1)$\mathrm{ln}x/\left(x-1\right)$ when x$x$ is close to 1$1$, without the loss of significant figures that occurs when lnx$\mathrm{ln}x$ and x1$x-1$ are computed separately.

### Bessel and Airy Functions

For computing the Bessel functions Jν(x)${J}_{\nu }\left(x\right)$, Yν(x)${Y}_{\nu }\left(x\right)$, Iν(x)${I}_{\nu }\left(x\right)$ and Kν(x)${K}_{\nu }\left(x\right)$ where x$x$ is real and ν = 0​ or ​1$\nu =0\text{​ or ​}1$, special functions are provided, which are much faster than the more general functions that allow a complex argument and arbitrary real ν0$\nu \ge 0$. Similarly, special functions are provided for computing the Airy functions and their derivatives Ai(x)$\mathrm{Ai}\left(x\right)$, Bi(x)$\mathrm{Bi}\left(x\right)$, Ai(x)${\mathrm{Ai}}^{\prime }\left(x\right)$, Bi(x)${\mathrm{Bi}}^{\prime }\left(x\right)$ for a real argument which are much faster than the functions for complex arguments.

### Confluent Hypergeometric Function F 1 1

Two functions are provided for the confluent hypergeometric function 1F1 ${}_{1}F_{1}$. Both return values for 1F1 (a ; b ; x) ${}_{1}F_{1}\left(a;b;x\right)$ where parameters a$a$ and b$b$, and argument x$x$, are all real, but one variant works in a scaled form designed to avoid unnecessary loss of precision. The unscaled function nag_specfun_1f1_real (s22ba) is easier to use and should be chosen in the first instance, changing to the scaled function nag_specfun_1f1_real_scaled (s22bb) only if problems are encountered.

## Functionality Index

 Airy function,
 Ai, real argument,
 scalar nag_specfun_airy_ai_real (s17ag)
 vectorized nag_specfun_airy_ai_real_vector (s17au)
 Ai or Ai ′ , complex argument, optionally scaled nag_specfun_airy_ai_complex (s17dg)
 Ai ′ , real argument,
 scalar nag_specfun_airy_ai_deriv (s17aj)
 vectorized nag_specfun_airy_ai_deriv_vector (s17aw)
 Bi, real argument,
 scalar nag_specfun_airy_bi_real (s17ah)
 vectorized nag_specfun_airy_bi_real_vector (s17av)
 Bi or Bi ′ , complex argument, optionally scaled nag_specfun_airy_bi_complex (s17dh)
 Bi ′ , real argument,
 scalar nag_specfun_airy_bi_deriv (s17ak)
 vectorized nag_specfun_airy_bi_deriv_vector (s17ax)
 Arccos,
 inverse circular cosine nag_specfun_arccos (s09ab)
 Arccosh,
 inverse hyperbolic cosine nag_specfun_arccosh (s11ac)
 Arcsin,
 inverse circular sine nag_specfun_arcsin (s09aa)
 Arcsinh,
 inverse hyperbolic sine nag_specfun_arcsinh (s11ab)
 Arctanh,
 inverse hyperbolic tangent nag_specfun_arctanh (s11aa)
 Bessel function,
 I0, real argument,
 scalar nag_specfun_bessel_i0_real (s18ae)
 vectorized nag_specfun_bessel_i0_real_vector (s18as)
 I1, real argument,
 scalar nag_specfun_bessel_i1_real (s18af)
 vectorized nag_specfun_bessel_i1_real_vector (s18at)
 Iν, complex argument, optionally scaled nag_specfun_bessel_i_complex (s18de)
 J0, real argument,
 scalar nag_specfun_bessel_j0_real (s17ae)
 vectorized nag_specfun_bessel_j0_real_vector (s17as)
 J1, real argument,
 scalar nag_specfun_bessel_j1_real (s17af)
 vectorized nag_specfun_bessel_j1_real_vector (s17at)
 Jα ± n(z), complex argument nag_specfun_bessel_j_seq_complex (s18gk)
 Jν, complex argument, optionally scaled nag_specfun_bessel_j_complex (s17de)
 K0, real argument,
 scalar nag_specfun_bessel_k0_real (s18ac)
 vectorized nag_specfun_bessel_k0_real_vector (s18aq)
 K1, real argument,
 vectorized nag_specfun_bessel_k1_real_vector (s18ar)
 Kν, complex argument, optionally scaled nag_specfun_bessel_k_complex (s18dc)
 Y0, real argument,
 scalar nag_specfun_bessel_y0_real (s17ac)
 vectorized nag_specfun_bessel_y0_real_vector (s17aq)
 Y1, real argument,
 vectorized nag_specfun_bessel_y1_real_vector (s17ar)
 Yν, complex argument, optionally scaled nag_specfun_bessel_y_complex (s17dc)
 beta function,
 incomplete nag_specfun_beta_incomplete (s14cc)
 Complement of the Cumulative Normal distribution nag_specfun_compcdf_normal (s15ac)
 Complement of the Error function,
 complex argument, scaled nag_specfun_erfc_complex (s15dd)
 real argument, scaled nag_specfun_erfcx_real (s15ag)
 Cosine,
 hyperbolic nag_specfun_cosh (s10ac)
 Cosine Integral nag_specfun_integral_cos (s13ac)
 Cumulative Normal distribution function nag_specfun_cdf_normal (s15ab)
 Dawson's Integral nag_specfun_dawson (s15af)
 Elliptic functions, Jacobian, sn, cn, dn,
 complex argument nag_specfun_jacellip_complex (s21cb)
 real argument nag_specfun_jacellip_real (s21ca)
 Elliptic integral,
 general,
 of 2nd kind, F(z , k ′  , a , b) nag_specfun_ellipint_general_2 (s21da)
 Legendre form,
 complete of 1st kind, K(m) nag_specfun_ellipint_complete_1 (s21bh)
 complete of 2nd kind, E (m) nag_specfun_ellipint_complete_2 (s21bj)
 of 1st kind, F(ϕ | m) nag_specfun_ellipint_legendre_1 (s21be)
 of 2nd kind, E (ϕ ∣ m) nag_specfun_ellipint_legendre_2 (s21bf)
 of 3rd kind, Π (n ; ϕ ∣ m) nag_specfun_ellipint_legendre_3 (s21bg)
 symmetrised,
 degenerate of 1st kind, RC nag_specfun_ellipint_symm_1_degen (s21ba)
 of 1st kind, RF nag_specfun_ellipint_symm_1 (s21bb)
 of 2nd kind, RD nag_specfun_ellipint_symm_2 (s21bc)
 of 3rd kind, RJ nag_specfun_ellipint_symm_3 (s21bd)
 Erf,
 real argument nag_specfun_erf_real (s15ae)
 Erfc,
 complex argument, scaled nag_specfun_erfc_complex (s15dd)
 erfcx,
 real argument nag_specfun_erfcx_real (s15ag)
 Exponential,
 complex nag_specfun_exp_complex (s01ea)
 Exponential Integral nag_specfun_integral_exp (s13aa)
 Fresnel integral,
 C,
 vectorized nag_specfun_fresnel_c_vector (s20ar)
 S,
 scalar nag_specfun_fresnel_s (s20ac)
 vectorized nag_specfun_fresnel_s_vector (s20aq)
 Gamma function nag_specfun_gamma (s14aa)
 Gamma function,
 incomplete nag_specfun_gamma_incomplete (s14ba)
 Generalized factorial function nag_specfun_gamma (s14aa)
 Hankel function Hν(1) or Hν(2),
 complex argument, optionally scaled nag_specfun_hankel_complex (s17dl)
 Hypergeometric functions,
 1F1 (a ; b ; x) , confluent, real argument nag_specfun_1f1_real (s22ba)
 1F1(a ; b ; x), confluent, real argument, scaled form nag_specfun_1f1_real_scaled (s22bb)
 Jacobian theta functions θk(x , q),
 real argument nag_specfun_jactheta_real (s21cc)
 Kelvin function,
 bei x,
 scalar nag_specfun_kelvin_bei (s19ab)
 vectorized nag_specfun_kelvin_bei_vector (s19ap)
 ber x,
 scalar nag_specfun_kelvin_ber (s19aa)
 vectorized nag_specfun_kelvin_ber_vector (s19an)
 kei x,
 vectorized nag_specfun_kelvin_kei_vector (s19ar)
 ker x,
 scalar nag_specfun_kelvin_ker (s19ac)
 vectorized nag_specfun_kelvin_ker_vector (s19aq)
 Legendre functions of 1st kind Pnm(x), Pnm(x) nag_specfun_legendre_p (s22aa)
 Logarithm of 1 + x nag_specfun_log_shifted (s01ba)
 Logarithm of beta function,
 real nag_specfun_beta_log_real (s14cb)
 Logarithm of gamma function,
 complex nag_specfun_gamma_log_complex (s14ag)
 real nag_specfun_gamma_log_real (s14ab)
 real, scaled nag_specfun_gamma_log_scaled_real (s14ah)
 Option Pricing,
 American option, Bjerksund and Stensland option price nag_specfun_opt_amer_bs_price (s30qc)
 Asian option, geometric continuous average rate price nag_specfun_opt_asian_geom_price (s30sa)
 Asian option, geometric continuous average rate price with Greeks nag_specfun_opt_asian_geom_greeks (s30sb)
 binary asset-or-nothing option price nag_specfun_opt_binary_aon_price (s30cc)
 binary asset-or-nothing option price with Greeks nag_specfun_opt_binary_aon_greeks (s30cd)
 binary cash-or-nothing option price nag_specfun_opt_binary_con_price (s30ca)
 binary cash-or-nothing option price with Greeks nag_specfun_opt_binary_con_greeks (s30cb)
 Black–Scholes–Merton option price nag_specfun_opt_bsm_price (s30aa)
 Black–Scholes–Merton option price with Greeks nag_specfun_opt_bsm_greeks (s30ab)
 European option, option prices, using Merton jump-diffusion model nag_specfun_opt_jumpdiff_merton_price (s30ja)
 European option, option price with Greeks, using Merton jump-diffusion model nag_specfun_opt_jumpdiff_merton_greeks (s30jb)
 floating-strike lookback option price nag_specfun_opt_lookback_fls_price (s30ba)
 floating-strike lookback option price with Greeks nag_specfun_opt_lookback_fls_greeks (s30bb)
 Heston's model option price nag_specfun_opt_heston_price (s30na)
 Heston's model option price with Greeks nag_specfun_opt_heston_greeks (s30nb)
 standard barrier option price nag_specfun_opt_barrier_std_price (s30fa)
 Polygamma function,
 ψ(n)(x), real x nag_specfun_psi_deriv_real (s14ae)
 ψ(n)(z), complex z nag_specfun_psi_deriv_complex (s14af)
 Psi function nag_specfun_polygamma (s14ac)
 Psi function derivatives, scaled nag_specfun_polygamma_deriv (s14ad)
 Scaled modified Bessel function(s),
 e − (x)I0(x), real argument,
 scalar nag_specfun_bessel_i0_scaled (s18ce)
 vectorized nag_specfun_bessel_i0_scaled_vector (s18cs)
 e − (x)I1(x), real argument,
 scalar nag_specfun_bessel_i1_scaled (s18cf)
 vectorized nag_specfun_bessel_i1_scaled_vector (s18ct)
 ex K0 (x), real argument,
 scalar nag_specfun_bessel_k0_scaled (s18cc)
 vectorized nag_specfun_bessel_k0_scaled_vector (s18cq)
 ex K1 (x), real argument,
 scalar nag_specfun_bessel_k1_scaled (s18cd)
 vectorized nag_specfun_bessel_k1_scaled_vector (s18cr)
 Sine,
 hyperbolic nag_specfun_sinh (s10ab)
 Tangent,
 circular nag_specfun_tan (s07aa)
 hyperbolic nag_specfun_tanh (s10aa)
 Zeros of Bessel functions Jα(x), Jα ′ (x), Yα(x), Yα ′ (x),
 scalar nag_specfun_bessel_zeros (s17al)

## References

Abramowitz M and Stegun I A (1972) Handbook of Mathematical Functions (3rd Edition) Dover Publications
Carlson B C (1965) On computing elliptic integrals and functions J. Math. Phys. 44 36–51
Carlson B C (1977a) Special Functions of Applied Mathematics Academic Press
Carlson B C (1977b) Elliptic integrals of the first kind SIAM J. Math. Anal. 8 231–242
Clenshaw C W (1962) Chebyshev Series for Mathematical Functions Mathematical tables HMSO
Fox L and Parker I B (1968) Chebyshev Polynomials in Numerical Analysis Oxford University Press
Haug E G (2007) The Complete Guide to Option Pricing Formulas (2nd Edition) McGraw-Hill
Joshi M S (2003) The Concepts and Practice of Mathematical Finance Cambridge University Press
Pearson J (2009) Computation of hypergeometric functions MSc Dissertation, Mathematical Institute, University of Oxford
Schonfelder J L (1976) The production of special function routines for a multi-machine library Softw. Pract. Exper. 6(1)

Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013