Edit on GitHub

Functions

SpecialFunctions.erf — Function
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).

External links: DLMF, Wikipedia.

See also: erfc(x), erfcx(x), erfi(x), dawson(x), erfinv(x), erfcinv(x).

Implementation by

  • Float32/Float64: C standard math library libm.
  • BigFloat: C library for multiple-precision floating-point MPFR
source
SpecialFunctions.erfc — Function
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$.

External links: DLMF, Wikipedia.

See also: erf(x).

Implementation by

  • Float32/Float64: C standard math library libm.
  • BigFloat: C library for multiple-precision floating-point MPFR
source
SpecialFunctions.erfcx — Function
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).

External links: DLMF, Wikipedia.

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.logerfc — Function
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$.

External links: Wikipedia.

See also: erfcx(x).

Implementation

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

source
SpecialFunctions.logerfcx — Function
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$.

External links: Wikipedia.

See also: erfcx(x).

Implementation

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

source
SpecialFunctions.erfi — Function
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} \, .\]

External links: Wikipedia.

See also: erf(x).

Implementation by

  • Float32/Float64: C standard math library libm.
source
SpecialFunctions.dawson — Function
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$.

External links: DLMF, Wikipedia.

See also: erfi(x).

Implementation by

  • Float32/Float64: C standard math library libm.
source
SpecialFunctions.erfinv — Function
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} \, .\]

External links: Wikipedia.

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.erfcinv — Function
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} \, .\]

External links: Wikipedia.

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.expint — Function
expint(z)
expint(ν, z)

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

External links: DLMF, Wikipedia

source
SpecialFunctions.expinti — Function
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.expintx — Function
expintx(z)
expintx(ν, z)

Computes the scaled exponential integral $\exp(z) \operatorname{E}_\nu(z) = e^z \int_1^\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.sinint — Function
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} \,.\]

External links: DLMF, Wikipedia.

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.cosint — Function
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.

External links: DLMF, Wikipedia.

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.sphericalbesselj — Function
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.sphericalbessely — Function
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.besselhx — Function
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.

External links: DLMF, Wikipedia

See also: besselh

source
SpecialFunctions.jinc — Function
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.

External links: Wikipedia

source
SpecialFunctions.ellipk — Function
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] \, .\]

External links: DLMF, Wikipedia.

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.ellipe — Function
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] \, .\]

External links: DLMF, Wikipedia.

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.zeta — Function
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}.\]

External links: Wikipedia

source
SpecialFunctions.gamma — Method
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.

External links: DLMF, Wikipedia.

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.loggamma — Method
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.gamma — Method
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)$.

External links: DLMF, Wikipedia

source
SpecialFunctions.loggamma — Method
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_inc — Function
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.

External links: DLMF, Wikipedia

See also gamma(z), gamma_inc_inv(a,p,q)

source
SpecialFunctions.beta_inc — Function
beta_inc(a, b, x, y=1-x)

Return a tuple $(I_{x}(a,b), 1-I_{x}(a,b))$ where $I_{x}(a,b)$ is the regularized incomplete beta function given by

\[I_{x}(a,b) = \frac{1}{B(a,b)} \int_{0}^{x} t^{a-1}(1-t)^{b-1} dt,\]

where $B(a,b) = \Gamma(a)\Gamma(b)/\Gamma(a+b)$.

External links: DLMF, Wikipedia

See also: beta_inc_inv

source
SpecialFunctions.beta — Function
beta(x, y)

Euler integral of the first kind $\operatorname{B}(x,y) = \Gamma(x)\Gamma(y)/\Gamma(x+y)$.

source
SpecialFunctions.logabsbinomial — Function
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))).

source