Functions

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} e^{-t} \, \mathrm{d}t & \text{for} \quad \Re(z) > 0 \end{cases}\]

and by analytic continuation in the whole complex plane.

Examples

julia> gamma(0)
Inf

julia> gamma(1)
1.0

julia> gamma(2)
1.0

julia> gamma(0.5)^2 ≈ π
true

julia> gamma(4 + 1) == prod(1:4) == factorial(4)
true

External links: DLMF 5.2.1, 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

gamma(a,x)

Returns the upper incomplete gamma function

\[\Gamma(a,x) = \int_x^\infty t^{a-1} e^{-t} \mathrm{d}t\]

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 8.2.2, Wikipedia

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πi$ (i.e. it may employ a different branch cut).

See also logabsgamma for real x.

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} \mathrm{d}t\]

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

logabsgamma(x)

Compute the logarithm of absolute value of gamma for Real x and returns a tuple (log(abs(gamma(x))), sign(gamma(x))).

See also loggamma.

loggamma1p(x)

Compute $\log(\Gamma(1+x))$ for -1 < x <= 1.

logfactorial(x)

Compute the logarithmic factorial of a nonnegative integer x. Equivalent to loggamma of x + 1, but loggamma extends this function to non-integer x.

digamma(x)

Compute the digamma function of x (the logarithmic derivative of gamma(x)).

invdigamma(x)

Compute the inverse digamma function of x.

trigamma(x)

Compute the trigamma function of x (the logarithmic second derivative of gamma(x)).

polygamma(m, x)

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

External links: Wikipedia

See also: gamma(z)

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} \mathrm{d}t.\]

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

\[Q(a,x) = \frac{1}{\Gamma (a)} \int_{x}^{\infty} e^{-t}t^{a-1} \mathrm{d}t.\]

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 8.2.4, Wikipedia

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

gamma_inc_inv(a,p,q)

Inverts the gamma_inc(a,x) function, by computing x given a,p,q in $P(a,x) = p$ and $Q(a,x) = q$.

External links: DLMF 8.2.4, Wikipedia

See also: gamma_inc(a,x,ind).

loggammadiv(a,b)

Computes $\log(\Gamma(b)/\Gamma(a+b))$ when b >= 8

gammax(x)

\[\operatorname{gammax}(x) = \begin{cases} e^{\operatorname{stirling}(x)}\quad\quad\quad \text{if} \quad x>0,\\ \dfrac{\Gamma(x)}{\sqrt{2 \pi}e^{-x + (x-0.5)\log(x)}},\quad \text{if} \quad x\leq 0. \end{cases}\]

rgammax(a,x)

Evaluation of $1/\Gamma(a) e^{-x} x^{a}$. Based on DRCOMP from the NSWC library.

rgamma1pm1(a)

Computation of $1/\Gamma(a+1) - 1$ for -0.5<=a<=1.5: $1/\Gamma (a+1) - 1$.

Uses the relation gamma(a+1) = a*gamma(a).

gamma_inc_asym(a, x, ind)

Compute $P(a,x)$ using asymptotic expansion given by:

\[R(a,x)/a * (1 + \sum_{n=1}^{N-1}(a_{n}/x^{n} + \Theta _{n}a_{n}/x^{n}))\]

where R(a,x) = rgammax(a,x). Used when 1 <= a <= BIG and x >= x0.

External links: DLMF 8.11.2

See also: gamma_inc(a,x,ind)

gamma_inc_cf(a, x, ind)

Computes $P(a,x)$ by continued fraction expansion given by:

\[R(a,x) \cdot \cfrac{1}{1-\cfrac{z}{a+1+\cfrac{z}{a+2-\cfrac{(a+1)z}{a+3+\cfrac{2z}{a+4-\cfrac{(a+2)z}{a+5+\cfrac{3z}{a+6-\ddots}}}}}}}\]

Used when 1 <= a <= BIG and x < x0.

External links: DLMF 8.9.2

See also: gamma_inc(a,x,ind)

gamma_inc_fsum(a,x)

Compute using Finite Sums for $Q(a,x)$ when a >= 1 && 2a is integer. Used when a <= x <= x0 and a = n/2. Refer Eqn (14) in the paper.

See also: gamma_inc(a,x,ind)

gamma_inc_minimax(a,x,z)

Compute $P(a,x)$ using minimax approximations given by :

\[1/2 * \operatorname{erfc}(\sqrt{y}) - e^{-y}/\sqrt{2\pi a} ⋅ T(a,\lambda)\]

where

\[T(a,\lambda) = \sum_{0}^{N} c_{k}(z)a^{-k}\]

