Edit on GitHub

Functions

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

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

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). http://dx.doi.org/10.1007/BF02142806, https://link.springer.com/article/10.1007/BF02142806.

Note: the second zero of $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). http://dx.doi.org/10.1007/BF02142806, https://link.springer.com/article/10.1007/BF02142806.

Note: the second zero of $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.polygamma — Function
polygamma(m, x)

Compute the polygamma function of order m of argument x (the (m+1)th derivative of the logarithm of gamma(x))

source
SpecialFunctions.airyaix — Function
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.airyaiprimex — Function
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.airybix — Function
airybix(x)

Scaled Airy function of the second kind $\operatorname{Bi}(x) e^{- \left| \operatorname{Re} \left( \frac{2}{3} x \sqrt{x} \right) \right|}$.

source
SpecialFunctions.airybiprimex — Function
airybiprimex(x)

Scaled derivative of the Airy function of the second kind $\operatorname{Bi}'(x) e^{- \left| \operatorname{Re} \left( \frac{2}{3} x \sqrt{x} \right) \right|}$.

source
SpecialFunctions.besseljx — Function
besseljx(nu, x)

Scaled Bessel function of the first kind of order nu, $J_\nu(x) e^{- | \operatorname{Im}(x) |}$.

source
SpecialFunctions.besselyx — Function
besselyx(nu, x)

Scaled Bessel function of the second kind of order nu, $Y_\nu(x) e^{- | \operatorname{Im}(x) |}$.

source
SpecialFunctions.besselh — Function
besselh(nu, [k=1,] x)

Bessel function of the third kind of order nu (the Hankel function). k is either 1 or 2, selecting hankelh1 or hankelh2, respectively. k defaults to 1 if it is omitted. (See also besselhx for an exponentially scaled variant.)

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.

source
SpecialFunctions.besselix — Function
besselix(nu, x)

Scaled modified Bessel function of the first kind of order nu, $I_\nu(x) e^{- | \operatorname{Re}(x) |}$.

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

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

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.zeta — Function
zeta(s, z)

Generalized zeta function $\zeta(s, z)$, defined by the sum $\sum_{k=0}^\infty ((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}$. For $z=1$, it yields the Riemann zeta function $\zeta(s)$.

source
zeta(s)

Riemann zeta function $\zeta(s)$.

source
SpecialFunctions.gamma — Function
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_inc — Function
gamma_inc(a,x,IND)

DLMF: https://dlmf.nist.gov/8.2#E4 , https://dlmf.nist.gov/8.2#E5 Wikipedia: https://en.wikipedia.org/wiki/Incompletegammafunction IND –> Accuracy desired ; IND=0 means 14 significant digits accuracy , IND=1 means 6 significant digit and IND=2 means only 3 digit accuracy suffices. gammainc(a,x) or P(a,x) is the Incomplete gamma function ratio given by : ``1/\Gamma (a) \int{0}^{x} e^{-t}t^{a-1} dt$gamma_q(a,x) or Q(a,x) is the Incomplete gamma function ratio given by : 1 - P(a,x) ->$1/\Gamma (a) \int{x}^{\infty} e^{-t}t^{a-1} dt`` Returns a tuple (gammainc, gammaq) where gammainc + gamma_q = 1.0

source
SpecialFunctions.gamma_inc_inv — Function
gamma_inc_inv(a,p,q)

DLMF: https://dlmf.nist.gov/8.2#E4 , https://dlmf.nist.gov/8.2#E5 Wiki: https://en.wikipedia.org/wiki/Incompletegammafunction

gammainc(a,x) or (P(a,x),Q(a,x)) is the Incomplete gamma function ratio given by : ``1/\Gamma (a) \int{0}^{x} e^{-t}t^{a-1} dt`gamma_inc_inv(a,p,q) inverts the gamma_inc function, by computingxgivena,p,q` in P(a,x)=p and Q(a,x)=q.

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

Computes Incomplete Beta Function Ratios given by:

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

and $I_{y}(a,b) = 1.0 - I_{x}(a,b)$. given $B(a,b) = 1/G(a,b) = \Gamma(a)\Gamma(b)/\Gamma(a+b)$ and $x+y = 1$.

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