Edit on GitHub

Functions

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.beta — Method
beta(x, y)

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

source
SpecialFunctions.beta_inc — Method
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_inc_cont_fraction — Method
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, Wikipedia

See also: beta_inc

Implementation

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

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

source
SpecialFunctions.beta_inc_power_series — Method
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 (1 + a\sum_{j=1}^{\infty}((1-b)(2-b)...(j-b)/j!(a+j)) x^{j})\]

External links: DLMF, Wikipedia

See also: beta_inc

Implementation

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

source
SpecialFunctions.coeff1 — Method

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

source
SpecialFunctions.coeff2 — Method

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

source
SpecialFunctions.coeff3 — Method

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

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.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.ellipe — Method
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.ellipk — Method
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.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.erfcinv — Method
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.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.erfinv — Method
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.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 — Method
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.faddeeva — Function
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)$.

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

Examples

Special Values:

julia> gamma(0)
Inf

julia> gamma(1)
1.0

julia> gamma(2)
1.0

$\Gamma(0.5)^2 = \pi$

julia> gamma(0.5)^2
3.1415926535897936

julia> gamma(0.5)^2 ≈ π
true

For integer n: $\Gamma(n+1) = prod(1:n) = factorial(n)$

julia> gamma(4+1)
24.0

julia> prod(1:4)  # == 1*2*3*4
24

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

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.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(a,x)=\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.gamma_inc_asym — Method
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

See also: gamma_inc(a,x,ind)

source
SpecialFunctions.gamma_inc_cf — Method
gamma_inc_cf(a, x, ind)

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

\[R(a,x) * \frac{1}{1-\frac{z}{a+1+\frac{z}{a+2-\frac{(a+1)z}{a+3+\frac{2z}{a+4-\frac{(a+2)z}{a+5+\frac{3z}{a+6-\dots}}}}}}}\]

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

External links: DLMF

See also: gamma_inc(a,x,ind)

source
SpecialFunctions.gamma_inc_inv_alarge — Method
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} + ...\]

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.

source
SpecialFunctions.gamma_inc_inv_psmall — Method
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}$.

source
SpecialFunctions.gamma_inc_inv_qsmall — Method
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}$.

source
SpecialFunctions.gamma_inc_minimax — Method
gamma_inc_minimax(a,x,z)

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

\[1/2 * 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 .

External links: DLMF

See also: gamma_inc(a,x,ind)

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

source
SpecialFunctions.gamma_inc_temme — Method
gamma_inc_temme(a, x, z)

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

\[1/2 * 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

See also: gamma_inc(a,x,ind)

source
SpecialFunctions.gamma_inc_temme_1 — Method
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 <= e_{0}/\sqrt{a}$.

External links: DLMF

source
SpecialFunctions.gammax — Method
gammax(x)

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

source
SpecialFunctions.jinc — Method
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.logabsbinomial — Method
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
SpecialFunctions.logerf — Method
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).

source
SpecialFunctions.logerfc — Method
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 — Method
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.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.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.ncF — Method
ncF(x,v1,v2,lambda)

Compute CDF of noncentral F distribution given by:

\[F(x, v1, v2; lambda) = I_{v1*x/(v1*x + v2)}(v1/2, v2/2; \lambda)\]

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

Wikipedia: https://en.wikipedia.org/wiki/Noncentral_F-distribution

source
SpecialFunctions.ncbeta — Method
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 >= 54$ : modification in Chattamvelli(1997) in ncbeta_poisson(a,b,lambda,x) by using both forward and backward recurrences.

source
SpecialFunctions.ncbeta_poisson — Method
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}\]

source
SpecialFunctions.ncbeta_tail — Method
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 https://dlmf.nist.gov/8.17.21.

source
SpecialFunctions.rgamma1pm1 — Method

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

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.sphericalbesselj — Method
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 — Method
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.zeta — Method
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