# Functions

SpecialFunctions.erfFunction
erf(x)

Compute the error function of $x$, defined by

$\operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x \exp(-t^2) \; \mathrm{d}t \quad \text{for} \quad x \in \mathbb{C} \, .$
erf(x, y)

Accurate version of erf(y) - erf(x) (for real arguments only).

Implementation by

• Float32/Float64: C standard math library libm.
• BigFloat: C library for multiple-precision floating-point MPFR
source
SpecialFunctions.erfcFunction
erfc(x)

Compute the complementary error function of $x$, defined by

$\operatorname{erfc}(x) = 1 - \operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_x^\infty \exp(-t^2) \; \mathrm{d}t \quad \text{for} \quad x \in \mathbb{C} \, .$

This is the accurate version of 1-erf(x) for large $x$.

See also: erf(x).

Implementation by

• Float32/Float64: C standard math library libm.
• BigFloat: C library for multiple-precision floating-point MPFR
source
SpecialFunctions.erfcxFunction
erfcx(x)

Compute the scaled complementary error function of $x$, defined by

$\operatorname{erfcx}(x) = e^{x^2} \operatorname{erfc}(x) \quad \text{for} \quad x \in \mathbb{C} \, .$

This is the accurate version of $e^{x^2} \operatorname{erfc}(x)$ for large $x$. Note also that $\operatorname{erfcx}(-ix)$ computes the Faddeeva function w(x).

See also: erfc(x).

Implementation by

• Float32/Float64: C standard math library libm.
• BigFloat: MPFR has an open TODO item for this function until then, we use DLMF 7.12.1 for the tail.
source
SpecialFunctions.logerfcFunction
logerfc(x)

Compute the natural logarithm of the complementary error function of $x$, that is

$\operatorname{logerfc}(x) = \operatorname{ln}(\operatorname{erfc}(x)) \quad \text{for} \quad x \in \mathbb{R} \, .$

This is the accurate version of $\operatorname{ln}(\operatorname{erfc}(x))$ for large $x$.

See also: erfcx(x).

Implementation

Based on the erfc(x) and erfcx(x) functions. Currently only implemented for Float32, Float64, and BigFloat.

source
SpecialFunctions.logerfcxFunction
logerfcx(x)

Compute the natural logarithm of the scaled complementary error function of $x$, that is

$\operatorname{logerfcx}(x) = \operatorname{ln}(\operatorname{erfcx}(x)) \quad \text{for} \quad x \in \mathbb{R} \, .$

This is the accurate version of $\operatorname{ln}(\operatorname{erfcx}(x))$ for large and negative $x$.

See also: erfcx(x).

Implementation

Based on the erfc(x) and erfcx(x) functions. Currently only implemented for Float32, Float64, and BigFloat.

source
SpecialFunctions.erfiFunction
erfi(x)

Compute the imaginary error function of $x$, defined by

$\operatorname{erfi}(x) = -i \operatorname{erf}(ix) \quad \text{for} \quad x \in \mathbb{C} \, .$

See also: erf(x).

Implementation by

• Float32/Float64: C standard math library libm.
source
SpecialFunctions.dawsonFunction
dawson(x)

Compute the Dawson function (scaled imaginary error function) of $x$, defined by

$\operatorname{dawson}(x) = \frac{\sqrt{\pi}}{2} e^{-x^2} \operatorname{erfi}(x) \quad \text{for} \quad x \in \mathbb{C} \, .$

This is the accurate version of $\frac{\sqrt{\pi}}{2} e^{-x^2} \operatorname{erfi}(x)$ for large $x$.

See also: erfi(x).

Implementation by

• Float32/Float64: C standard math library libm.
source
SpecialFunctions.erfinvFunction
erfinv(x)

Compute the inverse error function of a real $x$, that is

$\operatorname{erfinv}(x) = \operatorname{erf}^{-1}(x) \quad \text{for} \quad x \in \mathbb{R} \, .$

See also: erf(x).

Implementation

Using the rational approximants tabulated in:

J. M. Blair, C. A. Edwards, and J. H. Johnson, "Rational Chebyshev approximations for the inverse of the error function", Math. Comp. 30, pp. 827–830 (1976). https://doi.org/10.1090/S0025-5718-1976-0421040-7, http://www.jstor.org/stable/2005402