This is a higher accuracy approximation of Temme expansion, which deals with the region near a ≈ x with a large. Refer Appendix F in the paper for the extensive set of coefficients calculated using Brent's multiple precision arithmetic(set at 50 digits) in

Brent, R. P. A Fortran multiple-precision arithmetic package, ACM Trans. Math. Softw. 4(1978), 57-70, doi: 10.1145/355769.355775.

External links: DLMF 8.12.8

See also: gamma_inc(a,x,ind)

gamma_inc_taylor(a, x, ind)

Compute $P(a,x)$ using Taylor Series for P/R given by:

\[R(a,x)/a * (1 + \sum_{n=1}^{\infty} x^{n}/((a+1)(a+2)...(a+n)))\]

Used when 1 <= a <= BIG and x <= max{a, ln 10}.

External links: DLMF 8.11.2

See also: gamma_inc(a,x,ind)

gamma_inc_taylor_x(a,x,ind)

Computes $P(a,x)$ based on Taylor expansion of $P(a,x)/x^a$ given by:

\[J = -a * \sum_{1}^{\infty} (-x)^{n}/((a+n)n!)\]

and $P(a,x)/x^a$ is given by:

\[(1 - J) / \Gamma(a+1)\]

resulting from term-by-term integration of gamma_inc(a,x,ind). This is used when a < 1 and x < 1.1 - Refer Eqn (9) in the paper.

See also: gamma_inc(a,x,ind)

gamma_inc_temme(a, x, z)

Compute $P(a,x)$ using Temme's expansion given by :

\[1/2 * \operatorname{erfc}(\sqrt{y}) - e^{-y}/\sqrt{2\pi a} ⋅ T(a,\lambda)\]

where

\[T(a,\lambda) = \sum_{0}^{N} c_{k}(z)a^{-k}\]

This mainly solves the problem near the region when a ≈ x with a large, and is a lower accuracy version of the minimax approximation.

External links: DLMF 8.12.8

See also: gamma_inc(a,x,ind)

gamma_inc_temme_1(a, x, z, ind)

Computes $P(a,x)$ using simplified Temme expansion near $y=0$ by:

\[E(y) - (1 - y)/\sqrt{2\pi a} ⋅ T(a,\lambda)\]

where

\[E(y) = 1/2 - (1 - y/3) ⋅ \sqrt{y/\pi}\]

Used instead of it's previous function when $\sigma \leq e_{0}/\sqrt{a}$.

External links: DLMF 8.12.8

gamma_inc_inv_psmall(a, logr)

Compute x0 - initial approximation when p is small. Here we invert the series in Eqn (2.20) in the paper and write the inversion problem as:

\[x = r\left[1 + a\sum_{k=1}^{\infty}\frac{(-x)^{n}}{(a+n)n!}\right]^{-1/a},\]

where $r = (p\Gamma(1+a))^{1/a}$ Inverting this relation we obtain $x = r + \sum_{k=2}^{\infty} c_{k} r^{k}$.

gamma_inc_inv_qsmall(a, q, qgammaxa)

Compute x0 - initial approximation when q is small from $e^{-x_{0}} x_{0}^{a} = q \Gamma(a)$. Asymptotic expansions Eqn (2.29) in the paper is used here and higher approximations are obtained using

\[x \sim x_{0} - L + b \sum_{k=1}^{\infty} d_{k}/x_{0}^{k}\]

where $b = 1-a$, $L = \ln{x_0}$.

gamma_inc_inv_alarge(a, minpq, pcase)

Compute x0 - initial approximation when a is large. The inversion problem is rewritten as:

\[0.5 \operatorname{erfc}(\eta \sqrt{a/2}) + R_{a}(\eta) = q\]

For large values of a we can write: $\eta(q,a) = \eta_{0}(q,a) + \epsilon(\eta_{0},a)$ and it is possible to expand:

\[\epsilon(\eta_{0},a) = \epsilon_{1}(\eta_{0},a)/a + \epsilon_{2}(\eta_{0},a)/a^{2} + \epsilon_{3}(\eta_{0},a)/a^{3} + \cdots\]

which is calculated by coeff1, coeff2 and coeff3 functions below. This returns a tuple (x0,fp), where fp is computed since it's an approximation for the coefficient after inverting the original power series.

auxgam(x)

Compute function g in $1/\Gamma(x+1) = 1 + x (x-1) g(x)$, -1 <= x <= 1.

beta(x, y)

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

logbeta(x, y)

Natural logarithm of the beta function $\log(|\operatorname{B}(x,y)|)$.

See also logabsbeta.

logabsbeta(x, y)

