Functions

SpecialFunctions.erfcFunction
erfc(x)

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

\[\operatorname{erfc}(x) = 1 - \operatorname{erf}(x) = \frac{2}{\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.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).

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.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$.

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.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$.

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

External links: Wikipedia.

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$.

External links: DLMF, Wikipedia.

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

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). http://dx.doi.org/10.1090/S0025-5718-1976-0421040-7, http://www.jstor.org/stable/2005402

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

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). http://dx.doi.org/10.1090/S0025-5718-1976-0421040-7, http://www.jstor.org/stable/2005402

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

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). http://dx.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.

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). http://dx.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.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.

External links: DLMF, Wikipedia

See also: besselh

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

External links: DLMF, Wikipedia.

See also: ellipe(m).

Arguments

  • m: parameter $m$ is in relation to elliptic modulus $k$ by $k^2=m$ and 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

Also suggested in this paper that we should consider domain only from $(-\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] \, .\]

External links: DLMF, Wikipedia.

See also: ellipk(m).

Arguments

  • m: parameter $m$ is in relation to elliptic modulus $k$ by $k^2=m$ and 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

Also suggested in this paper that we should consider domain only from $(-\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}.\]

External links: Wikipedia

source
SpecialFunctions.gammaFunction
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).

Implementation by

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

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.\]

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_incFunction
beta_inc(a,b,x)

Returns a tuple $(I_{x}(a,b),1.0-I_{x}(a,b))$ where the Regularized Incomplete Beta Function is 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(a,b,p,q)

source
SpecialFunctions.betaFunction
beta(x, y)

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

source
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))).

source