combined with Newton iterations for BigFloat.

source
SpecialFunctions.erfcinvFunction
erfcinv(x)

Compute the inverse error complementary function of a real $x$, that is

$\operatorname{erfcinv}(x) = \operatorname{erfc}^{-1}(x) \quad \text{for} \quad x \in \mathbb{R} \, .$

See also: erfc(x).

Implementation

Using the rational approximants tabulated in:

J. M. Blair, C. A. Edwards, and J. H. Johnson, "Rational Chebyshev approximations for the inverse of the error function", Math. Comp. 30, pp. 827–830 (1976). https://doi.org/10.1090/S0025-5718-1976-0421040-7, http://www.jstor.org/stable/2005402

combined with Newton iterations for BigFloat.

source
SpecialFunctions.expintFunction
expint(z)
expint(ν, z)

Computes the exponential integral $\operatorname{E}_\nu(z) = \int_0^\infty \frac{e^{-zt}}{t^\nu} dt$. If $\nu$ is not specified, $\nu=1$ is used. Arbitrary complex $\nu$ and $z$ are supported.

source
SpecialFunctions.expintiFunction
expinti(x::Real)

Computes the exponential integral function $\operatorname{Ei}(x) = \int_{-\infty}^x \frac{e^t}{t} dt$, which is equivalent to $-\Re[\operatorname{E}_1(-x)]$ where $\operatorname{E}_1$ is the expint function.

source
SpecialFunctions.expintxFunction
expintx(z)
expintx(ν, z)

Computes the scaled exponential integral $\exp(z) \operatorname{E}_\nu(z) = e^z \int_0^\infty \frac{e^{-zt}}{t^\nu} dt$. If $\nu$ is not specified, $\nu=1$ is used. Arbitrary complex $\nu$ and $z$ are supported.

See also: expint(ν, z)

source
SpecialFunctions.sinintFunction
sinint(x)

Compute the sine integral function of $x$, defined by

$\operatorname{Si}(x) := \int_0^x \frac{\sin t}{t} \, \mathrm{d}t \quad \text{for} \quad x \in \mathbb{R} \,.$

See also: cosint(x).

Implementation

Using the rational approximants tabulated in:

A.J. MacLeod, "Rational approximations, software and test methods for sine and cosine integrals", Numer. Algor. 12, pp. 259–272 (1996). https://doi.org/10.1007/BF02142806, https://link.springer.com/article/10.1007/BF02142806.

Note: the second zero of $\text{Ci}(x)$ has a typo that is fixed: $r_1 = 3.38418 0422\mathbf{8} 51186 42639 78511 46402$ in the article, but is in fact: $r_1 = 3.38418 0422\mathbf{5} 51186 42639 78511 46402$.

source
SpecialFunctions.cosintFunction
cosint(x)

Compute the cosine integral function of $x$, defined by

$\operatorname{Ci}(x) := \gamma + \log x + \int_0^x \frac{\cos (t) - 1}{t} \, \mathrm{d}t \quad \text{for} \quad x > 0 \,,$

where $\gamma$ is the Euler-Mascheroni constant.

See also: sinint(x).

Implementation

Using the rational approximants tabulated in:

A.J. MacLeod, "Rational approximations, software and test methods for sine and cosine integrals", Numer. Algor. 12, pp. 259–272 (1996). https://doi.org/10.1007/BF02142806, https://link.springer.com/article/10.1007/BF02142806.

Note: the second zero of $\text{Ci}(x)$ has a typo that is fixed: $r_1 = 3.38418 0422\mathbf{8} 51186 42639 78511 46402$ in the article, but is in fact: $r_1 = 3.38418 0422\mathbf{5} 51186 42639 78511 46402$.

source
SpecialFunctions.airyaixFunction
airyaix(x)

Scaled Airy function of the first kind $\operatorname{Ai}(x) e^{\frac{2}{3} x \sqrt{x}}$. Throws DomainError for negative Real arguments.

source
SpecialFunctions.airyaiprimexFunction
airyaiprimex(x)

Scaled derivative of the Airy function of the first kind $\operatorname{Ai}'(x) e^{\frac{2}{3} x \sqrt{x}}$. Throws DomainError for negative Real arguments.

source
SpecialFunctions.sphericalbesseljFunction
sphericalbesselj(nu, x)

