Functions
SpecialFunctions.airyai — Functionairyai(x)Airy function of the first kind $\operatorname{Ai}(x)$.
External links: DLMF, Wikipedia
See also: airyaix, airyaiprime, airybi
SpecialFunctions.airyaiprime — Functionairyaiprime(x)Derivative of the Airy function of the first kind $\operatorname{Ai}'(x)$.
External links: DLMF, Wikipedia
See also: airyaiprimex, airyai, airybi
SpecialFunctions.airyaiprimex — Functionairyaiprimex(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, Wikipedia
See also: airyaiprime, airyai, airybi
SpecialFunctions.airyaix — Functionairyaix(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, Wikipedia
See also: airyai, airyaiprime, airybi
SpecialFunctions.airybi — Functionairybi(x)Airy function of the second kind $\operatorname{Bi}(x)$.
External links: DLMF, Wikipedia
See also: airybix, airybiprime, airyai
SpecialFunctions.airybiprime — Functionairybiprime(x)Derivative of the Airy function of the second kind $\operatorname{Bi}'(x)$.
External links: DLMF, Wikipedia
See also: airybiprimex, airybi, airyai
SpecialFunctions.airybiprimex — Functionairybiprimex(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, Wikipedia
See also: airybiprime, airybi, airyai
SpecialFunctions.airybix — Functionairybix(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, Wikipedia
See also: airybi, airybiprime, airyai
SpecialFunctions.auxgam — Methodauxgam(x)Compute function g in $1/\Gamma(x+1) = 1+x*(x-1)*g(x)$, -1 <= x <= 1.
SpecialFunctions.besselh — Functionbesselh(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 and DLMF, Wikipedia
See also: besselhx for an exponentially scaled variant.
SpecialFunctions.besselhx — Functionbesselhx(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
SpecialFunctions.besseli — Methodbesseli(nu, x)Modified Bessel function of the first kind of order nu, $I_\nu(x)$.
External links: DLMF, Wikipedia
See also: besselix(nu,x), besselj(nu,x), besselk(nu,x)
SpecialFunctions.besselix — Methodbesselix(nu, x)Scaled modified Bessel function of the first kind of order nu, $I_\nu(x) e^{- | \operatorname{Re}(x) |}$.
External links: DLMF, Wikipedia
See also: besseli(nu,x), besselj(nu,x), besselk(nu,x)
SpecialFunctions.besselj — Methodbesselj(nu, x)Bessel function of the first kind of order nu, $J_\nu(x)$.
External links: DLMF, Wikipedia
See also: besseljx(nu,x), besseli(nu,x), besselk(nu,x)
SpecialFunctions.besselj0 — Methodbesselj0(x)Bessel function of the first kind of order 0, $J_0(x)$.
External links: DLMF, Wikipedia
See also: besselj(nu,x)
SpecialFunctions.besselj1 — Methodbesselj1(x)Bessel function of the first kind of order 1, $J_1(x)$.
External links: DLMF, Wikipedia
See also: besselj(nu,x)
SpecialFunctions.besseljx — Methodbesseljx(nu, x)Scaled Bessel function of the first kind of order nu, $J_\nu(x) e^{- | \operatorname{Im}(x) |}$.
External links: DLMF, Wikipedia
See also: besselj(nu,x), besseli(nu,x), besselk(nu,x)
SpecialFunctions.besselk — Methodbesselk(nu, x)Modified Bessel function of the second kind of order nu, $K_\nu(x)$.
External links: DLMF, Wikipedia
See also: See also: besselkx(nu,x), besseli(nu,x), besselj(nu,x)
SpecialFunctions.besselkx — Methodbesselkx(nu, x)Scaled modified Bessel function of the second kind of order nu, $K_\nu(x) e^x$.
External links: DLMF, Wikipedia
See also: besselk(nu,x), besseli(nu,x), besselj(nu,x)
SpecialFunctions.bessely — Methodbessely(nu, x)Bessel function of the second kind of order nu, $Y_\nu(x)$.
External links: DLMF, Wikipedia
See also besselyx(nu,x) for a scaled variant.
SpecialFunctions.bessely0 — Methodbessely0(x)Bessel function of the second kind of order 0, $Y_0(x)$.
External links: DLMF, Wikipedia
See also: bessely(nu,x)
SpecialFunctions.bessely1 — Methodbessely1(x)Bessel function of the second kind of order 1, $Y_1(x)$.
External links: DLMF, Wikipedia
See also: bessely(nu,x)
SpecialFunctions.besselyx — Methodbesselyx(nu, x)Scaled Bessel function of the second kind of order nu, $Y_\nu(x) e^{- | \operatorname{Im}(x) |}$.
External links: DLMF, Wikipedia
See also bessely(nu,x)
SpecialFunctions.beta — Methodbeta(x, y)Euler integral of the first kind $\operatorname{B}(x,y) = \Gamma(x)\Gamma(y)/\Gamma(x+y)$.
SpecialFunctions.beta_inc — Methodbeta_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
SpecialFunctions.beta_inc_asymptotic_asymmetric — Methodbeta_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, Wikipedia
See also: beta_inc
SpecialFunctions.beta_inc_asymptotic_symmetric — Methodbeta_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, Wikipedia
See also: beta_inc
Implementation
BASYM(A,B,LAMBDA,EPS) from Didonato and Morris (1982)
SpecialFunctions.beta_inc_cont_fraction — Methodbeta_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)
SpecialFunctions.beta_inc_diff — Methodbeta_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.
External links: DLMF, Wikipedia
See also: beta_inc
SpecialFunctions.beta_inc_inv — Methodbeta_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
SpecialFunctions.beta_inc_power_series — Methodbeta_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)
SpecialFunctions.beta_inc_power_series1 — Methodbeta_inc_power_series1(a,b,x,epps)Another variant of BPSER(A,B,X,EPS).
External links: DLMF, Wikipedia
See also: beta_inc
Implementation
APSER(A,B,X,EPS) from Didonato and Morris (1982)
SpecialFunctions.beta_inc_power_series2 — Methodbeta_inc_power_series2(a,b,x,epps)Variant of BPSER(A,B,X,EPS).
External links: DLMF, Wikipedia
See also: beta_inc
Implementation
FPSER(A,B,X,EPS) from Didonato and Morris (1982)
SpecialFunctions.beta_integrand — Functionbeta_integrand(a, b, x, y, mu=0.0)Compute $e^{mu} * x^{a}y^{b}/B(a,b)$
SpecialFunctions.chepolsum — Methodchepolsum(n,x,a)Computes a series of Chebyshev Polynomials given by: a[1]/2 + a[2]T1(x) + .... + a[n]T{n-1}(X).
SpecialFunctions.coeff1 — MethodComputing 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
SpecialFunctions.coeff2 — MethodComputing 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
SpecialFunctions.coeff3 — MethodComputing 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
SpecialFunctions.cosint — Functioncosint(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$.
SpecialFunctions.dawson — Functiondawson(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.
SpecialFunctions.digamma — Methoddigamma(x)Compute the digamma function of x (the logarithmic derivative of gamma(x)).
SpecialFunctions.ellipe — Methodellipe(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]$.
SpecialFunctions.ellipk — Methodellipk(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]$.
SpecialFunctions.erf — Functionerf(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
SpecialFunctions.erfc — Functionerfc(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
SpecialFunctions.erfcinv — Methoderfcinv(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.
SpecialFunctions.erfcx — Functionerfcx(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.
SpecialFunctions.erfi — Functionerfi(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.
SpecialFunctions.erfinv — Methoderfinv(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.
SpecialFunctions.esum — Methodesum(mu,x)Compute $e^{mu+x}$
SpecialFunctions.eta — Methodeta(s)Dirichlet eta function $\eta(s) = \sum^\infty_{n=1}(-1)^{n-1}/n^{s}$.
SpecialFunctions.expint — Functionexpint(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.
SpecialFunctions.expinti — Methodexpinti(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.
SpecialFunctions.expintx — Functionexpintx(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)
SpecialFunctions.faddeeva — Functionfaddeeva(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)$.
SpecialFunctions.gamma — Methodgamma(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)$.
SpecialFunctions.gamma — Methodgamma(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
SpecialFunctions.gamma_inc — Functiongamma_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)
SpecialFunctions.gamma_inc_asym — Methodgamma_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)
SpecialFunctions.gamma_inc_cf — Methodgamma_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)
SpecialFunctions.gamma_inc_fsum — Methodgamma_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)
SpecialFunctions.gamma_inc_inv — Methodgamma_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, Wikipedia
See also: gamma_inc(a,x,ind).
SpecialFunctions.gamma_inc_inv_alarge — Methodgamma_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.
SpecialFunctions.gamma_inc_inv_psmall — Methodgamma_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}$.
SpecialFunctions.gamma_inc_inv_qsmall — Methodgamma_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}$.
SpecialFunctions.gamma_inc_minimax — Methodgamma_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)
SpecialFunctions.gamma_inc_taylor — Methodgamma_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
See also: gamma_inc(a,x,ind)
SpecialFunctions.gamma_inc_taylor_x — Methodgamma_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)
SpecialFunctions.gamma_inc_temme — Methodgamma_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)
SpecialFunctions.gamma_inc_temme_1 — Methodgamma_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
SpecialFunctions.gammax — Methodgammax(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}\]
SpecialFunctions.hankelh1 — Methodhankelh1(nu, x)Bessel function of the third kind of order nu, $H^{(1)}_\nu(x)$.
External links: DLMF, Wikipedia
See also: hankelh1x
SpecialFunctions.hankelh1x — Methodhankelh1x(nu, x)Scaled Bessel function of the third kind of order nu, $H^{(1)}_\nu(x) e^{-x i}$.
External links: DLMF, Wikipedia
See also: hankelh1
SpecialFunctions.hankelh2 — Methodhankelh2(nu, x)Bessel function of the third kind of order nu, $H^{(2)}_\nu(x)$.
External links: DLMF, Wikipedia
See also: hankelh2x(nu,x)
SpecialFunctions.hankelh2x — Methodhankelh2x(nu, x)Scaled Bessel function of the third kind of order nu, $H^{(2)}_\nu(x) e^{x i}$.
External links: DLMF, Wikipedia
See also: hankelh2(nu,x)
SpecialFunctions.invdigamma — Methodinvdigamma(x)Compute the inverse digamma function of x.
SpecialFunctions.jinc — Methodjinc(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
SpecialFunctions.lambdaeta — Methodlambdaeta(eta)Compute the value of $\lambda$ satisfying $\eta^{2}/2 = \lambda-1-\log{\lambda}$.
SpecialFunctions.logabsbeta — Methodlogabsbeta(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.
SpecialFunctions.logabsbinomial — Methodlogabsbinomial(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))).
SpecialFunctions.logabsgamma — Methodlogabsgamma(x)Compute the logarithm of absolute value of gamma for Real xand returns a tuple (log(abs(gamma(x))), sign(gamma(x))).
See also loggamma.
SpecialFunctions.logbeta — Methodlogbeta(x, y)Natural logarithm of the beta function $\log(|\operatorname{B}(x,y)|)$.
See also logabsbeta.
SpecialFunctions.logerf — Methodlogerf(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).
SpecialFunctions.logerfc — Methodlogerfc(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.
SpecialFunctions.logerfcx — Methodlogerfcx(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.
SpecialFunctions.logfactorial — Methodlogfactorial(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.
SpecialFunctions.loggamma — Methodloggamma(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).
SpecialFunctions.loggamma — Methodloggamma(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.
SpecialFunctions.loggamma1p — Methodloggamma1p(x)Compute $log(\Gamma(1+x))$ for -1 < x <= 1.
SpecialFunctions.loggammadiv — Methodloggammadiv(a,b)Computes $log(\Gamma(b)/\Gamma(a+b))$ when b >= 8
SpecialFunctions.ncF — MethodncF(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
SpecialFunctions.ncbeta — Methodncbeta(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.
SpecialFunctions.ncbeta_poisson — Methodncbeta_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}\]
SpecialFunctions.ncbeta_tail — Methodncbeta_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.
SpecialFunctions.polygamma — Methodpolygamma(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)
SpecialFunctions.rgamma1pm1 — Methodrgamma1pm1(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).
SpecialFunctions.rgammax — Methodrgammax(a,x)Evaluation of $1/\Gamma(a) e^{-x} x^{a}$. Based on DRCOMP from the NSWC library.
SpecialFunctions.sinint — Functionsinint(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$.
SpecialFunctions.sphericalbesselj — Methodsphericalbesselj(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.
SpecialFunctions.sphericalbessely — Methodsphericalbessely(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.
SpecialFunctions.stirling_corr — Methodstirling_corr(a0,b0)Compute stirling(a0) + stirling(b0) - stirling(a0 + b0) for a0, b0 >= 8
SpecialFunctions.stirling_error — Methodstirling_error(x)Compute $\ln{\Gamma(x)} - (x-0.5)*\ln{x} + x - \ln{(2\pi)}/2$. Adapted from stirling in IncgamFI.
SpecialFunctions.trigamma — Methodtrigamma(x)Compute the trigamma function of x (the logarithmic second derivative of gamma(x)).
SpecialFunctions.zeta — Methodzeta(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
SpecialFunctions.zeta — Methodzeta(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