Compute the natural logarithm of the absolute value of the beta function, returning a tuple (log(abs(beta(x,y))), sign(beta(x,y)))

See also logbeta.

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

External links Base.binomial

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} \mathrm{d}t,\]

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

External links: DLMF 8.17.1, Wikipedia

See also: beta_inc_inv

beta_inc_inv(a, b, p, q=1-p)

Return a tuple (x, 1-x) where x satisfies $I_{x}(a, b) = p$, i.e., x is the inverse of the regularized incomplete beta function $I_{x}(a, b)$.

See also: beta_inc

ncbeta(a,b,lambda,x)

Compute the CDF of the noncentral beta distribution given by

\[I_{x}(a,b; \lambda) = \sum_{j=0}^{\infty} q(\lambda/2,j) I_{x}(a+j,b;0)\]

For $\lambda < 54$ : algorithm suggested by Lenth(1987) in ncbeta_tail(a,b,lambda,x). Else for $\lambda \geq 54$: modification in Chattamvelli(1997) in ncbeta_poisson(a,b,lambda,x) by using both forward and backward recurrences.

ncbeta_poisson(a,b,lambda,x)

Compute CDF of noncentral beta if lambda >= 54 using: First $\lambda/2$ is calculated and the Poisson term is calculated using $P(j-1) = j/\lambda P(j)$ and $P(j+1) = \lambda/(j+1) P(j)$. Then backward recurrences are used until either the Poisson weights fall below errmax or iterlo is reached.

\[I_{x}(a+j-1,b) = I_{x}(a+j,b) + \Gamma(a+b+j-1)/\Gamma(a+j)\Gamma(b)x^{a+j-1}(1-x)^{b}\]

Then forward recurrences are used until error bound falls below errmax.

\[I_{x}(a+j+1,b) = I_{x}(a+j,b) - \Gamma(a+b+j)/\Gamma(a+j)\Gamma(b)x^{a+j}(1-x)^{b}\]

ncbeta_tail(x,a,b,lambda)

Compute tail of the noncentral beta distribution. Uses the recursive relation

\[I_{x}(a,b+1;0) = I_{x}(a,b;0) - \Gamma(a+b)/\Gamma(a+1)\Gamma(b) x^a (1-x)^b\]

and $\Gamma(a+1) = a\Gamma(a)$ given in DLMF 8.17.21.

beta_inc_power_series1(a,b,x,epps)

Another variant of BPSER(A,B,X,EPS).

External links: DLMF 8.17.22, Wikipedia

See also: beta_inc

Implementation

APSER(A,B,X,EPS) from Didonato and Morris (1982)

beta_inc_power_series2(a,b,x,epps)

Variant of BPSER(A,B,X,EPS).

External links: DLMF 8.17.22, Wikipedia

See also: beta_inc

Implementation

FPSER(A,B,X,EPS) from Didonato and Morris (1982)

beta_inc_asymptotic_asymmetric(a, b, x, y, w, epps)

Evaluation of $I_{x}(a,b)$ using asymptotic expansion. It is assumed a >= 15 and b <= 1, and epps is tolerance used.

External links: DLMF 8.17.22, Wikipedia

See also: beta_inc

beta_inc_asymptotic_symmetric(a,b,lambda,epps)

Compute $I_{x}(a,b)$ using asymptotic expansion for a, b >= 15. It is assumed that $\lambda = (a+b)*y - b$.

External links: DLMF 8.17.22, Wikipedia

See also: beta_inc

Implementation

BASYM(A,B,LAMBDA,EPS) from Didonato and Morris (1982)

beta_inc_power_series(a, b, x, epps)

Computes $I_x(a,b)$ using power series:

\[I_{x}(a,b) = G(a,b) x^{a}/a \left[1 + a \sum_{j=1}^{\infty} ((1-b)(2-b)\dots(j-b)/j!(a+j)) x^{j}\right]\]

External links: DLMF 8.17.22, Wikipedia

See also: beta_inc

Implementation

BPSER(A,B,X,EPS) from Didonato and Morris (1982)

beta_inc_diff(a, b, x, y, n, epps)

Compute $I_{x}(a,b) - I_{x}(a+n,b)$ where n is positive integer and epps is tolerance. A generalised version of DLMF 8.17.20.

External links: DLMF 8.17.20, Wikipedia

See also: beta_inc

beta_inc_cont_fraction(a,b,x,y,lambda,epps)

Compute $I_{x}(a,b)$ using continued fraction expansion when a, b > 1. It is assumed that $\lambda = (a+b)*y - b$

External links: DLMF 8.17.22, Wikipedia

See also: beta_inc

