diff --git a/docs/make.jl b/docs/make.jl index 62c5f6c1..ff36a922 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -8,6 +8,8 @@ makedocs(modules=[SpecialFunctions], format = Documenter.HTML(; assets = String[]), pages=["Home" => "index.md", "Overview" => "functions_overview.md", - "Reference" => "functions_list.md"]) + "Reference" => "functions_list.md"], + #warnonly=[:missing_docs], + ) deploydocs(repo="github.com/JuliaMath/SpecialFunctions.jl.git") diff --git a/docs/src/functions_list.md b/docs/src/functions_list.md index da32182e..dc595591 100644 --- a/docs/src/functions_list.md +++ b/docs/src/functions_list.md @@ -4,7 +4,150 @@ CurrentModule = SpecialFunctions ``` -```@autodocs -Modules = [SpecialFunctions] -Order = [:module, :type, :function] + +## Gamma Function + +```@docs +gamma +loggamma +logabsgamma +loggamma1p +logfactorial +digamma +invdigamma +trigamma +polygamma +gamma_inc +gamma_inc_inv +loggammadiv +gammax +rgammax +rgamma1pm1 +gamma_inc_asym +gamma_inc_cf +gamma_inc_fsum +gamma_inc_minimax +gamma_inc_taylor +gamma_inc_taylor_x +gamma_inc_temme +gamma_inc_temme_1 +gamma_inc_inv_psmall +gamma_inc_inv_qsmall +gamma_inc_inv_alarge +auxgam +``` + +### Beta Function +```@docs +beta +logbeta +logabsbeta +logabsbinomial +beta_inc +beta_inc_inv +ncbeta +ncbeta_poisson +ncbeta_tail +beta_inc_power_series1 +beta_inc_power_series2 +beta_inc_asymptotic_asymmetric +beta_inc_asymptotic_symmetric +beta_inc_power_series +beta_inc_diff +beta_inc_cont_fraction +beta_integrand +``` + +### Utilities + +```@docs +chepolsum +lambdaeta +stirling_corr +stirling_error +esum +coeff1 +coeff2 +coeff3 +ncF +``` + +## Exponential and Trigonometric Integrals + +```@docs +expint +expinti +expintx +sinint +cosint +``` + +## Error Functions, Dawson’s and Fresnel Integrals + +```@docs +erf +erf(::Real, ::Real) +erfc +logerf +erfcinv +erfcx +logerfc +logerfcx +erfi +erfinv +dawson +faddeeva +``` + +## Airy and Related Functions + +```@docs +airyai +airyaiprime +airybi +airybiprime +airyaix +airyaiprimex +airybix +airybiprimex +``` + +## Bessel Functions + +```@docs +besselj +besselj0 +besselj1 +besseljx +sphericalbesselj +bessely +bessely0 +bessely1 +besselyx +sphericalbessely +besselh +besselhx +hankelh1 +hankelh1x +hankelh2 +hankelh2x +besseli +besselix +besselk +besselkx +jinc +``` + +## Elliptic Integrals + +```@docs +ellipk +ellipe +``` + +## Zeta and Related Functions + +```@docs +eta +zeta ``` diff --git a/docs/src/functions_overview.md b/docs/src/functions_overview.md index d9f5bd1f..ca22f5ac 100644 --- a/docs/src/functions_overview.md +++ b/docs/src/functions_overview.md @@ -1,102 +1,120 @@ # Functions -Here the *Special Functions* are listed according to the structure of [NIST Digital Library of Mathematical Functions](https://dlmf.nist.gov/). +Here the *Special Functions* are listed according to the structure of [NIST Digital Library of Mathematical Functions (DLMF)](https://dlmf.nist.gov/). ## Gamma Function -- [Gamma Function - DLMF](https://dlmf.nist.gov/5) +[Gamma Function - DLMF](https://dlmf.nist.gov/5) | Function | Description | |:-------- |:----------- | -| [`gamma(z)`](@ref SpecialFunctions.gamma(::Number)) | [gamma function](https://en.wikipedia.org/wiki/Gamma_function) ``\Gamma(z)`` | -| [`loggamma(x)`](@ref SpecialFunctions.loggamma(::Number)) | accurate `log(gamma(x))` for large `x` | -| [`logabsgamma(x)`](@ref SpecialFunctions.logabsgamma) | accurate `log(abs(gamma(x)))` for large `x` | -| [`logfactorial(x)`](@ref SpecialFunctions.logfactorial) | accurate `log(factorial(x))` for large `x`; same as `loggamma(x+1)` for `x > 1`, zero otherwise | -| [`digamma(x)`](@ref SpecialFunctions.digamma) | [digamma function](https://en.wikipedia.org/wiki/Digamma_function) (i.e. the derivative of `loggamma` at `x`) | -| [`invdigamma(x)`](@ref SpecialFunctions.invdigamma) | [invdigamma function](http://bariskurt.com/calculating-the-inverse-of-digamma-function/) (i.e. inverse of `digamma` function at `x` using fixed-point iteration algorithm) | -| [`trigamma(x)`](@ref SpecialFunctions.trigamma) | [trigamma function](https://en.wikipedia.org/wiki/Trigamma_function) (i.e the logarithmic second derivative of `gamma` at `x`) | -| [`polygamma(m,x)`](@ref SpecialFunctions.polygamma) | [polygamma function](https://en.wikipedia.org/wiki/Polygamma_function) (i.e the (m+1)-th derivative of the `loggamma` function at `x`) | -| [`gamma(a,z)`](@ref SpecialFunctions.gamma(::Number,::Number)) | [upper incomplete gamma function ``\Gamma(a,z)``](https://en.wikipedia.org/wiki/Incomplete_gamma_function) | -| [`loggamma(a,z)`](@ref SpecialFunctions.loggamma(::Number,::Number)) | accurate `log(gamma(a,x))` for large arguments | -| [`gamma_inc(a,x,IND)`](@ref SpecialFunctions.gamma_inc) | [incomplete gamma function ratio P(a,x) and Q(a,x)](https://en.wikipedia.org/wiki/Incomplete_gamma_function) (i.e evaluates P(a,x) and Q(a,x)for accuracy specified by IND and returns tuple (p,q)) | -| [`gamma_inc_inv(a,p,q)`](@ref SpecialFunctions.gamma_inc_inv) | [inverse of incomplete gamma function ratio P(a,x) and Q(a,x)](https://en.wikipedia.org/wiki/Incomplete_gamma_function) (i.e evaluates x given P(a,x)=p and Q(a,x)=q | -| [`beta(x,y)`](@ref SpecialFunctions.beta) | [beta function](https://en.wikipedia.org/wiki/Beta_function) at `x,y` | -| [`logbeta(x,y)`](@ref SpecialFunctions.logbeta) | accurate `log(beta(x,y))` for large `x` or `y` | -| [`logabsbeta(x,y)`](@ref SpecialFunctions.logabsbeta) | accurate `log(abs(beta(x,y)))` for large `x` or `y` | -| [`logabsbinomial(x,y)`](@ref SpecialFunctions.logabsbinomial) | accurate `log(abs(binomial(n,k)))` for large `n` and `k` near `n/2` | -| [`beta_inc(a,b,x,y)`](@ref SpecialFunctions.beta_inc) | [incomplete beta function ratio Ix(a,b) and Iy(a,b)](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) (i.e evaluates Ix(a,b) and Iy(a,b) and returns tuple (p,q)) | -| [`beta_inc_inv(a,b,p,q)`](@ref SpecialFunctions.beta_inc_inv) | Inverse of the incomplete beta function (i.e evaluates x given ``I_{x}(a, b) = p``) | - - -## [Exponential and Trigonometric Integrals](https://dlmf.nist.gov/6) +| [`gamma(z)`](@ref SpecialFunctions.gamma(::Number)) | [gamma function](https://en.wikipedia.org/wiki/Gamma_function) ``\Gamma(z)`` | +| [`loggamma(x)`](@ref SpecialFunctions.loggamma(::Number)) | accurate `log(gamma(x))` for large `x` | +| [`logabsgamma(x)`](@ref SpecialFunctions.logabsgamma) | accurate `log(abs(gamma(x)))` for large `x` | +| [`logfactorial(x)`](@ref SpecialFunctions.logfactorial) | accurate `log(factorial(x))` for large `x`; same as `loggamma(x+1)` for `x > 1`, zero otherwise | +| [`digamma(x)`](@ref SpecialFunctions.digamma) | [digamma function](https://en.wikipedia.org/wiki/Digamma_function) (i.e. the derivative of `loggamma` at `x`) | +| [`invdigamma(x)`](@ref SpecialFunctions.invdigamma) | [invdigamma function](http://bariskurt.com/calculating-the-inverse-of-digamma-function/) (i.e. inverse of `digamma` function at `x` using fixed-point iteration algorithm) | +| [`trigamma(x)`](@ref SpecialFunctions.trigamma) | [trigamma function](https://en.wikipedia.org/wiki/Trigamma_function) (i.e the logarithmic second derivative of `gamma` at `x`) | +| [`polygamma(m,x)`](@ref SpecialFunctions.polygamma) | [polygamma function](https://en.wikipedia.org/wiki/Polygamma_function) (i.e the (m+1)-th derivative of the `loggamma` function at `x`) | +| [`gamma(a,z)`](@ref SpecialFunctions.gamma(::Number,::Number)) | [upper incomplete gamma function ``\Gamma(a,z)``](https://en.wikipedia.org/wiki/Incomplete_gamma_function) | +| [`loggamma(a,z)`](@ref SpecialFunctions.loggamma(::Number,::Number)) | accurate `log(gamma(a,x))` for large arguments | +| [`gamma_inc(a,x,IND)`](@ref SpecialFunctions.gamma_inc) | [incomplete gamma function ratio P(a,x) and Q(a,x)](https://en.wikipedia.org/wiki/Incomplete_gamma_function) (i.e evaluates P(a,x) and Q(a,x) for accuracy specified by IND and returns tuple (p,q)) | +| [`gamma_inc_inv(a,p,q)`](@ref SpecialFunctions.gamma_inc_inv) | [inverse of incomplete gamma function ratio P(a,x) and Q(a,x)](https://en.wikipedia.org/wiki/Incomplete_gamma_function) (i.e evaluates x given P(a,x)=p and Q(a,x)=q) | +| [`beta(x,y)`](@ref SpecialFunctions.beta) | [beta function](https://en.wikipedia.org/wiki/Beta_function) at `x,y` | +| [`logbeta(x,y)`](@ref SpecialFunctions.logbeta) | accurate `log(beta(x,y))` for large `x` or `y` | +| [`logabsbeta(x,y)`](@ref SpecialFunctions.logabsbeta) | accurate `log(abs(beta(x,y)))` for large `x` or `y` | +| [`logabsbinomial(x,y)`](@ref SpecialFunctions.logabsbinomial) | accurate `log(abs(binomial(n,k)))` for large `n` and `k` near `n/2` | +| [`beta_inc(a,b,x,y)`](@ref SpecialFunctions.beta_inc) | [incomplete beta function ratio Ix(a,b) and Iy(a,b)](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) (i.e evaluates Ix(a,b) and Iy(a,b) and returns tuple (p,q)) | +| [`beta_inc_inv(a,b,p,q)`](@ref SpecialFunctions.beta_inc_inv) | Inverse of the incomplete beta function (i.e evaluates x given ``I_{x}(a, b) = p``) | + + +## Exponential and Trigonometric Integrals + +[Exponential and Trigonometric Integrals - DLMF](https://dlmf.nist.gov/6) + | Function | Description | |:-------- |:----------- | -| [`expint(ν, z)`](@ref SpecialFunctions.expint) | [exponential integral](https://en.wikipedia.org/wiki/Exponential_integral) ``\operatorname{E}_\nu(z)`` | -| [`expinti(x)`](@ref SpecialFunctions.expinti) | [exponential integral](https://en.wikipedia.org/wiki/Exponential_integral) ``\operatorname{Ei}(x)`` | -| [`expintx(x)`](@ref SpecialFunctions.expintx) | scaled [exponential integral](https://en.wikipedia.org/wiki/Exponential_integral) ``e^z \operatorname{E}_\nu(z)`` | -| [`sinint(x)`](@ref SpecialFunctions.sinint) | [sine integral](https://en.wikipedia.org/wiki/Trigonometric_integral#Sine_integral) ``\operatorname{Si}(x)`` | -| [`cosint(x)`](@ref SpecialFunctions.cosint) | [cosine integral](https://en.wikipedia.org/wiki/Trigonometric_integral#Cosine_integral) ``\operatorname{Ci}(x)`` | +| [`expint(ν, z)`](@ref SpecialFunctions.expint) | [exponential integral](https://en.wikipedia.org/wiki/Exponential_integral) ``\operatorname{E}_\nu(z)`` | +| [`expinti(x)`](@ref SpecialFunctions.expinti) | [exponential integral](https://en.wikipedia.org/wiki/Exponential_integral) ``\operatorname{Ei}(x)`` | +| [`expintx(x)`](@ref SpecialFunctions.expintx) | scaled [exponential integral](https://en.wikipedia.org/wiki/Exponential_integral) ``e^z \operatorname{E}_\nu(z)`` | +| [`sinint(x)`](@ref SpecialFunctions.sinint) | [sine integral](https://en.wikipedia.org/wiki/Trigonometric_integral#Sine_integral) ``\operatorname{Si}(x)`` | +| [`cosint(x)`](@ref SpecialFunctions.cosint) | [cosine integral](https://en.wikipedia.org/wiki/Trigonometric_integral#Cosine_integral) ``\operatorname{Ci}(x)`` | -## [Error Functions, Dawson’s and Fresnel Integrals](https://dlmf.nist.gov/7) +## Error Functions, Dawson’s and Fresnel Integrals + +[Error Functions, Dawson’s and Fresnel Integrals - DLMF](https://dlmf.nist.gov/7) + | Function | Description | |:-------- |:----------- | -| [`erf(x)`](@ref SpecialFunctions.erf) | [error function](https://en.wikipedia.org/wiki/Error_function) at ``x`` | -| [`erf(x,y)`](@ref SpecialFunctions.erf) | accurate version of ``\operatorname{erf}(y) - \operatorname{erf}(x)`` | -| [`erfc(x)`](@ref SpecialFunctions.erfc) | complementary error function, i.e. the accurate version of ``1-\operatorname{erf}(x)`` for large ``x`` | -| [`erfcinv(x)`](@ref SpecialFunctions.erfcinv) | inverse function to [`erfc()`](@ref SpecialFunctions.erfc) | -| [`erfcx(x)`](@ref SpecialFunctions.erfcx) | scaled complementary error function, i.e. accurate ``e^{x^2} \operatorname{erfc}(x)`` for large ``x`` | -| [`logerfc(x)`](@ref SpecialFunctions.logerfc) | log of the complementary error function, i.e. accurate ``\operatorname{ln}(\operatorname{erfc}(x))`` for large ``x`` | +| [`erf(x)`](@ref SpecialFunctions.erf) | [error function](https://en.wikipedia.org/wiki/Error_function) at ``x`` | +| [`erf(x,y)`](@ref SpecialFunctions.erf) | accurate version of ``\operatorname{erf}(y) - \operatorname{erf}(x)`` | +| [`erfc(x)`](@ref SpecialFunctions.erfc) | complementary error function, i.e. the accurate version of ``1-\operatorname{erf}(x)`` for large ``x`` | +| [`erfcinv(x)`](@ref SpecialFunctions.erfcinv) | inverse function to [`erfc()`](@ref SpecialFunctions.erfc) | +| [`erfcx(x)`](@ref SpecialFunctions.erfcx) | scaled complementary error function, i.e. accurate ``e^{x^2} \operatorname{erfc}(x)`` for large ``x`` | +| [`logerfc(x)`](@ref SpecialFunctions.logerfc) | log of the complementary error function, i.e. accurate ``\operatorname{ln}(\operatorname{erfc}(x))`` for large ``x`` | | [`logerfcx(x)`](@ref SpecialFunctions.logerfcx) | log of the scaled complementary error function, i.e. accurate ``\operatorname{ln}(\operatorname{erfcx}(x))`` for large negative ``x`` | -| [`erfi(x)`](@ref SpecialFunctions.erfi) | imaginary error function defined as ``-i \operatorname{erf}(ix)`` | -| [`erfinv(x)`](@ref SpecialFunctions.erfinv) | inverse function to [`erf()`](@ref SpecialFunctions.erf) | -| [`dawson(x)`](@ref SpecialFunctions.dawson) | scaled imaginary error function, a.k.a. Dawson function, i.e. accurate ``\frac{\sqrt{\pi}}{2} e^{-x^2} \operatorname{erfi}(x)`` for large ``x`` | +| [`erfi(x)`](@ref SpecialFunctions.erfi) | imaginary error function defined as ``-i \operatorname{erf}(ix)`` | +| [`erfinv(x)`](@ref SpecialFunctions.erfinv) | inverse function to [`erf()`](@ref SpecialFunctions.erf) | +| [`dawson(x)`](@ref SpecialFunctions.dawson) | scaled imaginary error function, a.k.a. Dawson function, i.e. accurate ``\frac{\sqrt{\pi}}{2} e^{-x^2} \operatorname{erfi}(x)`` for large ``x`` | | [`faddeeva(x)`](@ref SpecialFunctions.faddeeva) | [Faddeeva function](https://en.wikipedia.org/wiki/Faddeeva_function), equivalent to ``\operatorname{erfcx}(-ix)`` | -## [Airy and Related Functions](https://dlmf.nist.gov/9) +## Airy and Related Functions + +[Airy and Related Functions - DLMF](https://dlmf.nist.gov/9) + | Function | Description | |:-------- |:----------- | -| [`airyai(z)`](@ref SpecialFunctions.airyai) | [Airy Ai function](https://en.wikipedia.org/wiki/Airy_function) at `z` | -| [`airyaiprime(z)`](@ref SpecialFunctions.airyaiprime) | derivative of the Airy Ai function at `z` | -| [`airybi(z)`](@ref SpecialFunctions.airybi) | [Airy Bi function](https://en.wikipedia.org/wiki/Airy_function) at `z` | -| [`airybiprime(z)`](@ref SpecialFunctions.airybiprime) | derivative of the Airy Bi function at `z` | -| [`airyaix(z)`](@ref SpecialFunctions.airyaix), [`airyaiprimex(z)`](@ref SpecialFunctions.airyaiprimex), [`airybix(z)`](@ref SpecialFunctions.airybix), [`airybiprimex(z)`](@ref SpecialFunctions.airybiprimex) | scaled Airy Ai function and `k`th derivatives at `z` | +| [`airyai(z)`](@ref SpecialFunctions.airyai) | [Airy Ai function](https://en.wikipedia.org/wiki/Airy_function) at `z` | +| [`airyaiprime(z)`](@ref SpecialFunctions.airyaiprime) | derivative of the Airy Ai function at `z` | +| [`airybi(z)`](@ref SpecialFunctions.airybi) | [Airy Bi function](https://en.wikipedia.org/wiki/Airy_function) at `z` | +| [`airybiprime(z)`](@ref SpecialFunctions.airybiprime) | derivative of the Airy Bi function at `z` | +| [`airyaix(z)`](@ref SpecialFunctions.airyaix), [`airyaiprimex(z)`](@ref SpecialFunctions.airyaiprimex), [`airybix(z)`](@ref SpecialFunctions.airybix), [`airybiprimex(z)`](@ref SpecialFunctions.airybiprimex) | scaled Airy Ai function and derivative at `z` | + +## Bessel Functions + +[Bessel Functions - DLMF](https://dlmf.nist.gov/10) -## [Bessel Functions](https://dlmf.nist.gov/10) | Function | Description | |:-------- |:----------- | -| [`besselj(nu,z)`](@ref SpecialFunctions.besselj) | [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the first kind of order `nu` at `z` | -| [`besselj0(z)`](@ref SpecialFunctions.besselj0) | `besselj(0,z)` | -| [`besselj1(z)`](@ref SpecialFunctions.besselj1) | `besselj(1,z)` | -| [`besseljx(nu,z)`](@ref SpecialFunctions.besseljx) | scaled Bessel function of the first kind of order `nu` at `z` | -| [`sphericalbesselj(nu,z)`](@ref SpecialFunctions.sphericalbesselj) | [Spherical Bessel function](https://en.wikipedia.org/wiki/Bessel_function#Spherical_Bessel_functions:_jn,_yn) of the first kind of order `nu` at `z` | -| [`bessely(nu,z)`](@ref SpecialFunctions.bessely) | [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the second kind of order `nu` at `z` | -| [`bessely0(z)`](@ref SpecialFunctions.bessely0) | `bessely(0,z)` | -| [`bessely1(z)`](@ref SpecialFunctions.bessely1) | `bessely(1,z)` | -| [`besselyx(nu,z)`](@ref SpecialFunctions.besselyx) | scaled Bessel function of the second kind of order `nu` at `z` | -| [`sphericalbessely(nu,z)`](@ref SpecialFunctions.sphericalbessely) | [Spherical Bessel function](https://en.wikipedia.org/wiki/Bessel_function#Spherical_Bessel_functions:_jn,_yn) of the second kind of order `nu` at `z` | -| [`besselh(nu,k,z)`](@ref SpecialFunctions.besselh) | [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the third kind (a.k.a. Hankel function) of order `nu` at `z`; `k` must be either `1` or `2` | -| [`hankelh1(nu,z)`](@ref SpecialFunctions.hankelh1) | `besselh(nu, 1, z)` | -| [`hankelh1x(nu,z)`](@ref SpecialFunctions.hankelh1x) | scaled `besselh(nu, 1, z)` | -| [`hankelh2(nu,z)`](@ref SpecialFunctions.hankelh2) | `besselh(nu, 2, z)` | -| [`hankelh2x(nu,z)`](@ref SpecialFunctions.hankelh2x) | scaled `besselh(nu, 2, z)` | -| [`besseli(nu,z)`](@ref SpecialFunctions.besseli) | modified [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the first kind of order `nu` at `z` | -| [`besselix(nu,z)`](@ref SpecialFunctions.besselix) | scaled modified Bessel function of the first kind of order `nu` at `z` | -| [`besselk(nu,z)`](@ref SpecialFunctions.besselk) | modified [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the second kind of order `nu` at `z` | -| [`besselkx(nu,z)`](@ref SpecialFunctions.besselkx) | scaled modified Bessel function of the second kind of order `nu` at `z` | -| [`jinc(x)`](@ref SpecialFunctions.jinc) | scaled [Bessel function of the first kind divided by `x`](https://en.wikipedia.org/wiki/Sombrero_function). A.k.a. sombrero or besinc | - - -## [Elliptic Integrals](https://dlmf.nist.gov/19) +| [`besselj(nu,z)`](@ref SpecialFunctions.besselj) | [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the first kind of order `nu` at `z` | +| [`besselj0(z)`](@ref SpecialFunctions.besselj0) | `besselj(0,z)` | +| [`besselj1(z)`](@ref SpecialFunctions.besselj1) | `besselj(1,z)` | +| [`besseljx(nu,z)`](@ref SpecialFunctions.besseljx) | scaled Bessel function of the first kind of order `nu` at `z` | +| [`sphericalbesselj(nu,z)`](@ref SpecialFunctions.sphericalbesselj) | [Spherical Bessel function](https://en.wikipedia.org/wiki/Bessel_function#Spherical_Bessel_functions:_jn,_yn) of the first kind of order `nu` at `z` | +| [`bessely(nu,z)`](@ref SpecialFunctions.bessely) | [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the second kind of order `nu` at `z` | +| [`bessely0(z)`](@ref SpecialFunctions.bessely0) | `bessely(0,z)` | +| [`bessely1(z)`](@ref SpecialFunctions.bessely1) | `bessely(1,z)` | +| [`besselyx(nu,z)`](@ref SpecialFunctions.besselyx) | scaled Bessel function of the second kind of order `nu` at `z` | +| [`sphericalbessely(nu,z)`](@ref SpecialFunctions.sphericalbessely) | [Spherical Bessel function](https://en.wikipedia.org/wiki/Bessel_function#Spherical_Bessel_functions:_jn,_yn) of the second kind of order `nu` at `z` | +| [`besselh(nu,k,z)`](@ref SpecialFunctions.besselh) | [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the third kind (a.k.a. Hankel function) of order `nu` at `z`; `k` must be either `1` or `2` | +| [`hankelh1(nu,z)`](@ref SpecialFunctions.hankelh1) | `besselh(nu, 1, z)` | +| [`hankelh1x(nu,z)`](@ref SpecialFunctions.hankelh1x) | scaled `besselh(nu, 1, z)` | +| [`hankelh2(nu,z)`](@ref SpecialFunctions.hankelh2) | `besselh(nu, 2, z)` | +| [`hankelh2x(nu,z)`](@ref SpecialFunctions.hankelh2x) | scaled `besselh(nu, 2, z)` | +| [`besseli(nu,z)`](@ref SpecialFunctions.besseli) | modified [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the first kind of order `nu` at `z` | +| [`besselix(nu,z)`](@ref SpecialFunctions.besselix) | scaled modified Bessel function of the first kind of order `nu` at `z` | +| [`besselk(nu,z)`](@ref SpecialFunctions.besselk) | modified [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the second kind of order `nu` at `z` | +| [`besselkx(nu,z)`](@ref SpecialFunctions.besselkx) | scaled modified Bessel function of the second kind of order `nu` at `z` | +| [`jinc(x)`](@ref SpecialFunctions.jinc) | scaled [Bessel function of the first kind divided by `x`](https://en.wikipedia.org/wiki/Sombrero_function). A.k.a. sombrero or besinc | + + +## Elliptic Integrals + +[Elliptic Integrals - DLMF](https://dlmf.nist.gov/19) + | Function | Description | |:-------- |:----------- | -| [`ellipk(m)`](@ref SpecialFunctions.ellipk) | [complete elliptic integral of 1st kind](https://en.wikipedia.org/wiki/Elliptic_integral#Notational_variants) ``K(m)`` | -| [`ellipe(m)`](@ref SpecialFunctions.ellipe) | [complete elliptic integral of 2nd kind](https://en.wikipedia.org/wiki/Elliptic_integral#Complete_elliptic_integral_of_the_second_kind) ``E(m)`` | +| [`ellipk(m)`](@ref SpecialFunctions.ellipk) | [complete elliptic integral of 1st kind](https://en.wikipedia.org/wiki/Elliptic_integral#Notational_variants) ``K(m)`` | +| [`ellipe(m)`](@ref SpecialFunctions.ellipe) | [complete elliptic integral of 2nd kind](https://en.wikipedia.org/wiki/Elliptic_integral#Complete_elliptic_integral_of_the_second_kind) ``E(m)`` | + + +## Zeta and Related Functions +[Zeta and Related Functions - DLMF](https://dlmf.nist.gov/25) -## [Zeta and Related Functions](https://dlmf.nist.gov/25) | Function | Description | |:-------- |:----------- | -| [`eta(x)`](@ref SpecialFunctions.eta) | [Dirichlet eta function](https://en.wikipedia.org/wiki/Dirichlet_eta_function) at `x` | -| [`zeta(x)`](@ref SpecialFunctions.zeta) | [Riemann zeta function](https://en.wikipedia.org/wiki/Riemann_zeta_function) at `x` | +| [`eta(x)`](@ref SpecialFunctions.eta) | [Dirichlet eta function](https://en.wikipedia.org/wiki/Dirichlet_eta_function) at `x` | +| [`zeta(x)`](@ref SpecialFunctions.zeta) | [Riemann zeta function](https://en.wikipedia.org/wiki/Riemann_zeta_function) at `x` | diff --git a/docs/src/index.md b/docs/src/index.md index 89be4bcd..a31c771d 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -15,5 +15,5 @@ The latest version of the package is available for Julia versions 1.3 and up. To install it, run the following at the Julia REPL: ```julia -Pkg.add("SpecialFunctions") +import Pkg; Pkg.add("SpecialFunctions") ``` diff --git a/src/bessel.jl b/src/bessel.jl index e398b9da..a068c722 100644 --- a/src/bessel.jl +++ b/src/bessel.jl @@ -61,98 +61,118 @@ function _biry(z::Complex{Float64}, id::Int32, kode::Int32) end -""" +@doc raw""" airyai(x) -Airy function of the first kind ``\\operatorname{Ai}(x)``. +Airy function of the first kind ``\operatorname{Ai}(x)``. -External links: [DLMF](https://dlmf.nist.gov/9.2), [Wikipedia](https://en.wikipedia.org/wiki/Airy_function) +External links: +[DLMF 9.2](https://dlmf.nist.gov/9.2), +[Wikipedia](https://en.wikipedia.org/wiki/Airy_function) See also: [`airyaix`](@ref), [`airyaiprime`](@ref), [`airybi`](@ref) """ function airyai end airyai(z::Complex{Float64}) = _airy(z, Int32(0), Int32(1)) -""" +@doc raw""" airyaiprime(x) -Derivative of the Airy function of the first kind ``\\operatorname{Ai}'(x)``. +Derivative of the Airy function of the first kind ``\operatorname{Ai}'(x)``. -External links: [DLMF](https://dlmf.nist.gov/9.2), [Wikipedia](https://en.wikipedia.org/wiki/Airy_function) +External links: +[DLMF 9.2](https://dlmf.nist.gov/9.2), +[Wikipedia](https://en.wikipedia.org/wiki/Airy_function) See also: [`airyaiprimex`](@ref), [`airyai`](@ref), [`airybi`](@ref) """ function airyaiprime end airyaiprime(z::Complex{Float64}) = _airy(z, Int32(1), Int32(1)) -""" +@doc raw""" airybi(x) -Airy function of the second kind ``\\operatorname{Bi}(x)``. +Airy function of the second kind ``\operatorname{Bi}(x)``. -External links: [DLMF](https://dlmf.nist.gov/9.2), [Wikipedia](https://en.wikipedia.org/wiki/Airy_function) +External links: +[DLMF 8.2](https://dlmf.nist.gov/9.2), +[Wikipedia](https://en.wikipedia.org/wiki/Airy_function) See also: [`airybix`](@ref), [`airybiprime`](@ref), [`airyai`](@ref) """ function airybi end airybi(z::Complex{Float64}) = _biry(z, Int32(0), Int32(1)) -""" +@doc raw""" airybiprime(x) -Derivative of the Airy function of the second kind ``\\operatorname{Bi}'(x)``. +Derivative of the Airy function of the second kind ``\operatorname{Bi}'(x)``. -External links: [DLMF](https://dlmf.nist.gov/9.2), [Wikipedia](https://en.wikipedia.org/wiki/Airy_function) +External links: +[DLMF 9.2](https://dlmf.nist.gov/9.2), +[Wikipedia](https://en.wikipedia.org/wiki/Airy_function) See also: [`airybiprimex`](@ref), [`airybi`](@ref), [`airyai`](@ref) """ function airybiprime end airybiprime(z::Complex{Float64}) = _biry(z, Int32(1), Int32(1)) -""" +@doc raw""" 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. +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](https://dlmf.nist.gov/9.2), [Wikipedia](https://en.wikipedia.org/wiki/Airy_function) +External links: +[DLMF 9.2](https://dlmf.nist.gov/9.2), +[Wikipedia](https://en.wikipedia.org/wiki/Airy_function) See also: [`airyai`](@ref), [`airyaiprime`](@ref), [`airybi`](@ref) """ function airyaix end airyaix(z::Complex{Float64}) = _airy(z, Int32(0), Int32(2)) -""" +@doc raw""" 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. +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](https://dlmf.nist.gov/9.2), [Wikipedia](https://en.wikipedia.org/wiki/Airy_function) +External links: +[DLMF 9.2](https://dlmf.nist.gov/9.2), +[Wikipedia](https://en.wikipedia.org/wiki/Airy_function) See also: [`airyaiprime`](@ref), [`airyai`](@ref), [`airybi`](@ref) """ function airyaiprimex end airyaiprimex(z::Complex{Float64}) = _airy(z, Int32(1), Int32(2)) -""" +@doc raw""" 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|}``. +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](https://dlmf.nist.gov/9.2), [Wikipedia](https://en.wikipedia.org/wiki/Airy_function) +External links: +[DLMF 9.2](https://dlmf.nist.gov/9.2), +[Wikipedia](https://en.wikipedia.org/wiki/Airy_function) See also: [`airybi`](@ref), [`airybiprime`](@ref), [`airyai`](@ref) """ function airybix end airybix(z::Complex{Float64}) = _biry(z, Int32(0), Int32(2)) -""" +@doc raw""" 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|}``. +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](https://dlmf.nist.gov/9.2), [Wikipedia](https://en.wikipedia.org/wiki/Airy_function) +External links: +[DLMF 9.2](https://dlmf.nist.gov/9.2), +[Wikipedia](https://en.wikipedia.org/wiki/Airy_function) See also: [`airybiprime`](@ref), [`airybi`](@ref), [`airyai`](@ref) """ @@ -300,7 +320,9 @@ Bessel function of the third kind of order `nu` (the Hankel function). `k` is ei selecting [`hankelh1`](@ref) or [`hankelh2`](@ref), respectively. `k` defaults to 1 if it is omitted. -External links: [DLMF](https://dlmf.nist.gov/10.2.5) and [DLMF](https://dlmf.nist.gov/10.2.6), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1) +External links: +[DLMF 10.2.5](https://dlmf.nist.gov/10.2.5) and [DLMF 10.2.6](https://dlmf.nist.gov/10.2.6), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1) See also: [`besselhx`](@ref) for an exponentially scaled variant. """ @@ -326,20 +348,22 @@ function besselh(nu::Float64, k::Integer, x::Float64) end end -""" +@doc raw""" besselhx(nu, [k=1,] z) -Compute the scaled Hankel function ``\\exp(∓iz) H_ν^{(k)}(z)``, where +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 +proportional to ``\exp(∓iz)/\sqrt{z}`` for large ``|z|``, and so the [`besselh`](@ref) 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](https://dlmf.nist.gov/10.2), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1) +External links: +[DLMF 10.2](https://dlmf.nist.gov/10.2), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1) See also: [`besselh`](@ref) """ @@ -438,14 +462,19 @@ function besselyx(nu::Float64, z::Complex{Float64}) end end -""" +@doc raw""" besseli(nu, x) -Modified Bessel function of the first kind of order `nu`, ``I_\\nu(x)``. +Modified Bessel function of the first kind of order `nu`, ``I_\nu(x)``. -External links: [DLMF](https://dlmf.nist.gov/10.25.2), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I%CE%B1,_K%CE%B1) +External links: +[DLMF 10.25.2](https://dlmf.nist.gov/10.25.2), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I%CE%B1,_K%CE%B1) -See also: [`besselix(nu,x)`](@ref SpecialFunctions.besselix), [`besselj(nu,x)`](@ref SpecialFunctions.besselj), [`besselk(nu,x)`](@ref SpecialFunctions.besselk) +See also: +[`besselix(nu,x)`](@ref SpecialFunctions.besselix), +[`besselj(nu,x)`](@ref SpecialFunctions.besselj), +[`besselk(nu,x)`](@ref SpecialFunctions.besselk) """ function besseli(nu::Real, x::AbstractFloat) if x < 0 && !isinteger(nu) @@ -454,14 +483,19 @@ function besseli(nu::Real, x::AbstractFloat) real(besseli(float(nu), complex(x))) end -""" +@doc raw""" besselix(nu, x) -Scaled modified Bessel function of the first kind of order `nu`, ``I_\\nu(x) e^{- | \\operatorname{Re}(x) |}``. +Scaled modified Bessel function of the first kind of order `nu`, ``I_\nu(x) e^{- | \operatorname{Re}(x) |}``. -External links: [DLMF](https://dlmf.nist.gov/10.25.2), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I%CE%B1,_K%CE%B1) +External links: +[DLMF 10.25.2](https://dlmf.nist.gov/10.25.2), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I%CE%B1,_K%CE%B1) -See also: [`besseli(nu,x)`](@ref SpecialFunctions.besseli), [`besselj(nu,x)`](@ref SpecialFunctions.besselj), [`besselk(nu,x)`](@ref SpecialFunctions.besselk) +See also: +[`besseli(nu,x)`](@ref SpecialFunctions.besseli), +[`besselj(nu,x)`](@ref SpecialFunctions.besselj), +[`besselk(nu,x)`](@ref SpecialFunctions.besselk) """ function besselix(nu::Real, x::AbstractFloat) if x < 0 && !isinteger(nu) @@ -470,14 +504,19 @@ function besselix(nu::Real, x::AbstractFloat) real(besselix(float(nu), complex(x))) end -""" +@doc raw""" besselj(nu, x) -Bessel function of the first kind of order `nu`, ``J_\\nu(x)``. +Bessel function of the first kind of order `nu`, ``J_\nu(x)``. -External links: [DLMF](https://dlmf.nist.gov/10.2.2), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind:_J%CE%B1) +External links: +[DLMF 10.2.2](https://dlmf.nist.gov/10.2.2), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind:_J%CE%B1) -See also: [`besseljx(nu,x)`](@ref SpecialFunctions.besseljx), [`besseli(nu,x)`](@ref SpecialFunctions.besseli), [`besselk(nu,x)`](@ref SpecialFunctions.besselk) +See also: +[`besseljx(nu,x)`](@ref SpecialFunctions.besseljx), +[`besseli(nu,x)`](@ref SpecialFunctions.besseli), +[`besselk(nu,x)`](@ref SpecialFunctions.besselk) """ function besselj(nu::Real, x::AbstractFloat) if isinteger(nu) @@ -490,14 +529,19 @@ function besselj(nu::Real, x::AbstractFloat) real(besselj(float(nu), complex(x))) end -""" +@doc raw""" besseljx(nu, x) -Scaled Bessel function of the first kind of order `nu`, ``J_\\nu(x) e^{- | \\operatorname{Im}(x) |}``. +Scaled Bessel function of the first kind of order `nu`, ``J_\nu(x) e^{- | \operatorname{Im}(x) |}``. -External links: [DLMF](https://dlmf.nist.gov/10.2.2), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind:_J%CE%B1) +External links: +[DLMF 10.2.2](https://dlmf.nist.gov/10.2.2), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind:_J%CE%B1) -See also: [`besselj(nu,x)`](@ref SpecialFunctions.besselj), [`besseli(nu,x)`](@ref SpecialFunctions.besseli), [`besselk(nu,x)`](@ref SpecialFunctions.besselk) +See also: +[`besselj(nu,x)`](@ref SpecialFunctions.besselj), +[`besseli(nu,x)`](@ref SpecialFunctions.besseli), +[`besselk(nu,x)`](@ref SpecialFunctions.besselk) """ function besseljx(nu::Real, x::AbstractFloat) if x < 0 && !isinteger(nu) @@ -506,14 +550,19 @@ function besseljx(nu::Real, x::AbstractFloat) real(besseljx(float(nu), complex(x))) end -""" +@doc raw""" besselk(nu, x) -Modified Bessel function of the second kind of order `nu`, ``K_\\nu(x)``. +Modified Bessel function of the second kind of order `nu`, ``K_\nu(x)``. -External links: [DLMF](https://dlmf.nist.gov/10.25.3), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I%CE%B1,_K%CE%B1) +External links: +[DLMF 10.25.3](https://dlmf.nist.gov/10.25.3), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I%CE%B1,_K%CE%B1) -See also: See also: [`besselkx(nu,x)`](@ref SpecialFunctions.besselkx), [`besseli(nu,x)`](@ref SpecialFunctions.besseli), [`besselj(nu,x)`](@ref SpecialFunctions.besselj) +See also: +[`besselkx(nu,x)`](@ref SpecialFunctions.besselkx), +[`besseli(nu,x)`](@ref SpecialFunctions.besseli), +[`besselj(nu,x)`](@ref SpecialFunctions.besselj) """ function besselk(nu::Real, x::AbstractFloat) if x < 0 @@ -524,14 +573,19 @@ function besselk(nu::Real, x::AbstractFloat) real(besselk(float(nu), complex(x))) end -""" +@doc raw""" besselkx(nu, x) -Scaled modified Bessel function of the second kind of order `nu`, ``K_\\nu(x) e^x``. +Scaled modified Bessel function of the second kind of order `nu`, ``K_\nu(x) e^x``. -External links: [DLMF](https://dlmf.nist.gov/10.25.3), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I%CE%B1,_K%CE%B1) +External links: +[DLMF 10.25.3](https://dlmf.nist.gov/10.25.3), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I%CE%B1,_K%CE%B1) -See also: [`besselk(nu,x)`](@ref SpecialFunctions.besselk), [`besseli(nu,x)`](@ref SpecialFunctions.besseli), [`besselj(nu,x)`](@ref SpecialFunctions.besselj) +See also: +[`besselk(nu,x)`](@ref SpecialFunctions.besselk), +[`besseli(nu,x)`](@ref SpecialFunctions.besseli), +[`besselj(nu,x)`](@ref SpecialFunctions.besselj) """ function besselkx(nu::Real, x::AbstractFloat) if x < 0 @@ -547,7 +601,9 @@ end Bessel function of the second kind of order `nu`, ``Y_\\nu(x)``. -External links: [DLMF](https://dlmf.nist.gov/10.2.3), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_second_kind:_Y%CE%B1) +External links: +[DLMF 10.2.3](https://dlmf.nist.gov/10.2.3), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_second_kind:_Y%CE%B1) See also [`besselyx(nu,x)`](@ref SpecialFunctions.besselyx) for a scaled variant. """ @@ -566,7 +622,9 @@ end Scaled Bessel function of the second kind of order `nu`, ``Y_\\nu(x) e^{- | \\operatorname{Im}(x) |}``. -External links: [DLMF](https://dlmf.nist.gov/10.2.3), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_second_kind:_Y%CE%B1) +External links: +[DLMF 10.2.3](https://dlmf.nist.gov/10.2.3), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_second_kind:_Y%CE%B1) See also [`bessely(nu,x)`](@ref SpecialFunctions.bessely) """ @@ -615,7 +673,9 @@ besselh(nu::Float32, k::Integer, x::Float32) = Complex{Float32}(besselh(Float64( Bessel function of the first kind of order 0, ``J_0(x)``. -External links: [DLMF](https://dlmf.nist.gov/10.2.2), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind:_J%CE%B1) +External links: +[DLMF 10.2.2](https://dlmf.nist.gov/10.2.2), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind:_J%CE%B1) See also: [`besselj(nu,x)`](@ref SpecialFunctions.besselj) """ @@ -630,7 +690,9 @@ end Bessel function of the first kind of order 1, ``J_1(x)``. -External links: [DLMF](https://dlmf.nist.gov/10.2.2), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind:_J%CE%B1) +External links: +[DLMF 10.2.2](https://dlmf.nist.gov/10.2.2), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind:_J%CE%B1) See also: [`besselj(nu,x)`](@ref SpecialFunctions.besselj) """ @@ -651,7 +713,9 @@ end Bessel function of the second kind of order 0, ``Y_0(x)``. -External links: [DLMF](https://dlmf.nist.gov/10.2.3), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_second_kind:_Y%CE%B1) +External links: +[DLMF 10.2.3](https://dlmf.nist.gov/10.2.3), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_second_kind:_Y%CE%B1) See also: [`bessely(nu,x)`](@ref SpecialFunctions.bessely) """ @@ -669,7 +733,9 @@ end Bessel function of the second kind of order 1, ``Y_1(x)``. -External links: [DLMF](https://dlmf.nist.gov/10.2.3), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_second_kind:_Y%CE%B1) +External links: +[DLMF 10.2.3](https://dlmf.nist.gov/10.2.3), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_second_kind:_Y%CE%B1) See also: [`bessely(nu,x)`](@ref SpecialFunctions.bessely) """ @@ -709,9 +775,9 @@ end """ 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. +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. """ sphericalbessely(nu, x::T) where {T} = √((float(T))(π)/2x) * bessely(nu + one(nu)/2, x) @@ -720,7 +786,9 @@ sphericalbessely(nu, x::T) where {T} = √((float(T))(π)/2x) * bessely(nu + one Bessel function of the third kind of order `nu`, ``H^{(1)}_\\nu(x)``. -External links: [DLMF](https://dlmf.nist.gov/10.2.5), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1) +External links: +[DLMF 10.2.5](https://dlmf.nist.gov/10.2.5), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1) See also: [`hankelh1x`](@ref SpecialFunctions.hankelh1x) """ @@ -731,7 +799,9 @@ hankelh1(nu, z) = besselh(nu, 1, z) Bessel function of the third kind of order `nu`, ``H^{(2)}_\\nu(x)``. -External links: [DLMF](https://dlmf.nist.gov/10.2.6), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1) +External links: +[DLMF 10.2.6](https://dlmf.nist.gov/10.2.6), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1) See also: [`hankelh2x(nu,x)`](@ref SpecialFunctions.hankelh2x) """ @@ -742,28 +812,35 @@ hankelh2(nu, z) = besselh(nu, 2, z) Scaled Bessel function of the third kind of order `nu`, ``H^{(1)}_\\nu(x) e^{-x i}``. -External links: [DLMF](https://dlmf.nist.gov/10.2.5), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1) +External links: +[DLMF 10.2.5](https://dlmf.nist.gov/10.2.5), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1) See also: [`hankelh1`](@ref SpecialFunctions.hankelh1) """ hankelh1x(nu, z) = besselhx(nu, 1, z) -""" +@doc raw""" hankelh2x(nu, x) -Scaled Bessel function of the third kind of order `nu`, ``H^{(2)}_\\nu(x) e^{x i}``. +Scaled Bessel function of the third kind of order `nu`, ``H^{(2)}_\nu(x) e^{x i}``. -External links: [DLMF](https://dlmf.nist.gov/10.2.6), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1) +External links: +[DLMF 10.2.6](https://dlmf.nist.gov/10.2.6), +[Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1) See also: [`hankelh2(nu,x)`](@ref SpecialFunctions.hankelh2) """ hankelh2x(nu, z) = besselhx(nu, 2, z) -""" +@doc raw""" 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}``. +Following convention: +```math +\operatorname{jinc}{x} = \frac{2 J_1({\pi x})}{\pi x}. +``` Sometimes known as sombrero or besinc function. External links: [Wikipedia](https://en.wikipedia.org/wiki/Sombrero_function) diff --git a/src/beta_inc.jl b/src/beta_inc.jl index e972f6a5..c6f09e03 100644 --- a/src/beta_inc.jl +++ b/src/beta_inc.jl @@ -4,10 +4,10 @@ const exparg_n = log(nextfloat(floatmin(Float64))) const exparg_p = log(prevfloat(floatmax(Float64))) #COMPUTE log(gamma(b)/gamma(a+b)) when b >= 8 -""" +@doc raw""" loggammadiv(a,b) -Computes ``log(\\Gamma(b)/\\Gamma(a+b))`` when b >= 8 +Computes ``\log(\Gamma(b)/\Gamma(a+b))`` when `b >= 8` """ loggammadiv(a::Number, b::Number) = _loggammadiv(promote(float(a), float(b))...) @@ -48,8 +48,7 @@ end """ stirling_corr(a0,b0) -Compute stirling(a0) + stirling(b0) - stirling(a0 + b0) -for a0, b0 >= 8 +Compute `stirling(a0) + stirling(b0) - stirling(a0 + b0)` for `a0, b0 >= 8` """ function stirling_corr(a0::Float64, b0::Float64) a = min(a0, b0) @@ -74,10 +73,10 @@ function stirling_corr(a0::Float64, b0::Float64) return @horner(t, .833333333333333E-01, -.277777777760991E-02, .793650666825390E-03, -.595202931351870E-03, .837308034031215E-03, -.165322962780713E-02)/a + w end -""" +@doc raw""" esum(mu,x) -Compute ``e^{mu+x}`` +Compute ``e^{\mu+x}`` """ function esum(mu::Float64, x::Float64) if x > 0.0 @@ -93,10 +92,10 @@ function esum(mu::Float64, x::Float64) end end -""" +@doc raw""" beta_integrand(a, b, x, y, mu=0.0) -Compute ``e^{mu} * x^{a}y^{b}/B(a,b)`` +Compute ``e^{\mu} x^a y^b / B(a,b)`` """ function beta_integrand(a::Float64, b::Float64, x::Float64, y::Float64, mu::Float64=0.0) a0, b0 = minmax(a,b) @@ -178,13 +177,15 @@ function beta_integrand(a::Float64, b::Float64, x::Float64, y::Float64, mu::Floa end end -""" +@doc raw""" 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`` +It is assumed that ``\lambda = (a+b)*y - b`` -External links: [DLMF](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) +External links: +[DLMF 8.17.22](https://dlmf.nist.gov/8.17.22), +[Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) See also: [`beta_inc`](@ref) @@ -214,45 +215,47 @@ function beta_inc_cont_fraction(a::Float64, b::Float64, x::Float64, y::Float64, #CONT FRACTION while true - n += 1.0 - t = n/a - w = n*(b - n)*x - e = a/s - alpha = (p*(p+c0)*e*e)*(w*x) - e = (1.0 + t)/(c1 + 2*t) - beta = n + w/s +e*(c + n*yp1) - p = 1.0 + t - s += 2.0 - - #update an, bn, anp1, bnp1 - t = alpha*an + beta*anp1 - an = anp1 - anp1 = t - t = alpha*bn + beta*bnp1 - bn = bnp1 - bnp1 = t - - r0 = r - r = anp1/bnp1 - if abs(r - r0) <= epps*r - break - end - #rescale - an /= bnp1 - bn /= bnp1 - anp1 = r - bnp1 = 1.0 + n += 1.0 + t = n/a + w = n*(b - n)*x + e = a/s + alpha = (p*(p+c0)*e*e)*(w*x) + e = (1.0 + t)/(c1 + 2*t) + beta = n + w/s +e*(c + n*yp1) + p = 1.0 + t + s += 2.0 + + #update an, bn, anp1, bnp1 + t = alpha*an + beta*anp1 + an = anp1 + anp1 = t + t = alpha*bn + beta*bnp1 + bn = bnp1 + bnp1 = t + + r0 = r + r = anp1/bnp1 + if abs(r - r0) <= epps*r + break + end + #rescale + an /= bnp1 + bn /= bnp1 + anp1 = r + bnp1 = 1.0 end return ans*r end -""" +@doc raw""" 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``. +It is assumed that ``\lambda = (a+b)*y - b``. -External links: [DLMF](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) +External links: +[DLMF 8.17.22](https://dlmf.nist.gov/8.17.22), +[Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) See also: [`beta_inc`](@ref) @@ -353,7 +356,9 @@ end 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](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) +External links: +[DLMF 8.17.22](https://dlmf.nist.gov/8.17.22), +[Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) See also: [`beta_inc`](@ref) """ @@ -431,7 +436,9 @@ end Variant of `BPSER(A,B,X,EPS)`. -External links: [DLMF](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) +External links: +[DLMF 8.17.22](https://dlmf.nist.gov/8.17.22), +[Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) See also: [`beta_inc`](@ref) @@ -475,7 +482,9 @@ end Another variant of `BPSER(A,B,X,EPS)`. -External links: [DLMF](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) +External links: +[DLMF 8.17.22](https://dlmf.nist.gov/8.17.22), +[Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) See also: [`beta_inc`](@ref) @@ -511,14 +520,16 @@ function beta_inc_power_series1(a::Float64, b::Float64, x::Float64, epps::Float6 end #B .LE. 1 OR B*X .LE. 0.7 -""" +@doc raw""" beta_inc_power_series(a, b, x, epps) -Computes ``I_x(a,b)`` using power series : +Computes ``I_x(a,b)`` using power series: ```math -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}) +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](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) +External links: +[DLMF 8.17.22](https://dlmf.nist.gov/8.17.22), +[Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) See also: [`beta_inc`](@ref) @@ -616,9 +627,11 @@ end 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](https://dlmf.nist.gov/8.17.20). +A generalised version of [DLMF 8.17.20](https://dlmf.nist.gov/8.17.20). -External links: [DLMF](https://dlmf.nist.gov/8.17.20), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) +External links: +[DLMF 8.17.20](https://dlmf.nist.gov/8.17.20), +[Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) See also: [`beta_inc`](@ref) """ @@ -714,17 +727,19 @@ end #DLMF : https://dlmf.nist.gov/8.17#E1 #Wikipedia : https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function -""" +@doc raw""" 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 ```math -I_{x}(a,b) = \\frac{1}{B(a,b)} \\int_{0}^{x} t^{a-1}(1-t)^{b-1} dt, +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)``. +where ``B(a,b) = \Gamma(a)\Gamma(b)/\Gamma(a+b)``. -External links: [DLMF](https://dlmf.nist.gov/8.17.1), [Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) +External links: +[DLMF 8.17.1](https://dlmf.nist.gov/8.17.1), +[Wikipedia](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) See also: [`beta_inc_inv`](@ref) """ @@ -1008,9 +1023,9 @@ function _beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64=1-p) prev = max(sq, fpu) end - adj = p_approx - tx = x - adj - while prev <= (sq = adj^2) || tx < 0.0 || tx > 1.0 + adj = p_approx + tx = x - adj + while prev <= (sq = adj^2) || tx < 0.0 || tx > 1.0 adj /= 3.0 tx = x - adj end @@ -1038,4 +1053,3 @@ function _beta_inc_inv(a::T, b::T, p::T, q::T) where {T<:Union{Float16, Float32} x, y = _beta_inc_inv(Float64(a), Float64(b), Float64(p), Float64(q)) T(x), T(y) end - diff --git a/src/betanc.jl b/src/betanc.jl index 7c13480e..7b44f28d 100644 --- a/src/betanc.jl +++ b/src/betanc.jl @@ -4,15 +4,15 @@ const errmax = 1e-15 #Russell Lenth, Algorithm AS 226: Computing Noncentral Beta Probabilities, #Applied Statistics,Volume 36, Number 2, 1987, pages 241-244 -""" - ncbeta_tail(x,a,b,lambda) +@doc raw""" + ncbeta_tail(x,a,b,lambda) Compute tail of the noncentral beta distribution. Uses the recursive relation ```math -I_{x}(a,b+1;0) = I_{x}(a,b;0) - \\Gamma(a+b)/\\Gamma(a+1)\\Gamma(b)x^{a}(1-x)^{b} +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. +and ``\Gamma(a+1) = a\Gamma(a)`` given in [DLMF 8.17.21](https://dlmf.nist.gov/8.17.21). """ function ncbeta_tail(a::Float64, b::Float64, lambda::Float64, x::Float64) if x <= 0.0 @@ -51,18 +51,20 @@ function ncbeta_tail(a::Float64, b::Float64, lambda::Float64, x::Float64) return ans end -""" +@doc raw""" 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. +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. ```math -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} +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`. ```math -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} +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} ``` """ function ncbeta_poisson(a::Float64, b::Float64, lambda::Float64, x::Float64) @@ -140,15 +142,16 @@ end #R Chattamvelli, R Shanmugam, Algorithm AS 310: Computing the Non-central Beta Distribution Function, #Applied Statistics, Volume 46, Number 1, 1997, pages 146-156 -""" - ncbeta(a,b,lambda,x) +@doc raw""" + ncbeta(a,b,lambda,x) Compute the CDF of the noncentral beta distribution given by ```math -I_{x}(a,b;\\lambda ) = \\sum_{j=0}^{\\infty}q(\\lambda/2,j)I_{x}(a+j,b;0) +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. +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. """ function ncbeta(a::Float64, b::Float64, lambda::Float64, x::Float64) ans = x @@ -165,27 +168,28 @@ function ncbeta(a::Float64, b::Float64, lambda::Float64, x::Float64) end end -""" +@doc raw""" ncF(x,v1,v2,lambda) Compute CDF of noncentral F distribution given by: ```math -F(x, v1, v2; lambda) = I_{v1*x/(v1*x + v2)}(v1/2, v2/2; \\lambda) +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. +where ``I_{x}(a,b; \lambda)`` is the noncentral beta function computed above. -Wikipedia: https://en.wikipedia.org/wiki/Noncentral_F-distribution +External links: +[Wikipedia](https://en.wikipedia.org/wiki/Noncentral_F-distribution) """ function ncF(x::Float64, v1::Float64, v2::Float64, lambda::Float64) return ncbeta(v1/2, v2/2, lambda, (v1*x)/(v1*x + v2)) end function ncbeta(a::T,b::T,lambda::T,x::T) where {T<:Union{Float16,Float32}} - T.(ncbeta(Float64(a),Float64(b),Float64(lambda),Float64(x))) + T.(ncbeta(Float64(a),Float64(b),Float64(lambda),Float64(x))) end function ncF(x::T,v1::T,v2::T,lambda::T) where {T<:Union{Float16,Float32}} - T.(ncF(Float64(x),Float64(v1),Float64(v2),Float64(lambda))) + T.(ncF(Float64(x),Float64(v1),Float64(v2),Float64(lambda))) end ncbeta(a::Real,b::Real,lambda::Real,x::Real) = ncbeta(promote(float(a),float(b),float(lambda),float(x))...) diff --git a/src/ellip.jl b/src/ellip.jl index 5065f60d..7636fcd7 100644 --- a/src/ellip.jl +++ b/src/ellip.jl @@ -15,13 +15,16 @@ Computes Complete Elliptic Integral of 1st kind ``K(m)`` for parameter ``m`` giv \quad \text{for} \quad m \in \left( -\infty, 1 \right] \, . ``` -External links: [DLMF](https://dlmf.nist.gov/19.2.8), [Wikipedia](https://en.wikipedia.org/wiki/Elliptic_integral#Notational_variants). +External links: +[DLMF 19.2.8](https://dlmf.nist.gov/19.2.8), +[Wikipedia](https://en.wikipedia.org/wiki/Elliptic_integral#Notational_variants). See also: [`ellipe(m)`](@ref SpecialFunctions.ellipe). # 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``. +- `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 @@ -189,16 +192,19 @@ Computes Complete Elliptic Integral of 2nd kind ``E(m)`` for parameter ``m`` giv \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] \, . +\quad \text{for} \quad m \in \left( -\infty, 1 \right] . ``` -External links: [DLMF](https://dlmf.nist.gov/19.2.8), [Wikipedia](https://en.wikipedia.org/wiki/Elliptic_integral#Complete_elliptic_integral_of_the_second_kind). +External links: +[DLMF 19.2.8](https://dlmf.nist.gov/19.2.8), +[Wikipedia](https://en.wikipedia.org/wiki/Elliptic_integral#Complete_elliptic_integral_of_the_second_kind). See also: [`ellipk(m)`](@ref SpecialFunctions.ellipk). # 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``. +- `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 diff --git a/src/erf.jl b/src/erf.jl index dd04e904..c878cf5d 100644 --- a/src/erf.jl +++ b/src/erf.jl @@ -86,13 +86,12 @@ Compute the error function of ``x``, defined by \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](https://dlmf.nist.gov/7.2.1), [Wikipedia](https://en.wikipedia.org/wiki/Error_function). +External links: +[DLMF 7.2.1](https://dlmf.nist.gov/7.2.1), +[Wikipedia](https://en.wikipedia.org/wiki/Error_function). -See also: [`erfc(x)`](@ref erfc), [`erfcx(x)`](@ref erfcx), +See also: +[`erfc(x)`](@ref erfc), [`erfcx(x)`](@ref erfcx), [`erfi(x)`](@ref erfi), [`dawson(x)`](@ref dawson), [`erfinv(x)`](@ref erfinv), [`erfcinv(x)`](@ref erfcinv). @@ -101,8 +100,12 @@ See also: [`erfc(x)`](@ref erfc), [`erfcx(x)`](@ref erfcx), [libm](https://en.wikipedia.org/wiki/C_mathematical_functions#libm). - `BigFloat`: C library for multiple-precision floating-point [MPFR](https://www.mpfr.org/) """ -erf +function erf end +""" + erf(x, y) +Accurate version of `erf(y) - erf(x)` (for real arguments only). +""" function erf(x::Real, y::Real) if abs(x) ≤ 1/√2 && abs(y) ≤ 1/√2 erf(y) - erf(x) @@ -129,7 +132,7 @@ Compute the complementary error function of ``x``, defined by This is the accurate version of `1-erf(x)` for large ``x``. -External links: [DLMF](https://dlmf.nist.gov/7.2.2), +External links: [DLMF 7.2.2](https://dlmf.nist.gov/7.2.2), [Wikipedia](https://en.wikipedia.org/wiki/Error_function#Complementary_error_function). See also: [`erf(x)`](@ref erf). @@ -155,7 +158,7 @@ Compute the scaled complementary error function of ``x``, defined by 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](https://dlmf.nist.gov/7.2.3), +External links: [DLMF 7.2.3](https://dlmf.nist.gov/7.2.3), [Wikipedia](https://en.wikipedia.org/wiki/Error_function#Complementary_error_function). See also: [`erfc(x)`](@ref erfc). @@ -204,7 +207,7 @@ Compute the Dawson function (scaled imaginary error function) of ``x``, defined This is the accurate version of ``\frac{\sqrt{\pi}}{2} e^{-x^2} \operatorname{erfi}(x)`` for large ``x``. -External links: [DLMF](https://dlmf.nist.gov/7.2.5), +External links: [DLMF 7.2.5](https://dlmf.nist.gov/7.2.5), [Wikipedia](https://en.wikipedia.org/wiki/Dawson_function). See also: [`erfi(x)`](@ref erfi). @@ -215,13 +218,13 @@ See also: [`erfi(x)`](@ref erfi). """ dawson -""" +@doc raw""" faddeeva(z) Compute the Faddeeva function of complex `z`, defined by -``e^{-z^2} \\operatorname{erfc}(-iz)``. +``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)``. +is equivalent to ``\operatorname{erfcx}(-iz)``. """ faddeeva @@ -563,11 +566,11 @@ end Compute the natural logarithm of the complementary error function of ``x``, that is ```math -\operatorname{logerfc}(x) = \operatorname{ln}(\operatorname{erfc}(x)) +\operatorname{logerfc}(x) = \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``. +This is the accurate version of ``\ln(\operatorname{erfc}(x))`` for large ``x``. External links: [Wikipedia](https://en.wikipedia.org/wiki/Error_function). @@ -590,11 +593,11 @@ end Compute the natural logarithm of the scaled complementary error function of ``x``, that is ```math -\operatorname{logerfcx}(x) = \operatorname{ln}(\operatorname{erfcx}(x)) +\operatorname{logerfcx}(x) = \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``. +This is the accurate version of ``\ln(\operatorname{erfcx}(x))`` for large and negative ``x``. External links: [Wikipedia](https://en.wikipedia.org/wiki/Error_function). @@ -615,7 +618,7 @@ end 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`. +`log(erf(x, y))`, which works for large `x`, `y`. External links: [Wikipedia](https://en.wikipedia.org/wiki/Error_function). diff --git a/src/expint.jl b/src/expint.jl index a0acd81d..63a86b19 100644 --- a/src/expint.jl +++ b/src/expint.jl @@ -58,7 +58,7 @@ function E₁_taylor_coefficients(::Type{T}, n::Integer) where {T<:Number} # iteratively compute the terms in the series, starting with k=1 term::T = 1 terms = T[-eulergamma, term] - for k=2:n + for k = 2:n term = -term * (k-1) / (k * k) push!(terms, term) end @@ -504,24 +504,33 @@ function _expint(ν::Number, z::Number, niter::Int=1000, ::Val{expscaled}=Val{fa end end -""" +@doc raw""" 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. +Computes the exponential integral +```math +\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](https://dlmf.nist.gov/8.19), [Wikipedia](https://en.wikipedia.org/wiki/Exponential_integral) +External links: +[DLMF 8.19](https://dlmf.nist.gov/8.19), +[Wikipedia](https://en.wikipedia.org/wiki/Exponential_integral) """ expint(ν::Number, z::Number, niter::Int=1000) = _expint(ν, z, niter, Val{false}()) -""" +@doc raw""" 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. +Computes the scaled exponential integral +```math +\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)`](@ref SpecialFunctions.expint) """ @@ -530,11 +539,15 @@ expintx(ν::Number, z::Number, niter::Int=1000) = _expint(ν, z, niter, Val{true ############################################################################## # expinti function Ei -""" +@doc raw""" 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. +Computes the exponential integral function +```math +\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. """ expinti(x::Real) = x ≤ 0 ? -expint(-x) : -real(expint(complex(-x))) diff --git a/src/gamma.jl b/src/gamma.jl index 227d9c37..5004404d 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -220,20 +220,22 @@ end # is equivalent to Mathematica's Zeta[s,z], and is equivalent to the # Hurwitz zeta function for real(z) > 0. -""" +@doc raw""" zeta(s, z) Generalized zeta function defined by ```math -\\zeta(s, z)=\\sum_{k=0}^\\infty \\frac{1}{((k+z)^2)^{s/2}}, +\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``, +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}``. +``\sum_{k=0}^\infty (k+z)^{-s}``. -The Riemann zeta function is recovered as ``\\zeta(s)=\\zeta(s,1)``. +The Riemann zeta function is recovered as ``\zeta(s) = \zeta(s,1)``. -External links: [Riemann zeta function](https://en.wikipedia.org/wiki/Riemann_zeta_function), [Hurwitz zeta function](https://en.wikipedia.org/wiki/Hurwitz_zeta_function) +External links: +[Riemann zeta function](https://en.wikipedia.org/wiki/Riemann_zeta_function), +[Hurwitz zeta function](https://en.wikipedia.org/wiki/Hurwitz_zeta_function) """ zeta(s::Number, z::Number) = _zeta(map(float, promote(s, z))...) @@ -335,8 +337,8 @@ end """ polygamma(m, x) -Compute the polygamma function of order `m` of argument `x` (the `(m+1)`th derivative of the -logarithm of `gamma(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](https://en.wikipedia.org/wiki/Polygamma_function) @@ -405,13 +407,13 @@ function _invdigamma(y::Float64) return x_new end -""" +@doc raw""" zeta(s) Riemann zeta function ```math -\\zeta(s)=\\sum_{n=1}^\\infty \\frac{1}{n^s}\\quad\\text{for}\\quad s\\in\\mathbb{C}. +\zeta(s) = \sum_{n=1}^\infty \frac{1}{n^s}\quad\text{for}\quad s\in\mathbb{C}. ``` External links: [Wikipedia](https://en.wikipedia.org/wiki/Riemann_zeta_function) @@ -483,10 +485,13 @@ function _zeta(x::BigFloat) return z end -""" +@doc raw""" eta(s) -Dirichlet eta function ``\\eta(s) = \\sum^\\infty_{n=1}(-1)^{n-1}/n^{s}``. +Dirichlet eta function +```math +\eta(s) = \sum_{n=1}^\infty (-1)^{n-1} / n^{s}. +``` """ eta(s::Number) = _eta(float(s)) @@ -551,7 +556,7 @@ Compute the gamma function for complex ``z``, defined by n! & \text{for} \quad z = n+1 \;, n = 0,1,2,\dots \\ - \int_0^\infty t^{z-1} {\mathrm e}^{-t} \, {\mathrm d}t + \int_0^\infty t^{z-1} e^{-t} \, \mathrm{d}t & \text{for} \quad \Re(z) > 0 \end{cases} ``` @@ -577,10 +582,11 @@ true ``` External links: -[DLMF](https://dlmf.nist.gov/5.2.1), +[DLMF 5.2.1](https://dlmf.nist.gov/5.2.1), [Wikipedia](https://en.wikipedia.org/wiki/Gamma_function). -See also: [`loggamma(z)`](@ref SpecialFunctions.loggamma) for ``\log \Gamma(z)`` and +See also: +[`loggamma(z)`](@ref SpecialFunctions.loggamma) for ``\log \Gamma(z)`` and [`gamma(a,z)`](@ref SpecialFunctions.gamma(::Number,::Number)) for the upper incomplete gamma function ``\Gamma(a,z)``. @@ -600,9 +606,9 @@ function gamma(n::Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64}) @inbounds return Float64(Base._fact_table64[n-1]) end -_gamma(x::Float64) = nan_dom_err(ccall((:tgamma, libopenlibm), Float64, (Float64,), x), x) -_gamma(x::Float32) = nan_dom_err(ccall((:tgammaf, libopenlibm), Float32, (Float32,), x), x) -_gamma(x::Float16) = Float16(_gamma(Float32(x))) +_gamma(x::Float64) = nan_dom_err(ccall((:tgamma, libopenlibm), Float64, (Float64,), x), x) +_gamma(x::Float32) = nan_dom_err(ccall((:tgammaf, libopenlibm), Float32, (Float32,), x), x) +_gamma(x::Float16) = Float16(_gamma(Float32(x))) function _gamma(x::BigFloat) isnan(x) && return x @@ -618,7 +624,7 @@ _gamma(z::Complex) = exp(loggamma(z)) logabsgamma(x) Compute the logarithm of absolute value of [`gamma`](@ref) for -[`Real`](@ref) `x` and returns a tuple `(log(abs(gamma(x))), sign(gamma(x)))`. +`Real` `x` and returns a tuple `(log(abs(gamma(x))), sign(gamma(x)))`. See also [`loggamma`](@ref). """ @@ -648,11 +654,11 @@ logfactorial(x::Integer) = x < 0 ? throw(DomainError(x, "`x` must be non-negativ """ loggamma(x) -Computes the logarithm of [`gamma`](@ref) for given `x`. If `x` is a `Real`, then it -throws a `DomainError` if `gamma(x)` is negative. +Computes the logarithm of [`gamma`](@ref) 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`` +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`](@ref) for real `x`. @@ -760,10 +766,10 @@ function loggamma_asymptotic(z::Complex{Float64}) 6.4102564102564102564102561e-03,-2.9550653594771241830065352e-02) end -""" +@doc raw""" beta(x, y) -Euler integral of the first kind ``\\operatorname{B}(x,y) = \\Gamma(x)\\Gamma(y)/\\Gamma(x+y)``. +Euler integral of the first kind ``\operatorname{B}(x,y) = \Gamma(x)\Gamma(y)/\Gamma(x+y)``. """ function beta(a::Number, b::Number) lab, sign = logabsbeta(a, b) @@ -778,10 +784,10 @@ if Base.MPFR.version() >= v"4.0.0" end end -""" +@doc raw""" logbeta(x, y) -Natural logarithm of the [`beta`](@ref) function ``\\log(|\\operatorname{B}(x,y)|)``. +Natural logarithm of the [`beta`](@ref) function ``\log(|\operatorname{B}(x,y)|)``. See also [`logabsbeta`](@ref). """ @@ -796,7 +802,8 @@ end """ logabsbeta(x, y) -Compute the natural logarithm of the absolute value of the [`beta`](@ref) function, returning a tuple `(log(abs(beta(x,y))), sign(beta(x,y)))` +Compute the natural logarithm of the absolute value of the [`beta`](@ref) +function, returning a tuple `(log(abs(beta(x,y))), sign(beta(x,y)))` See also [`logbeta`](@ref). """ @@ -838,10 +845,13 @@ logabsbeta(a::Number, b::Number) = loggamma(a) + loggamma(b) - loggamma(a + b), """ logabsbinomial(n, k) -Accurate natural logarithm of the absolute value of the [`binomial`](@ref) +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`](https://docs.julialang.org/en/v1/base/math/#Base.binomial) """ function logabsbinomial(n::T, k::T) where {T<:Integer} S = float(T) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index b196f8da..b992b19e 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -30,10 +30,11 @@ const d6 = [-.592166437353694E-03] const d70 = .344367606892378E-03 const d80 = -.652623918595309E-03 -""" - rgamma1pm1(a) +@doc raw""" + rgamma1pm1(a) + +Computation of ``1/\Gamma(a+1) - 1`` for `-0.5<=a<=1.5`: ``1/\Gamma (a+1) - 1``. -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)`. """ function rgamma1pm1(a::Float64) @@ -56,10 +57,10 @@ function rgamma1pm1(a::Float64) end end -""" +@doc raw""" rgammax(a,x) -Evaluation of ``1/\\Gamma(a) e^{-x} x^{a}``. Based on `DRCOMP` from the `NSWC` library. +Evaluation of ``1/\Gamma(a) e^{-x} x^{a}``. Based on `DRCOMP` from the `NSWC` library. """ function rgammax(a::Float64, x::Float64) if x == 0.0 @@ -81,10 +82,10 @@ function rgammax(a::Float64, x::Float64) end end -""" +@doc raw""" auxgam(x) -Compute function `g` in ``1/\\Gamma(x+1) = 1+x*(x-1)*g(x)``, `-1 <= x <= 1`. +Compute function `g` in ``1/\Gamma(x+1) = 1 + x (x-1) g(x)``, `-1 <= x <= 1`. """ function auxgam(x::Float64) @assert -1 <= x <= 1 @@ -96,10 +97,10 @@ function auxgam(x::Float64) end end -""" +@doc raw""" loggamma1p(x) -Compute ``log(\\Gamma(1+x))`` for `-1 < x <= 1`. +Compute ``\log(\Gamma(1+x))`` for `-1 < x <= 1`. """ function loggamma1p(x::Float64) @assert -1 < x <= 1 @@ -130,10 +131,10 @@ function chepolsum(x::Float64, a::Array{Float64,1}) end end -""" +@doc raw""" stirling_error(x) -Compute ``\\ln{\\Gamma(x)} - (x-0.5)*\\ln{x} + x - \\ln{(2\\pi)}/2``. Adapted from `stirling` in `IncgamFI`. +Compute ``\ln{\Gamma(x)} - (x-0.5)*\ln{x} + x - \ln{(2\pi)}/2``. Adapted from `stirling` in `IncgamFI`. """ function stirling_error(x::Float64) if x < floatmin(Float64)*1000.0 @@ -161,8 +162,9 @@ end gammax(x) ```math -\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. +\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} ``` """ @@ -176,10 +178,10 @@ function gammax(x::Float64) end end -""" +@doc raw""" lambdaeta(eta) -Compute the value of ``\\lambda`` satisfying ``\\eta^{2}/2 = \\lambda-1-\\log{\\lambda}``. +Compute the value of ``\lambda`` satisfying ``\eta^{2}/2 = \lambda-1-\log{\lambda}``. """ function lambdaeta(eta::Float64) s = eta*eta*0.5 @@ -217,10 +219,10 @@ function lambdaeta(eta::Float64) return la end -""" -Computing the first coefficient for the expansion : +@doc raw""" +Computing the first coefficient for the expansion: ```math -\\epsilon (\\eta_{0},a) = \\epsilon_{1} (\\eta_{0},a)/a + \\epsilon_{2} (\\eta_{0},a)/a^{2} + \\epsilon_{3} (\\eta_{0},a)/a^{3} +\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 """ @@ -241,10 +243,11 @@ function coeff1(eta::Float64) end return coeff1 end -""" -Computing the second coefficient for the expansion : + +@doc raw""" +Computing the second coefficient for the expansion: ```math -\\epsilon (\\eta_{0},a) = \\epsilon_{1} (\\eta_{0},a)/a + \\epsilon_{2} (\\eta_{0},a)/a^{2} + \\epsilon_{3} (\\eta_{0},a)/a^{3} +\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 """ @@ -289,10 +292,11 @@ function coeff2(eta::Float64) end return coeff2 end -""" -Computing the third and last coefficient for the expansion : + +@doc raw""" +Computing the third and last coefficient for the expansion: ```math -\\epsilon (\\eta_{0},a) = \\epsilon_{1} (\\eta_{0},a)/a + \\epsilon_{2} (\\eta_{0},a)/a^{2} + \\epsilon_{3} (\\eta_{0},a)/a^{3} +\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 """ @@ -367,16 +371,16 @@ function coeff3(eta::Float64) end -""" +@doc raw""" gamma_inc_cf(a, x, ind) -Computes ``P(a,x)`` by continued fraction expansion given by : +Computes ``P(a,x)`` by continued fraction expansion given by: ```math -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}}}}}}} +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](https://dlmf.nist.gov/8.9.2) +External links: [DLMF 8.9.2](https://dlmf.nist.gov/8.9.2) See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc) """ @@ -406,16 +410,17 @@ function gamma_inc_cf(a::Float64, x::Float64, ind::Integer) q = rgammax(a, x)*a2n return (1.0 - q, q) end -""" + +@doc raw""" gamma_inc_taylor(a, x, ind) -Compute ``P(a,x)`` using Taylor Series for P/R given by : +Compute ``P(a,x)`` using Taylor Series for P/R given by: ```math -R(a,x)/a * (1 + \\sum_{n=1}^{\\infty} x^{n}/((a+1)(a+2)...(a+n))) +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](https://dlmf.nist.gov/8.11.2) +External links: [DLMF 8.11.2](https://dlmf.nist.gov/8.11.2) See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc) """ @@ -428,14 +433,14 @@ function gamma_inc_taylor(a::Float64, x::Float64, ind::Integer) wk[1] = t loop = 2 for indx = 2:20 - apn += 1.0 - t *= x/apn - if t <= 1.0e-3 - loop = indx - flag = true - break - end - wk[indx] = t + apn += 1.0 + t *= x/apn + if t <= 1.0e-3 + loop = indx + flag = true + break + end + wk[indx] = t end if !flag loop = 20 @@ -443,29 +448,30 @@ function gamma_inc_taylor(a::Float64, x::Float64, ind::Integer) sm = t tol = 0.5*acc #tolerance while true - apn += 1.0 - t *= x/apn - sm += t - if t <= tol - break - end + apn += 1.0 + t *= x/apn + sm += t + if t <= tol + break + end end for j = loop-1:-1:1 - sm += wk[j] + sm += wk[j] end p = (rgammax(a, x)/a)*(1.0 + sm) return (p, 1.0 - p) end -""" + +@doc raw""" gamma_inc_asym(a, x, ind) Compute ``P(a,x)`` using asymptotic expansion given by: ```math -R(a,x)/a * (1 + \\sum_{n=1}^{N-1}(a_{n}/x^{n} + \\Theta _{n}a_{n}/x^{n})) +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](https://dlmf.nist.gov/8.11.2) +External links: [DLMF 8.11.2](https://dlmf.nist.gov/8.11.2) See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc) """ @@ -505,16 +511,17 @@ function gamma_inc_asym(a::Float64, x::Float64, ind::Integer) q = (rgammax(a, x)/x)*(1.0 + sm) return (1.0 - q, q) end -""" + +@doc raw""" gamma_inc_taylor_x(a,x,ind) -Computes ``P(a,x)`` based on Taylor expansion of ``P(a,x)/x**a`` given by: +Computes ``P(a,x)`` based on Taylor expansion of ``P(a,x)/x^a`` given by: ```math -J = -a * \\sum_{1}^{\\infty} (-x)^{n}/((a+n)n!) +J = -a * \sum_{1}^{\infty} (-x)^{n}/((a+n)n!) ``` -and ``P(a,x)/x**a`` is given by : +and ``P(a,x)/x^a`` is given by: ```math -(1 - J)/ \\Gamma(a+1) +(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. @@ -552,21 +559,24 @@ function gamma_inc_taylor_x(a::Float64, x::Float64, ind::Integer) end end -""" +@doc raw""" gamma_inc_minimax(a,x,z) Compute ``P(a,x)`` using minimax approximations given by : ```math -1/2 * erfc(\\sqrt{y}) - e^{-y}/\\sqrt{2\\pi*a}* T(a,\\lambda) +1/2 * \operatorname{erfc}(\sqrt{y}) - e^{-y}/\sqrt{2\pi a} ⋅ T(a,\lambda) ``` where ```math -T(a,\\lambda) = \\sum_{0}^{N} c_{k}(z)a^{-k} +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 . +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](https://dlmf.nist.gov/8.12.8) +External links: [DLMF 8.12.8](https://dlmf.nist.gov/8.12.8) See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc) """ @@ -616,19 +626,20 @@ function gamma_inc_minimax(a::Float64, x::Float64, z::Float64) end end -""" +@doc raw""" gamma_inc_temme(a, x, z) Compute ``P(a,x)`` using Temme's expansion given by : ```math -1/2 * erfc(\\sqrt{y}) - e^{-y}/\\sqrt{2\\pi*a}* T(a,\\lambda) +1/2 * \operatorname{erfc}(\sqrt{y}) - e^{-y}/\sqrt{2\pi a} ⋅ T(a,\lambda) ``` where ```math -T(a,\\lambda) = \\sum_{0}^{N} c_{k}(z)a^{-k} +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. +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](https://dlmf.nist.gov/8.12.8) +External links: [DLMF 8.12.8](https://dlmf.nist.gov/8.12.8) See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc) """ @@ -650,20 +661,20 @@ function gamma_inc_temme(a::Float64, x::Float64, z::Float64) end end -""" +@doc raw""" gamma_inc_temme_1(a, x, z, ind) -Computes ``P(a,x)`` using simplified Temme expansion near ``y=0`` by : +Computes ``P(a,x)`` using simplified Temme expansion near ``y=0`` by: ```math -E(y) - (1 - y)/\\sqrt{2\\pi*a} * T(a,\\lambda) +E(y) - (1 - y)/\sqrt{2\pi a} ⋅ T(a,\lambda) ``` where ```math -E(y) = 1/2 - (1 - y/3)*(\\sqrt(y/\\pi)) +E(y) = 1/2 - (1 - y/3) ⋅ \sqrt{y/\pi} ``` -Used instead of it's previous function when ``\\sigma <= e_{0}/\\sqrt{a}``. +Used instead of it's previous function when ``\sigma \leq e_{0}/\sqrt{a}``. -External links: [DLMF](https://dlmf.nist.gov/8.12.8) +External links: [DLMF 8.12.8](https://dlmf.nist.gov/8.12.8) """ function gamma_inc_temme_1(a::Float64, x::Float64, z::Float64, ind::Integer) iop = ind + 1 @@ -741,16 +752,16 @@ function gamma_inc_fsum(a::Float64, x::Float64) end -""" +@doc raw""" 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: ```math -x = r\\left[1 + a\\sum_{k=1}^{\\infty}\\frac{(-x)^{n}}{(a+n)n!}\\right]^{-1/a}, +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}``. +where ``r = (p\Gamma(1+a))^{1/a}`` +Inverting this relation we obtain ``x = r + \sum_{k=2}^{\infty} c_{k} r^{k}``. """ function gamma_inc_inv_psmall(a::Float64, logr::Float64) r = exp(logr) @@ -769,15 +780,15 @@ function gamma_inc_inv_psmall(a::Float64, logr::Float64) return x0 end -""" +@doc raw""" 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)``. +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 ```math -x \\sim x_{0} - L + b \\sum_{k=1}^{\\infty} d_{k}/x_{0}^{k} +x \sim x_{0} - L + b \sum_{k=1}^{\infty} d_{k}/x_{0}^{k} ``` -where ``b = 1-a``, ``L = \\ln{x_0}``. +where ``b = 1-a``, ``L = \ln{x_0}``. """ function gamma_inc_inv_qsmall(a::Float64, q::Float64, qgammaxa::Float64) b = 1.0 - a @@ -803,21 +814,22 @@ function gamma_inc_inv_qsmall(a::Float64, q::Float64, qgammaxa::Float64) return x0 end -""" +@doc raw""" gamma_inc_inv_alarge(a, minpq, pcase) Compute `x0` - initial approximation when `a` is large. -The inversion problem is rewritten as : +The inversion problem is rewritten as: ```math -0.5 \\operatorname{erfc}(\\eta \\sqrt{a/2}) + R_{a}(\\eta) = q +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)`` +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: ```math -\\epsilon(\\eta_{0},a) = \\epsilon_{1}(\\eta_{0},a)/a + \\epsilon_{2}(\\eta_{0},a)/a^{2} + \\epsilon_{3}(\\eta_{0},a)/a^{3} + ... +\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. +which is calculated by [`coeff1`](@ref), [`coeff2`](@ref) and [`coeff3`](@ref) 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. """ function gamma_inc_inv_alarge(a::Float64, minpq::Float64, pcase::Bool) r = erfcinv(2*minpq) @@ -833,27 +845,32 @@ end # Volume 12 Issue 4, Dec. 1986 Pages 377-393 # doi>10.1145/22721.23109 -""" +@doc raw""" 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 = P(a,x)`` is the Incomplete gamma function ratio given by: ```math -P(a,x)=\\frac{1}{\\Gamma (a)} \\int_{0}^{x} e^{-t}t^{a-1} dt. +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: +and ``q = Q(a,x)`` is the Incomplete gamma function ratio given by: ```math -Q(a,x)=\\frac{1}{\\Gamma (a)} \\int_{x}^{\\infty} e^{-t}t^{a-1} dt. +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)``. +``\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. +`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](https://dlmf.nist.gov/8.2.4), [Wikipedia](https://en.wikipedia.org/wiki/Incomplete_gamma_function) +External links: +[DLMF 8.2.4](https://dlmf.nist.gov/8.2.4), +[Wikipedia](https://en.wikipedia.org/wiki/Incomplete_gamma_function) -See also [`gamma(z)`](@ref SpecialFunctions.gamma), [`gamma_inc_inv(a,p,q)`](@ref SpecialFunctions.gamma_inc_inv) +See also +[`gamma(z)`](@ref SpecialFunctions.gamma), +[`gamma_inc_inv(a,p,q)`](@ref SpecialFunctions.gamma_inc_inv) """ gamma_inc(a::Real,x::Real,ind::Integer=0) = _gamma_inc(promote(float(a),float(x))...,ind) @@ -949,7 +966,10 @@ end function _gamma_inc(a::BigFloat,x::BigFloat,ind::Integer) #BigFloat version from GNU MPFR wrapped via ccall z = BigFloat() - ccall((:mpfr_gamma_inc, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Int32), z, a, x, ROUNDING_MODE[]) + ccall((:mpfr_gamma_inc, :libmpfr), + Int32, + (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Int32), + z, a, x, ROUNDING_MODE[]) q = z/gamma(a) return (1.0 - q, q) end @@ -981,9 +1001,12 @@ end """ 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``. +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](https://dlmf.nist.gov/8.2.4), [Wikipedia](https://en.wikipedia.org/wiki/Incomplete_gamma_function) +External links: +[DLMF 8.2.4](https://dlmf.nist.gov/8.2.4), +[Wikipedia](https://en.wikipedia.org/wiki/Incomplete_gamma_function) See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc). """ @@ -1096,23 +1119,25 @@ promotereal(x::Real, y::Number) = let (u, v) = promote(x, y); real(u), v; end promotereal(x::Number, y::Real) = let (u, v) = promote(x, y); u, real(v); end promotereal(x, y) = promote(x,y) -""" +@doc raw""" gamma(a,x) Returns the upper incomplete gamma function ```math -\\Gamma(a,x) = \\int_x^\\infty t^{a-1} e^{-t} dt \\, +\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)`](@ref) corresponds to ``\\Gamma(a) = \\Gamma(a,0)``. +(The ordinary gamma function [`gamma(x)`](@ref) corresponds to ``\Gamma(a) = \Gamma(a,0)``. See also the [`gamma_inc`](@ref) function to compute both the upper and lower -(``\\gamma(a,x)``) incomplete gamma functions scaled by ``\\Gamma(a)``. +(``\gamma(a,x)``) incomplete gamma functions scaled by ``\Gamma(a)``. -External links: [DLMF](https://dlmf.nist.gov/8.2.2), [Wikipedia](https://en.wikipedia.org/wiki/Incomplete_gamma_function) +External links: +[DLMF 8.2.2](https://dlmf.nist.gov/8.2.2), +[Wikipedia](https://en.wikipedia.org/wiki/Incomplete_gamma_function) """ -gamma(a::Number,x::Number) = _gamma(promotereal(float(a), float(x))...) -gamma(a::Integer,x::Number) = _gamma(a, float(x)) +gamma(a::Number, x::Number) = _gamma(promotereal(float(a), float(x))...) +gamma(a::Integer, x::Number) = _gamma(a, float(x)) function _gamma(a::Number, x::Number) if a isa Real && x isa Real && !isfinite(a*x) @@ -1155,17 +1180,17 @@ function _gamma_big(a::Real,x::Real) return z end -""" +@doc raw""" loggamma(a,x) Returns the log of the upper incomplete gamma function [`gamma(a,x)`](@ref): ```math -\\log \\Gamma(a,x) = \\log \\int_x^\\infty t^{a-1} e^{-t} dt \\, +\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`` +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)`](@ref). diff --git a/src/logabsgamma/e_lgamma_r.jl b/src/logabsgamma/e_lgamma_r.jl index 639cca25..313a8ccf 100644 --- a/src/logabsgamma/e_lgamma_r.jl +++ b/src/logabsgamma/e_lgamma_r.jl @@ -145,12 +145,12 @@ function _logabsgamma(x::Float64) y = x - c p1 = y * @evalpoly(y, -7.72156649015328655494e-02, 6.32827064025093366517e-01, 1.45492250137234768737, 9.77717527963372745603e-01, 2.28963728064692451092e-01, 1.33810918536787660377e-02) p2 = @evalpoly(y, 1.0, 2.45597793713041134822, 2.12848976379893395361, 7.69285150456672783825e-01, 1.04222645593369134254e-01, 3.21709242282423911810e-03) - r += muladd(y, -0.5, p1 / p2) + r += muladd(y, -0.5, p1 / p2) end elseif ix < 0x40200000 #= x < 8.0 =# i = round(x, RoundToZero) y = x - i - z = 1.0 + z = 1.0 p = 0.0 u = x while u ≥ 3.0 @@ -165,7 +165,7 @@ function _logabsgamma(x::Float64) z = 1.0 / x y = z * z w = muladd(z, @evalpoly(y, 8.33333333333329678849e-2, -2.77777777728775536470e-3, 7.93650558643019558500e-4, -5.95187557450339963135e-4, 8.36339918996282139126e-4, -1.63092934096575273989e-3), 4.18938533204672725052e-1) - r = muladd(x - 0.5, log(x) - 1.0, w) + r = muladd(x - 0.5, log(x) - 1.0, w) else #= 2^58 ≤ x ≤ Inf =# r = muladd(x, log(x), -x) end diff --git a/src/logabsgamma/e_lgammaf_r.jl b/src/logabsgamma/e_lgammaf_r.jl index 88389652..ee8208ff 100644 --- a/src/logabsgamma/e_lgammaf_r.jl +++ b/src/logabsgamma/e_lgammaf_r.jl @@ -89,7 +89,7 @@ function _logabsgamma(x::Float32) end p = y * @evalpoly(y, -7.7215664089f-2, 2.1498242021f-1, 3.2577878237f-1, 1.4635047317f-1, 2.6642270386f-2, 1.8402845599f-3, 3.1947532989f-5) q = @evalpoly(y, 1.0f0, 1.3920053244f0, 7.2193557024f-1, 1.7193385959f-1, 1.8645919859f-2, 7.7794247773f-4, 7.3266842264f-6) - r = log(z) + muladd(0.5f0, y, p / q) + r = log(z) + muladd(0.5f0, y, p / q) elseif ix < 0x5c800000 #= 8.0 ≤ x < 2^58 =# z = 1.0f0 / x y = z * z diff --git a/src/sincosint.jl b/src/sincosint.jl index 783ec6b0..c4f00efd 100644 --- a/src/sincosint.jl +++ b/src/sincosint.jl @@ -238,7 +238,8 @@ Compute the sine integral function of ``x``, defined by x \in \mathbb{R} \,. ``` -External links: [DLMF](https://dlmf.nist.gov/6.2.9), +External links: +[DLMF 6.2.9](https://dlmf.nist.gov/6.2.9), [Wikipedia](https://en.wikipedia.org/wiki/Trigonometric_integral#Sine_integral). See also: [`cosint(x)`](@ref SpecialFunctions.cosint). @@ -270,7 +271,8 @@ x > 0 \,, ``` where ``\gamma`` is the Euler-Mascheroni constant. -External links: [DLMF](https://dlmf.nist.gov/6.2.13), +External links: +[DLMF 6.2.13](https://dlmf.nist.gov/6.2.13), [Wikipedia](https://en.wikipedia.org/wiki/Trigonometric_integral#Cosine_integral). See also: [`sinint(x)`](@ref SpecialFunctions.sinint).