- :
`a`=**airy**`(`

¶`z`) - :
`a`=**airy**`(`

¶`k`,`z`) - :
`a`=**airy**`(`

¶`k`,`z`,`scale`) - :
`[`

`a`,`ierr`] =**airy**`(…)`

¶ -
Compute Airy functions of the first and second kind, and their derivatives.

K Function Scale factor (if

`scale`is true) --- -------- --------------------------------------- 0 Ai (Z) exp ((2/3) * Z * sqrt (Z)) 1 dAi(Z)/dZ exp ((2/3) * Z * sqrt (Z)) 2 Bi (Z) exp (-abs (real ((2/3) * Z * sqrt (Z)))) 3 dBi(Z)/dZ exp (-abs (real ((2/3) * Z * sqrt (Z))))The function call

`airy (`

is equivalent to`z`)`airy (0,`

.`z`)The optional third input

`scale`determines whether to apply scaling as described above. It is false by default.The result

`a`is the same size as`z`.The optional output

`ierr`contains the following status information and is the same size as the result.- Normal return.
- Input error, return
`NaN`

. - Overflow, return
`Inf`

. - Loss of significance by argument reduction results in less than half of machine accuracy.
- Loss of significance by argument reduction, output may be inaccurate.
- Error—no computation, algorithm termination condition not met,
return
`NaN`

.

- :
`J`=**besselj**`(`

¶`alpha`,`x`) - :
`J`=**besselj**`(`

¶`alpha`,`x`,`opt`) - :
`[`

`J`,`ierr`] =**besselj**`(…)`

¶ Compute Bessel functions of the first kind.

The order of the Bessel function

`alpha`must be real. The points for evaluation`x`may be complex.If the optional argument

`opt`is 1 or true, the result`J`is multiplied by`exp (-abs (imag (`

.`x`)))If

`alpha`is a scalar, the result is the same size as`x`. If`x`is a scalar, the result is the same size as`alpha`. If`alpha`is a row vector and`x`is a column vector, the result is a matrix with`length (`

rows and`x`)`length (`

columns. Otherwise,`alpha`)`alpha`and`x`must conform and the result will be the same size.If requested,

`ierr`contains the following status information and is the same size as the result.- Normal return.
- Input error, return
`NaN`

. - Overflow, return
`Inf`

. - Loss of significance by argument reduction results in less than half of machine accuracy.
- Loss of significance by argument reduction, output may be inaccurate.
- Error—no computation, algorithm termination condition not met, return
`NaN`

.

- :
`Y`=**bessely**`(`

¶`alpha`,`x`) - :
`Y`=**bessely**`(`

¶`alpha`,`x`,`opt`) - :
`[`

`Y`,`ierr`] =**bessely**`(…)`

¶ Compute Bessel functions of the second kind.

The order of the Bessel function

`alpha`must be real. The points for evaluation`x`may be complex.If the optional argument

`opt`is 1 or true, the result`Y`is multiplied by`exp (-abs (imag (`

.`x`)))If

`alpha`is a scalar, the result is the same size as`x`. If`x`is a scalar, the result is the same size as`alpha`. If`alpha`is a row vector and`x`is a column vector, the result is a matrix with`length (`

rows and`x`)`length (`

columns. Otherwise,`alpha`)`alpha`and`x`must conform and the result will be the same size.If requested,

`ierr`contains the following status information and is the same size as the result.- Normal return.
- Input error, return
`NaN`

. - Overflow, return
`Inf`

. - Loss of significance by argument reduction results in less than half of machine accuracy.
- Complete loss of significance by argument reduction, return
`NaN`

. - Error—no computation, algorithm termination condition not met, return
`NaN`

.

- :
`I`=**besseli**`(`

¶`alpha`,`x`) - :
`I`=**besseli**`(`

¶`alpha`,`x`,`opt`) - :
`[`

`I`,`ierr`] =**besseli**`(…)`

¶ Compute modified Bessel functions of the first kind.

The order of the Bessel function

`alpha`must be real. The points for evaluation`x`may be complex.If the optional argument

`opt`is 1 or true, the result`I`is multiplied by`exp (-abs (real (`

.`x`)))If

`alpha`is a scalar, the result is the same size as`x`. If`x`is a scalar, the result is the same size as`alpha`. If`alpha`is a row vector and`x`is a column vector, the result is a matrix with`length (`

rows and`x`)`length (`

columns. Otherwise,`alpha`)`alpha`and`x`must conform and the result will be the same size.If requested,

`ierr`contains the following status information and is the same size as the result.- Normal return.
- Input error, return
`NaN`