Spherical bessel function of the first kind at order nu, $j_ν(x)$. This is the non-singular solution to the radial part of the Helmholz equation in spherical coordinates.

source
SpecialFunctions.sphericalbesselyFunction
sphericalbessely(nu, x)

Spherical bessel function of the second kind at order nu, $y_ν(x)$. This is the singular solution to the radial part of the Helmholz equation in spherical coordinates. Sometimes known as a spherical Neumann function.

source
SpecialFunctions.besselhxFunction
besselhx(nu, [k=1,] z)

Compute the scaled Hankel function $\exp(∓iz) H_ν^{(k)}(z)$, where $k$ is 1 or 2, $H_ν^{(k)}(z)$ is besselh(nu, k, z), and $∓$ is $-$ for $k=1$ and $+$ for $k=2$. k defaults to 1 if it is omitted.

The reason for this function is that $H_ν^{(k)}(z)$ is asymptotically proportional to $\exp(∓iz)/\sqrt{z}$ for large $|z|$, and so the besselh function is susceptible to overflow or underflow when z has a large imaginary part. The besselhx function cancels this exponential factor (analytically), so it avoids these problems.

See also: besselh

source
SpecialFunctions.jincFunction
jinc(x)

Bessel function of the first kind divided by x. Following convention: $\operatorname{jinc}{x} = \frac{2 \cdot J_1{\pi x}}{\pi x}$. Sometimes known as sombrero or besinc function.

source
SpecialFunctions.ellipkFunction
ellipk(m)

Computes Complete Elliptic Integral of 1st kind $K(m)$ for parameter $m$ given by

$\operatorname{ellipk}(m) = K(m) = \int_0^{ \frac{\pi}{2} } \frac{1}{\sqrt{1 - m \sin^2 \theta}} \, \mathrm{d}\theta \quad \text{for} \quad m \in \left( -\infty, 1 \right] \, .$

See also: ellipe(m).

Arguments

• m: parameter $m$, restricted to the domain $(-\infty,1]$, is related to the elliptic modulus $k$ by $k^2=m$ and to the modular angle $\alpha$ by $k=\sin \alpha$.

Implementation

Using piecewise approximation polynomial as given in

'Fast Computation of Complete Elliptic Integrals and Jacobian Elliptic Functions', Fukushima, Toshio. (2014). F09-FastEI. Celest Mech Dyn Astr, DOI 10.1007/s10569-009-9228-z, https://pdfs.semanticscholar.org/8112/c1f56e833476b61fc54d41e194c962fbe647.pdf

For $m<0$, followed by

Fukushima, Toshio. (2014). 'Precise, compact, and fast computation of complete elliptic integrals by piecewise minimax rational function approximation'. Journal of Computational and Applied Mathematics. 282. DOI 10.13140/2.1.1946.6245., https://www.researchgate.net/publication/267330394

As suggested in this paper, the domain is restricted to $(-\infty,1]$.

source
SpecialFunctions.ellipeFunction
ellipe(m)

Computes Complete Elliptic Integral of 2nd kind $E(m)$ for parameter $m$ given by

$\operatorname{ellipe}(m) = E(m) = \int_0^{ \frac{\pi}{2} } \sqrt{1 - m \sin^2 \theta} \, \mathrm{d}\theta \quad \text{for} \quad m \in \left( -\infty, 1 \right] \, .$

See also: ellipk(m).

Arguments

• m: parameter $m$, restricted to the domain $(-\infty,1]$, is related to the elliptic modulus $k$ by $k^2=m$ and to the modular angle $\alpha$ by $k=\sin \alpha$.

Implementation

Using piecewise approximation polynomial as given in

'Fast Computation of Complete Elliptic Integrals and Jacobian Elliptic Functions', Fukushima, Toshio. (2014). F09-FastEI. Celest Mech Dyn Astr, DOI 10.1007/s10569-009-9228-z, https://pdfs.semanticscholar.org/8112/c1f56e833476b61fc54d41e194c962fbe647.pdf

For $m<0$, followed by

Fukushima, Toshio. (2014). 'Precise, compact, and fast computation of complete elliptic integrals by piecewise minimax rational function approximation'. Journal of Computational and Applied Mathematics. 282. DOI 10.13140/2.1.1946.6245., https://www.researchgate.net/publication/267330394