Implementation

BFRAC(A,B,X,Y,LAMBDA,EPS) from Didonato and Morris (1982)

beta_integrand(a, b, x, y, mu=0.0)

Compute $e^{\mu} x^a y^b / B(a,b)$

chepolsum(n,x,a)

Computes a series of Chebyshev Polynomials given by: a[1]/2 + a[2]T1(x) + .... + a[n]T{n-1}(X).

lambdaeta(eta)

Compute the value of $\lambda$ satisfying $\eta^{2}/2 = \lambda-1-\log{\lambda}$.

stirling_corr(a0,b0)

Compute stirling(a0) + stirling(b0) - stirling(a0 + b0) for a0, b0 >= 8

stirling_error(x)

Compute $\ln{\Gamma(x)} - (x-0.5)*\ln{x} + x - \ln{(2\pi)}/2$. Adapted from stirling in IncgamFI.

esum(mu,x)

Compute $e^{\mu+x}$

Computing the first coefficient for the expansion:

\[\epsilon (\eta_{0},a) = \epsilon_{1}(\eta_{0},a)/a + \epsilon_{2}(\eta_{0},a)/a^{2} + \epsilon_{3}(\eta_{0},a)/a^{3}\]

Refer Eqn (3.12) in the paper

Computing the second coefficient for the expansion:

\[\epsilon (\eta_{0},a) = \epsilon_{1} (\eta_{0},a)/a + \epsilon_{2} (\eta_{0},a)/a^{2} + \epsilon_{3} (\eta_{0},a)/a^{3}\]

Refer Eqn (3.12) in the paper

Computing the third and last coefficient for the expansion:

\[\epsilon(\eta_{0},a) = \epsilon_{1}(\eta_{0},a)/a + \epsilon_{2}(\eta_{0},a)/a^{2} + \epsilon_{3}(\eta_{0},a)/a^{3}\]

Refer Eqn (3.12) in the paper

ncF(x,v1,v2,lambda)

Compute CDF of noncentral F distribution given by:

\[F(x, v_1, v_2; \lambda) = I_{v_1 x/(v_1 x + v_2)}(v_1/2, v_2/2; \lambda)\]

where $I_{x}(a,b; \lambda)$ is the noncentral beta function computed above.

External links: Wikipedia

expint(z)
expint(ν, z)

Computes the exponential integral

\[\operatorname{E}_\nu(z) = \int_1^\infty \frac{e^{-zt}}{t^\nu} \mathrm{d}t.\]

If $\nu$ is not specified, $\nu=1$ is used. Arbitrary complex $\nu$ and $z$ are supported.

External links: DLMF 8.19, Wikipedia

expinti(x::Real)

Computes the exponential integral function

\[\operatorname{Ei}(x) = \int_{-\infty}^x \frac{e^t}{t} \mathrm{d}t,\]

which is equivalent to $-\Re[\operatorname{E}_1(-x)]$ where $\operatorname{E}_1$ is the expint 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} \mathrm{d}t.\]

If $\nu$ is not specified, $\nu = 1$ is used. Arbitrary complex $\nu$ and $z$ are supported.

See also: expint(ν, z)

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

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

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

External links: DLMF 7.2.1, 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

erf(x, y)

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

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 7.2.2, Wikipedia.

See also: erf(x).

Implementation by

• Float32/Float64: C standard math library libm.
• BigFloat: C library for multiple-precision floating-point MPFR

logerf(x, y)

Compute the natural logarithm of two-argument error function. This is an accurate version of log(erf(x, y)), which works for large x, y.

External links: Wikipedia.

See also: erf(x,y).

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.

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 7.2.3, 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.

logerfc(x)

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

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

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

External links: Wikipedia.

See also: erfcx(x).

Implementation

Based on the erfc(x) and erfcx(x) functions.

logerfcx(x)

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

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

This is the accurate version of $\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.

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.

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.

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 7.2.5, Wikipedia.

See also: erfi(x).

Implementation by

• Float32/Float64: C standard math library libm.

faddeeva(z)

Compute the Faddeeva function of complex z, defined by $e^{-z^2} \operatorname{erfc}(-iz)$. Note that this function, also named w (original Faddeeva package) or wofz (Scilab package), is equivalent to $\operatorname{erfcx}(-iz)$.

airyai(x)

Airy function of the first kind $\operatorname{Ai}(x)$.

External links: DLMF 9.2, Wikipedia

See also: airyaix, airyaiprime, airybi

airyaiprime(x)

Derivative of the Airy function of the first kind $\operatorname{Ai}'(x)$.

External links: DLMF 9.2, Wikipedia