. - Overflow, return
`Inf`

. - Loss of significance by argument reduction results in less than half of machine accuracy.
- Complete loss of significance by argument reduction, return
`NaN`

. - Error—no computation, algorithm termination condition not met, return
`NaN`

.

- :
`K`=**besselk**`(`

¶`alpha`,`x`) - :
`K`=**besselk**`(`

¶`alpha`,`x`,`opt`) - :
`[`

`K`,`ierr`] =**besselk**`(…)`

¶ -
Compute modified Bessel functions of the second kind.

The order of the Bessel function

`alpha`must be real. The points for evaluation`x`may be complex.If the optional argument

`opt`is 1 or true, the result`K`is multiplied by`exp (`

.`x`)`alpha`is a scalar, the result is the same size as`x`. If`x`is a scalar, the result is the same size as`alpha`. If`alpha`is a row vector and`x`is a column vector, the result is a matrix with`length (`

rows and`x`)`length (`

columns. Otherwise,`alpha`)`alpha`and`x`must conform and the result will be the same size.If requested,

`ierr`contains the following status information and is the same size as the result.- Normal return.
- Input error, return
`NaN`

. - Overflow, return
`Inf`

. - Loss of significance by argument reduction results in less than half of machine accuracy.
- Complete loss of significance by argument reduction, return
`NaN`

. - Error—no computation, algorithm termination condition not met, return
`NaN`

.

- :
`H`=**besselh**`(`

¶`alpha`,`x`) - :
`H`=**besselh**`(`

¶`alpha`,`k`,`x`) - :
`H`=**besselh**`(`

¶`alpha`,`k`,`x`,`opt`) - :
`[`

`H`,`ierr`] =**besselh**`(…)`

¶ Compute Bessel functions of the third kind (Hankel functions).

The order of the Bessel function

`alpha`must be real. The kind of Hankel function is specified by`k`and may be either first (`k`= 1) or second (`k`= 2). The default is Hankel functions of the first kind. The points for evaluation`x`may be complex.If the optional argument

`opt`is 1 or true, the result is multiplied by`exp (-I*`

for`x`)`k`= 1 or`exp (I*`

for`x`)`k`= 2.`alpha`is a scalar, the result is the same size as`x`. If`x`is a scalar, the result is the same size as`alpha`. If`alpha`is a row vector and`x`is a column vector, the result is a matrix with`length (`

rows and`x`)`length (`

columns. Otherwise,`alpha`)`alpha`and`x`must conform and the result will be the same size.If requested,

`ierr`contains the following status information and is the same size as the result.- Normal return.
- Input error, return
`NaN`

. - Overflow, return
`Inf`

. - Loss of significance by argument reduction results in less than half of machine accuracy.
- Complete loss of significance by argument reduction, return
`NaN`

. - Error—no computation, algorithm termination condition not met, return
`NaN`

.

- :
`y`=**beta**`(`

¶`a`,`b`) Compute the Beta function for real inputs

`a`and`b`.The Beta function definition is

beta (a, b) = gamma (a) * gamma (b) / gamma (a + b).

The Beta function can grow quite large and it is often more useful to work with the logarithm of the output rather than the function directly. See

`betaln`

, for computing the logarithm of the Beta function in an efficient manner.**See also:**betaln, betainc, betaincinv.

- :
`I`=**betainc**`(`

¶`x`,`a`,`b`) - :
`I`=**betainc**`(`

¶`x`,`a`,`b`,`tail`) Compute the incomplete beta function.

This is defined as

x / 1 | I_x (a, b) = ---------- | t^(a-1) (1-t)^(b-1) dt beta (a,b) | / 0

with real

`x`in the range [0,1]. The inputs`a`and`b`must be real and strictly positive (> 0). If one of the inputs is not a scalar then the other inputs must be scalar or of compatible dimensions.By default,

`tail`is`"lower"`

and the incomplete beta function integrated from 0 to`x`is computed. If`tail`is`"upper"`

then the complementary function integrated from`x`to 1 is calculated. The two choices are related bybetainc (

`x`,`a`,`b`,`"upper"`

) = 1 - betainc (`x`,`a`,`b`,`"lower"`

).`betainc`

uses a more sophisticated algorithm than subtraction to get numerically accurate results when the`"lower"`

value is small.Reference: A. Cuyt, V. Brevik Petersen, B. Verdonk, H. Waadeland, W.B. Jones, Handbook of Continued Fractions for Special Functions, ch. 18.

**See also:**beta, betaincinv, betaln.