As suggested in this paper, the domain is restricted to $(-\infty,1]$.

source
SpecialFunctions.zetaFunction
zeta(s, z)

Generalized zeta function defined by

$\zeta(s, z)=\sum_{k=0}^\infty \frac{1}{((k+z)^2)^{s/2}},$

where any term with $k+z=0$ is excluded. For $\Re z > 0$, this definition is equivalent to the Hurwitz zeta function $\sum_{k=0}^\infty (k+z)^{-s}$.

The Riemann zeta function is recovered as $\zeta(s)=\zeta(s,1)$.

External links: Riemann zeta function, Hurwitz zeta function

source
zeta(s)

Riemann zeta function

$\zeta(s)=\sum_{n=1}^\infty \frac{1}{n^s}\quad\text{for}\quad s\in\mathbb{C}.$

source
SpecialFunctions.gammaMethod
gamma(z)

Compute the gamma function for complex $z$, defined by

$\Gamma(z) := \begin{cases} n! & \text{for} \quad z = n+1 \;, n = 0,1,2,\dots \\ \int_0^\infty t^{z-1} {\mathrm e}^{-t} \, {\mathrm d}t & \text{for} \quad \Re(z) > 0 \end{cases}$

and by analytic continuation in the whole complex plane.

See also: loggamma(z) for $\log \Gamma(z)$ and gamma(a,z) for the upper incomplete gamma function $\Gamma(a,z)$.

Implementation by

• Float: C standard math library libm.
• Complex: by exp(loggamma(z)).
• BigFloat: C library for multiple-precision floating-point MPFR
source
SpecialFunctions.loggammaMethod
loggamma(x)

Computes the logarithm of gamma for given x. If x is a Real, then it throws a DomainError if gamma(x) is negative.

If x is complex, then exp(loggamma(x)) matches gamma(x) (up to floating-point error), but loggamma(x) may differ from log(gamma(x)) by an integer multiple of $2\pi i$ (i.e. it may employ a different branch cut).

See also logabsgamma for real x.

source
SpecialFunctions.gammaMethod
gamma(a,x)

Returns the upper incomplete gamma function

$\Gamma(a,x) = \int_x^\infty t^{a-1} e^{-t} dt \,$

supporting arbitrary real or complex a and x.

(The ordinary gamma function gamma(x) corresponds to $\Gamma(a) = \Gamma(a,0)$. See also the gamma_inc function to compute both the upper and lower ($\gamma(a,x)$) incomplete gamma functions scaled by $\Gamma(a)$.

source
SpecialFunctions.loggammaMethod
loggamma(a,x)

Returns the log of the upper incomplete gamma function gamma(a,x):

$\log \Gamma(a,x) = \log \int_x^\infty t^{a-1} e^{-t} dt \,$

supporting arbitrary real or complex a and x.

If a and/or x is complex, then exp(loggamma(a,x)) matches gamma(a,x) (up to floating-point error), but loggamma(a,x) may differ from log(gamma(a,x)) by an integer multiple of $2\pi i$ (i.e. it may employ a different branch cut).

See also loggamma(x).

source
SpecialFunctions.gamma_incFunction
gamma_inc(a,x,IND=0)

Returns a tuple $(p, q)$ where $p + q = 1$, and $p=P(a,x)$ is the Incomplete gamma function ratio given by:

$P(a,x)=\frac{1}{\Gamma (a)} \int_{0}^{x} e^{-t}t^{a-1} dt.$

and $q=Q(a,x)$ is the Incomplete gamma function ratio given by:

$Q(x,a)=\frac{1}{\Gamma (a)} \int_{x}^{\infty} e^{-t}t^{a-1} dt.$

In terms of these, the lower incomplete gamma function is $\gamma(a,x) = P(a,x) \Gamma(a)$ and the upper incomplete gamma function is $\Gamma(a,x) = Q(a,x) \Gamma(a)$.

IND ∈ [0,1,2] sets accuracy: IND=0 means 14 significant digits accuracy, IND=1 means 6 significant digit, and IND=2 means only 3 digit accuracy.

SpecialFunctions.logabsbinomialFunction
logabsbinomial(n, k)
Accurate natural logarithm of the absolute value of the binomial coefficient binomial(n, k) for large n and k near n/2.
Returns a tuple (log(abs(binomial(n,k))), sign(binomial(n,k))).