See also: airyaiprimex, airyai, airybi

airybi(x)

Airy function of the second kind $\operatorname{Bi}(x)$.

External links: DLMF 8.2, Wikipedia

See also: airybix, airybiprime, airyai

airybiprime(x)

Derivative of the Airy function of the second kind $\operatorname{Bi}'(x)$.

External links: DLMF 9.2, Wikipedia

See also: airybiprimex, airybi, airyai

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.

External links: DLMF 9.2, Wikipedia

See also: airyai, airyaiprime, airybi

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.

External links: DLMF 9.2, Wikipedia

See also: airyaiprime, airyai, airybi

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

External links: DLMF 9.2, Wikipedia

See also: airybi, airybiprime, airyai

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

External links: DLMF 9.2, Wikipedia

See also: airybiprime, airybi, airyai

besselj(nu, x)

Bessel function of the first kind of order nu, $J_\nu(x)$.

External links: DLMF 10.2.2, Wikipedia

See also: besseljx(nu,x), besseli(nu,x), besselk(nu,x)

besselj0(x)

Bessel function of the first kind of order 0, $J_0(x)$.

External links: DLMF 10.2.2, Wikipedia

See also: besselj(nu,x)

besselj1(x)

Bessel function of the first kind of order 1, $J_1(x)$.

External links: DLMF 10.2.2, Wikipedia

See also: besselj(nu,x)

besseljx(nu, x)

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

External links: DLMF 10.2.2, Wikipedia

See also: besselj(nu,x), besseli(nu,x), besselk(nu,x)

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.

bessely(nu, x)

Bessel function of the second kind of order nu, $Y_\nu(x)$.

External links: DLMF 10.2.3, Wikipedia

See also besselyx(nu,x) for a scaled variant.

bessely0(x)

Bessel function of the second kind of order 0, $Y_0(x)$.

External links: DLMF 10.2.3, Wikipedia

See also: bessely(nu,x)

bessely1(x)

Bessel function of the second kind of order 1, $Y_1(x)$.

External links: DLMF 10.2.3, Wikipedia

See also: bessely(nu,x)

besselyx(nu, x)

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

External links: DLMF 10.2.3, Wikipedia

See also bessely(nu,x)

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.

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.

External links: DLMF 10.2.5 and DLMF 10.2.6, Wikipedia

See also: besselhx for an exponentially scaled variant.

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 10.2, Wikipedia

See also: besselh

hankelh1(nu, x)

Bessel function of the third kind of order nu, $H^{(1)}_\nu(x)$.

External links: DLMF 10.2.5, Wikipedia

See also: hankelh1x

hankelh1x(nu, x)

Scaled Bessel function of the third kind of order nu, $H^{(1)}_\nu(x) e^{-x i}$.

External links: DLMF 10.2.5, Wikipedia

See also: hankelh1

hankelh2(nu, x)

Bessel function of the third kind of order nu, $H^{(2)}_\nu(x)$.

External links: DLMF 10.2.6, Wikipedia

See also: hankelh2x(nu,x)

hankelh2x(nu, x)

Scaled Bessel function of the third kind of order nu, $H^{(2)}_\nu(x) e^{x i}$.

External links: DLMF 10.2.6, Wikipedia

See also: hankelh2(nu,x)

besseli(nu, x)

Modified Bessel function of the first kind of order nu, $I_\nu(x)$.

External links: DLMF 10.25.2, Wikipedia

See also: besselix(nu,x), besselj(nu,x), besselk(nu,x)

besselix(nu, x)

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

External links: DLMF 10.25.2, Wikipedia

See also: besseli(nu,x), besselj(nu,x), besselk(nu,x)

besselk(nu, x)

Modified Bessel function of the second kind of order nu, $K_\nu(x)$.

External links: DLMF 10.25.3, Wikipedia

See also: besselkx(nu,x), besseli(nu,x), besselj(nu,x)

besselkx(nu, x)

Scaled modified Bessel function of the second kind of order nu, $K_\nu(x) e^x$.

External links: DLMF 10.25.3, Wikipedia

See also: besselk(nu,x), besseli(nu,x), besselj(nu,x)

jinc(x)

Bessel function of the first kind divided by x. Following convention:

\[\operatorname{jinc}{x} = \frac{2 J_1({\pi x})}{\pi x}.\]

Sometimes known as sombrero or besinc function.

External links: Wikipedia

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 19.2.8, 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]$.

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 19.2.8, 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]$.

eta(s)

Dirichlet eta function

\[\eta(s) = \sum_{n=1}^\infty (-1)^{n-1} / n^{s}.\]

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

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