- :
`x`=**betaincinv**`(`

¶`y`,`a`,`b`) - :
`x`=**betaincinv**`(`

¶`y`,`a`,`b`, "lower") - :
`x`=**betaincinv**`(`

¶`y`,`a`,`b`, "upper") Compute the inverse of the normalized incomplete beta function.

The normalized incomplete beta function is defined as

x / 1 | I_x (a, b) = ---------- | t^(a-1) (1-t)^(b-1) dt beta (a,b) | / 0

If two inputs are scalar, then

`betaincinv (`

is returned for each of the other inputs.`y`,`a`,`b`)If two or more inputs are not scalar, the sizes of them must agree, and

`betaincinv`

is applied element-by-element.The variable

`y`must be in the interval [0,1], while`a`and`b`must be real and strictly positive.By default,

`tail`is`"lower"`

and the inverse of the incomplete beta function integrated from 0 to`x`is computed. If`tail`is`"upper"`

then the complementary function integrated from`x`to 1 is inverted.The function is computed by standard Newton’s method, by solving

`y`- betainc (`x`,`a`,`b`) = 0

- :
`lnb`=**betaln**`(`

¶`a`,`b`) Compute the natural logarithm of the Beta function for real inputs

`a`and`b`.`betaln`

is defined asbetaln (a, b) = log (beta (a, b))

and is calculated in a way to reduce the occurrence of underflow.

The Beta function can grow quite large and it is often more useful to work with the logarithm of the output rather than the function directly.

**See also:**beta, betainc, betaincinv, gammaln.

- :
`b`=**bincoeff**`(`

¶`n`,`k`) Return the binomial coefficient of

`n`and`k`.The binomial coefficient is defined as

/ \ | n | n (n-1) (n-2) ... (n-k+1) | | = ------------------------- | k | k! \ /

For example:

bincoeff (5, 2) ⇒ 10

In most cases, the

`nchoosek`

function is faster for small scalar integer arguments. It also warns about loss of precision for big arguments.**See also:**nchoosek.

- :
`k`=**commutation_matrix**`(`

¶`m`,`n`) Return the commutation matrix K(m,n) which is the unique

`m`*`n`by`m`*`n`matrix such that*K(m,n) * vec(A) = vec(A')*for all*m*by*n*matrices*A*.If only one argument

`m`is given,*K(m,m)*is returned.See Magnus and Neudecker (1988), Matrix Differential Calculus with Applications in Statistics and Econometrics.

- :
`y`=**cosint**`(`

¶`x`) Compute the cosine integral function:

+oo / Ci (x) = - | (cos (t)) / t dt / x

An equivalent definition is

x / | cos (t) - 1 Ci (x) = gamma + log (x) + | ------------- dt | t / 0

Reference:

M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, 1964.

- :
`d`=**duplication_matrix**`(`

¶`n`) Return the duplication matrix

*Dn*which is the unique*N^2*-by-*N*(N+1)/2*matrix such that`Dn * vech (A) = vec (A)`

for all symmetric*N*-by-*N*matrices*A*.See Magnus and Neudecker (1988), Matrix Differential Calculus with Applications in Statistics and Econometrics.

- :
`v`=**dawson**`(`

¶`z`) Compute the Dawson (scaled imaginary error) function.

The Dawson function is defined as

(sqrt (pi) / 2) * exp (-z^2) * erfi (z)

- :
`[`

`sn`,`cn`,`dn`,`err`] =**ellipj**`(`

¶`u`,`m`) - :
`[`

`sn`,`cn`,`dn`,`err`] =**ellipj**`(`

¶`u`,`m`,`tol`) Compute the Jacobi elliptic functions

`sn`,`cn`, and`dn`of complex argument`u`and real parameter`m`.If

`m`is a scalar, the results are the same size as`u`. If`u`is a scalar, the results are the same size as`m`. If`u`is a column vector and`m`is a row vector, the results are matrices with`length (`

rows and`u`)`length (`

columns. Otherwise,`m`)`u`and`m`must conform in size and the results will be the same size as the inputs.The value of

`u`may be complex. The value of`m`must be 0 ≤`m`≤ 1.The optional input

`tol`is currently ignored (MATLAB uses this to allow faster, less accurate approximation).If requested,

`err`contains the following status information and is the same size as the result.- Normal return.
- Error—no computation, algorithm termination condition not met,
return
`NaN`

.

Reference: Milton Abramowitz and Irene A Stegun, Handbook of Mathematical Functions, Chapter 16 (Sections 16.4, 16.13, and 16.15), Dover, 1965.

**See also:**ellipke.

- :
`k`=**ellipke**`(`

¶`m`) - :
`k`=**ellipke**`(`

¶`m`,`tol`) - :
`[`

`k`,`e`] =**ellipke**`(…)`

¶ Compute complete elliptic integrals of the first K(

`m`) and second E(`m`) kind.`m`must be a scalar or real array with -Inf ≤`m`≤ 1.The optional input

`tol`controls the stopping tolerance of the algorithm and defaults to`eps (class (`

. The tolerance can be increased to compute a faster, less accurate approximation.`m`))When called with one output only elliptic integrals of the first kind are returned.

Mathematical Note:

Elliptic integrals of the first kind are defined as

1 / dt K (m) = | ------------------------------ / sqrt ((1 - t^2)*(1 - m*t^2)) 0

Elliptic integrals of the second kind are defined as

1 / sqrt (1 - m*t^2) E (m) = | ------------------ dt / sqrt (1 - t^2) 0

Reference: Milton Abramowitz and Irene A. Stegun, Handbook of Mathematical Functions, Chapter 17, Dover, 1965.

**See also:**ellipj.

- :
`v`=**erf**`(`

¶`z`) Compute the error function.

The error function is defined as

z 2 / erf (z) = --------- * | e^(-t^2) dt sqrt (pi) / t=0

- :
`v`=**erfc**`(`

¶`z`) Compute the complementary error function.

The complementary error function is defined as

`1 - erf (`

.`z`)

- :
`v`=**erfcx**`(`

¶`z`) Compute the scaled complementary error function.

The scaled complementary error function is defined as

exp (z^2) * erfc (z)

- :
`v`=**erfi**`(`

¶`z`) Compute the imaginary error function.

The imaginary error function is defined as

-i * erf (i*z)

- :
`y`=**erfinv**`(`

¶`x`) Compute the inverse error function.

The inverse error function is defined such that

erf (

`y`) ==`x`

- :
`y`=**erfcinv**`(`

¶`x`) Compute the inverse complementary error function.

The inverse complementary error function is defined such that

erfc (

`y`) ==`x`

- :
`y`=**expint**`(`

¶`x`) Compute the exponential integral.

The exponential integral is defined as:

+oo / | exp (-t) E_1 (x) = | -------- dt | t / x

Note: For compatibility, this function uses the MATLAB definition of the exponential integral. Most other sources refer to this particular value as

*E_1 (x)*, and the exponential integral as+oo / | exp (-t) Ei (x) = - | -------- dt | t / -x

The two definitions are related, for positive real values of

`x`, by`E_1 (-x) = -Ei (x) - i*pi`

.References:

M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, 1964.

N. Bleistein and R.A. Handelsman, Asymptotic expansions of integrals, 1986.

- :
`v`=**gamma**`(`

¶`z`) Compute the Gamma function.

The Gamma function is defined as

infinity / gamma (z) = | t^(z-1) exp (-t) dt. / t=0

Programming Note: The gamma function can grow quite large even for small input values. In many cases it may be preferable to use the natural logarithm of the gamma function (

`gammaln`

) in calculations to minimize loss of precision. The final result is then`exp (`

`result_using_gammaln`).

- :
`y`=**gammainc**`(`

¶`x`,`a`) - :
`y`=**gammainc**`(`

¶`x`,`a`,`tail`) Compute the normalized incomplete gamma function.

This is defined as

x 1 / gammainc (x, a) = --------- | exp (-t) t^(a-1) dt gamma (a) / t=0

with the limiting value of 1 as

`x`approaches infinity. The standard notation is*P(a,x)*, e.g., Abramowitz and Stegun (6.5.1).If

`a`is scalar, then`gammainc (`

is returned for each element of`x`,`a`)`x`and vice versa.If neither

`x`nor`a`is scalar then the sizes of`x`and`a`must agree, and`gammainc`

is applied element-by-element. The elements of`a`must be non-negative.By default,

`tail`is`"lower"`

and the incomplete gamma function integrated from 0 to`x`is computed. If`tail`is`"upper"`

then the complementary function integrated from`x`to infinity is calculated.If

`tail`is`"scaledlower"`

, then the lower incomplete gamma function is multiplied by*gamma(a+1)*exp(x)/(x^a)*. If`tail`is`"scaledupper"`

, then the upper incomplete gamma function is multiplied by the same quantity.References:

M. Abramowitz and I.A. Stegun, Handbook of mathematical functions, Dover publications, Inc., 1972.

W. Gautschi, A computational procedure for incomplete gamma functions, ACM Trans. Math Software, pp. 466–481, Vol 5, No. 4, 2012.

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in Fortran 77, ch. 6.2, Vol 1, 1992.

**See also:**gamma, gammaincinv, gammaln.

- :
`x`=**gammaincinv**`(`

¶`y`,`a`) - :
`x`=**gammaincinv**`(`

¶`y`,`a`,`tail`) Compute the inverse of the normalized incomplete gamma function.

The normalized incomplete gamma function is defined as

x 1 / gammainc (x, a) = --------- | exp (-t) t^(a-1) dt gamma (a) / t=0

and

`gammaincinv (gammainc (`

for each non-negative value of`x`,`a`),`a`) =`x``x`. If`a`is scalar then`gammaincinv (`

is returned for each element of`y`,`a`)`y`and vice versa.If neither

`y`nor`a`is scalar then the sizes of`y`and`a`must agree, and`gammaincinv`

is applied element-by-element. The variable`y`must be in the interval*[0,1]*while`a`must be real and positive.By default,

`tail`is`"lower"`

and the inverse of the incomplete gamma function integrated from 0 to`x`is computed. If`tail`is`"upper"`

, then the complementary function integrated from`x`to infinity is inverted.The function is computed with Newton’s method by solving

`y`- gammainc (`x`,`a`) = 0Reference: A. Gil, J. Segura, and N. M. Temme, Efficient and accurate algorithms for the computation and inversion of the incomplete gamma function ratios, SIAM J. Sci. Computing, pp. A2965–A2981, Vol 34, 2012.

- :
`l`=**legendre**`(`

¶`n`,`x`) - :
`l`=**legendre**`(`

¶`n`,`x`,`normalization`) Compute the associated Legendre function of degree

`n`and order`m`= 0 …`n`.The value

`n`must be a real non-negative integer.`x`is a vector with real-valued elements in the range [-1, 1].The optional argument

`normalization`may be one of`"unnorm"`

,`"sch"`

, or`"norm"`

. The default if no normalization is given is`"unnorm"`

.When the optional argument

`normalization`is`"unnorm"`

, compute the associated Legendre function of degree`n`and order`m`and return all values for`m`= 0 …`n`. The return value has one dimension more than`x`.The associated Legendre function of degree

`n`and order`m`:m m 2 m/2 d^m P(x) = (-1) * (1-x ) * ---- P(x) n dx^m n

with Legendre polynomial of degree

`n`:1 d^n 2 n P(x) = ------ [----(x - 1) ] n 2^n n! dx^n

`legendre (3, [-1.0, -0.9, -0.8])`

returns the matrix:x | -1.0 | -0.9 | -0.8 ------------------------------------ m=0 | -1.00000 | -0.47250 | -0.08000 m=1 | 0.00000 | -1.99420 | -1.98000 m=2 | 0.00000 | -2.56500 | -4.32000 m=3 | 0.00000 | -1.24229 | -3.24000

When the optional argument

`normalization`is`"sch"`

, compute the Schmidt semi-normalized associated Legendre function. The Schmidt semi-normalized associated Legendre function is related to the unnormalized Legendre functions by the following:For Legendre functions of degree

`n`and order 0:0 0 SP(x) = P(x) n n

For Legendre functions of degree n and order m:

m m m 2(n-m)! 0.5 SP(x) = P(x) * (-1) * [-------] n n (n+m)!

When the optional argument

`normalization`is`"norm"`

, compute the fully normalized associated Legendre function. The fully normalized associated Legendre function is related to the unnormalized associated Legendre functions by the following:For Legendre functions of degree

`n`and order`m`m m m (n+0.5)(n-m)! 0.5 NP(x) = P(x) * (-1) * [-------------] n n (n+m)!

- :
`y`=**gammaln**`(`

¶`x`) - :
`y`=**lgamma**`(`

¶`x`) Return the natural logarithm of the gamma function of

`x`.Programming Note:

`lgamma`

is an alias for`gammaln`

and either name can be used in Octave.

- :
`y`=**psi**`(`

¶`z`) - :
`y`=**psi**`(`

¶`k`,`z`) Compute the psi (polygamma) function.

The polygamma functions are the

`k`th derivative of the logarithm of the gamma function. If unspecified,`k`defaults to zero. A value of zero computes the digamma function, a value of 1, the trigamma function, and so on.The digamma function is defined:

psi (z) = d (log (gamma (z))) / dx

When computing the digamma function (when

`k`equals zero),`z`can have any value real or complex value. However, for polygamma functions (`k`higher than 0),`z`must be real and non-negative.