From 1db90d6102c460be5b15989d46f6e17c4283368b Mon Sep 17 00:00:00 2001 From: Abhro R <5664668+abhro@users.noreply.github.com> Date: Sat, 15 Mar 2025 15:41:15 -0400 Subject: [PATCH] Replace TeX escape sequences with unicode counterparts The default KaTeX renderer supports it, and it also clears up the docstrings when viewing in the REPL --- docs/src/functions_overview.md | 16 +++---- src/bessel.jl | 26 ++++++------ src/beta_inc.jl | 16 +++---- src/betanc.jl | 22 +++++----- src/ellip.jl | 20 ++++----- src/erf.jl | 8 ++-- src/expint.jl | 12 +++--- src/gamma.jl | 26 +++++------- src/gamma_inc.jl | 78 +++++++++++++++++----------------- src/sincosint.jl | 4 +- 10 files changed, 112 insertions(+), 116 deletions(-) diff --git a/docs/src/functions_overview.md b/docs/src/functions_overview.md index ca22f5ac..ced6da5f 100644 --- a/docs/src/functions_overview.md +++ b/docs/src/functions_overview.md @@ -8,7 +8,7 @@ Here the *Special Functions* are listed according to the structure of [NIST Digi | Function | Description | |:-------- |:----------- | -| [`gamma(z)`](@ref SpecialFunctions.gamma(::Number)) | [gamma function](https://en.wikipedia.org/wiki/Gamma_function) ``\Gamma(z)`` | +| [`gamma(z)`](@ref SpecialFunctions.gamma(::Number)) | [gamma function](https://en.wikipedia.org/wiki/Gamma_function) ``Γ(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 | @@ -16,7 +16,7 @@ Here the *Special Functions* are listed according to the structure of [NIST Digi | [`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) | +| [`gamma(a,z)`](@ref SpecialFunctions.gamma(::Number,::Number)) | [upper incomplete gamma function ``Γ(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) | @@ -34,11 +34,11 @@ Here the *Special Functions* are listed according to the structure of [NIST Digi | 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}_ν(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}_ν(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 @@ -56,7 +56,7 @@ Here the *Special Functions* are listed according to the structure of [NIST Digi | [`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`` | +| [`dawson(x)`](@ref SpecialFunctions.dawson) | scaled imaginary error function, a.k.a. Dawson function, i.e. accurate ``\frac{\sqrt{π}}{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)`` | diff --git a/src/bessel.jl b/src/bessel.jl index a068c722..62a72ce8 100644 --- a/src/bessel.jl +++ b/src/bessel.jl @@ -465,7 +465,7 @@ 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_ν(x)``. External links: [DLMF 10.25.2](https://dlmf.nist.gov/10.25.2), @@ -486,7 +486,7 @@ 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_ν(x) e^{- | \operatorname{Re}(x) |}``. External links: [DLMF 10.25.2](https://dlmf.nist.gov/10.25.2), @@ -507,7 +507,7 @@ 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_ν(x)``. External links: [DLMF 10.2.2](https://dlmf.nist.gov/10.2.2), @@ -532,7 +532,7 @@ 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_ν(x) e^{- | \operatorname{Im}(x) |}``. External links: [DLMF 10.2.2](https://dlmf.nist.gov/10.2.2), @@ -553,7 +553,7 @@ 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_ν(x)``. External links: [DLMF 10.25.3](https://dlmf.nist.gov/10.25.3), @@ -576,7 +576,7 @@ 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_ν(x) e^x``. External links: [DLMF 10.25.3](https://dlmf.nist.gov/10.25.3), @@ -599,7 +599,7 @@ end """ bessely(nu, x) -Bessel function of the second kind of order `nu`, ``Y_\\nu(x)``. +Bessel function of the second kind of order `nu`, ``Y_ν(x)``. External links: [DLMF 10.2.3](https://dlmf.nist.gov/10.2.3), @@ -620,7 +620,7 @@ end besselyx(nu, x) Scaled Bessel function of the second kind of order `nu`, -``Y_\\nu(x) e^{- | \\operatorname{Im}(x) |}``. +``Y_ν(x) e^{- | \\operatorname{Im}(x) |}``. External links: [DLMF 10.2.3](https://dlmf.nist.gov/10.2.3), @@ -784,7 +784,7 @@ sphericalbessely(nu, x::T) where {T} = √((float(T))(π)/2x) * bessely(nu + one """ hankelh1(nu, x) -Bessel function of the third kind of order `nu`, ``H^{(1)}_\\nu(x)``. +Bessel function of the third kind of order `nu`, ``H^{(1)}_ν(x)``. External links: [DLMF 10.2.5](https://dlmf.nist.gov/10.2.5), @@ -797,7 +797,7 @@ hankelh1(nu, z) = besselh(nu, 1, z) """ hankelh2(nu, x) -Bessel function of the third kind of order `nu`, ``H^{(2)}_\\nu(x)``. +Bessel function of the third kind of order `nu`, ``H^{(2)}_ν(x)``. External links: [DLMF 10.2.6](https://dlmf.nist.gov/10.2.6), @@ -810,7 +810,7 @@ hankelh2(nu, z) = besselh(nu, 2, z) """ hankelh1x(nu, x) -Scaled Bessel function of the third kind of order `nu`, ``H^{(1)}_\\nu(x) e^{-x i}``. +Scaled Bessel function of the third kind of order `nu`, ``H^{(1)}_ν(x) e^{-x i}``. External links: [DLMF 10.2.5](https://dlmf.nist.gov/10.2.5), @@ -823,7 +823,7 @@ 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)}_ν(x) e^{x i}``. External links: [DLMF 10.2.6](https://dlmf.nist.gov/10.2.6), @@ -839,7 +839,7 @@ hankelh2x(nu, z) = besselhx(nu, 2, z) Bessel function of the first kind divided by `x`. Following convention: ```math -\operatorname{jinc}{x} = \frac{2 J_1({\pi x})}{\pi x}. +\operatorname{jinc}{x} = \frac{2 J_1({π x})}{π x}. ``` Sometimes known as sombrero or besinc function. diff --git a/src/beta_inc.jl b/src/beta_inc.jl index c6f09e03..ffb6175d 100644 --- a/src/beta_inc.jl +++ b/src/beta_inc.jl @@ -7,7 +7,7 @@ const exparg_p = log(prevfloat(floatmax(Float64))) @doc raw""" loggammadiv(a,b) -Computes ``\log(\Gamma(b)/\Gamma(a+b))`` when `b >= 8` +Computes ``\log(Γ(b)/Γ(a+b))`` when `b >= 8` """ loggammadiv(a::Number, b::Number) = _loggammadiv(promote(float(a), float(b))...) @@ -76,7 +76,7 @@ end @doc raw""" esum(mu,x) -Compute ``e^{\mu+x}`` +Compute ``e^{μ+x}`` """ function esum(mu::Float64, x::Float64) if x > 0.0 @@ -95,7 +95,7 @@ end @doc raw""" beta_integrand(a, b, x, y, mu=0.0) -Compute ``e^{\mu} x^a y^b / B(a,b)`` +Compute ``e^{μ} 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) @@ -181,7 +181,7 @@ end 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 ``λ = (a+b)*y - b`` External links: [DLMF 8.17.22](https://dlmf.nist.gov/8.17.22), @@ -251,7 +251,7 @@ end 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 ``λ = (a+b)*y - b``. External links: [DLMF 8.17.22](https://dlmf.nist.gov/8.17.22), @@ -525,7 +525,7 @@ end Computes ``I_x(a,b)`` using power series: ```math -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] +I_{x}(a,b) = G(a,b) x^{a}/a \left[1 + a \sum_{j=1}^∞ ((1-b)(2-b)\dots(j-b)/j!(a+j)) x^{j}\right] ``` External links: [DLMF 8.17.22](https://dlmf.nist.gov/8.17.22), @@ -733,9 +733,9 @@ end 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} \mathrm{d}t, +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) = Γ(a)Γ(b)/Γ(a+b)``. External links: [DLMF 8.17.1](https://dlmf.nist.gov/8.17.1), diff --git a/src/betanc.jl b/src/betanc.jl index 7b44f28d..2e3a7d89 100644 --- a/src/betanc.jl +++ b/src/betanc.jl @@ -10,9 +10,9 @@ const errmax = 1e-15 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) - Γ(a+b)/Γ(a+1)Γ(b) x^a (1-x)^b ``` -and ``\Gamma(a+1) = a\Gamma(a)`` given in [DLMF 8.17.21](https://dlmf.nist.gov/8.17.21). +and ``Γ(a+1) = aΓ(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 @@ -55,16 +55,16 @@ end 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)``. +First ``λ/2`` is calculated and the Poisson term is calculated using +``P(j-1) = j/λ P(j)`` and ``P(j+1) = λ/(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) + Γ(a+b+j-1)/Γ(a+j)Γ(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) - Γ(a+b+j)/Γ(a+j)Γ(b) x^{a+j} (1-x)^{b} ``` """ function ncbeta_poisson(a::Float64, b::Float64, lambda::Float64, x::Float64) @@ -147,10 +147,10 @@ end 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; λ) = \sum_{j=0}^∞ q(λ/2,j) I_{x}(a+j,b;0) ``` -For ``\lambda < 54`` : algorithm suggested by Lenth(1987) in `ncbeta_tail(a,b,lambda,x)`. -Else for ``\lambda \geq 54``: modification in Chattamvelli(1997) in +For ``λ < 54`` : algorithm suggested by Lenth(1987) in `ncbeta_tail(a,b,lambda,x)`. +Else for ``λ ≥ 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) @@ -173,9 +173,9 @@ end Compute CDF of noncentral F distribution given by: ```math -F(x, v_1, v_2; \lambda) = I_{v_1 x/(v_1 x + v_2)}(v_1/2, v_2/2; \lambda) +F(x, v_1, v_2; λ) = I_{v_1 x/(v_1 x + v_2)}(v_1/2, v_2/2; λ) ``` -where ``I_{x}(a,b; \lambda)`` is the noncentral beta function computed above. +where ``I_{x}(a,b; λ)`` is the noncentral beta function computed above. External links: [Wikipedia](https://en.wikipedia.org/wiki/Noncentral_F-distribution) diff --git a/src/ellip.jl b/src/ellip.jl index 7636fcd7..95db36b7 100644 --- a/src/ellip.jl +++ b/src/ellip.jl @@ -11,8 +11,8 @@ Computes Complete Elliptic Integral of 1st kind ``K(m)`` for parameter ``m`` giv ```math \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] \, . += \int_0^{ \frac{π}{2} } \frac{1}{\sqrt{1 - m \sin^2 θ}} \, \mathrm{d}θ +\quad \text{for} \quad m \in \left( -∞, 1 \right] \, . ``` External links: @@ -22,9 +22,9 @@ External links: See also: [`ellipe(m)`](@ref SpecialFunctions.ellipe). # Arguments -- `m`: parameter ``m``, restricted to the domain ``(-\infty,1]``, is related to +- `m`: parameter ``m``, restricted to the domain ``(-∞,1]``, is related to the elliptic modulus ``k`` by ``k^2=m`` and to the modular angle - ``\alpha`` by ``k = \sin \alpha``. + ``α`` by ``k = \sin α``. # Implementation Using piecewise approximation polynomial as given in @@ -40,7 +40,7 @@ For ``m<0``, followed by > Journal of Computational and Applied Mathematics. 282. > DOI 10.13140/2.1.1946.6245., > -As suggested in this paper, the domain is restricted to ``(-\infty,1]``. +As suggested in this paper, the domain is restricted to ``(-∞,1]``. """ ellipk(m::Real) = _ellipk(float(m)) @@ -191,8 +191,8 @@ Computes Complete Elliptic Integral of 2nd kind ``E(m)`` for parameter ``m`` giv ```math \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] . += \int_0^{ \frac{π}{2} } \sqrt{1 - m \sin^2 θ} \, \mathrm{d}θ +\quad \text{for} \quad m \in \left( -∞, 1 \right] . ``` External links: @@ -202,9 +202,9 @@ External links: See also: [`ellipk(m)`](@ref SpecialFunctions.ellipk). # Arguments -- `m`: parameter ``m``, restricted to the domain ``(-\infty,1]``, is related to +- `m`: parameter ``m``, restricted to the domain ``(-∞,1]``, is related to the elliptic modulus ``k`` by ``k^2=m`` and to the modular angle - ``\alpha`` by ``k=\sin \alpha``. + ``α`` by ``k = \sin α``. # Implementation Using piecewise approximation polynomial as given in @@ -220,7 +220,7 @@ For ``m<0``, followed by > Journal of Computational and Applied Mathematics. 282. > DOI 10.13140/2.1.1946.6245., > -As suggested in this paper, the domain is restricted to ``(-\infty,1]``. +As suggested in this paper, the domain is restricted to ``(-∞,1]``. """ ellipe(m::Real) = _ellipe(float(m)) diff --git a/src/erf.jl b/src/erf.jl index c878cf5d..a3b4581b 100644 --- a/src/erf.jl +++ b/src/erf.jl @@ -82,7 +82,7 @@ end Compute the error function of ``x``, defined by ```math -\operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x \exp(-t^2) \; \mathrm{d}t +\operatorname{erf}(x) = \frac{2}{\sqrt{π}} \int_0^x \exp(-t^2) \; \mathrm{d}t \quad \text{for} \quad x \in \mathbb{C} \, . ``` @@ -126,7 +126,7 @@ Compute the complementary error function of ``x``, defined by ```math \operatorname{erfc}(x) = 1 - \operatorname{erf}(x) -= \frac{2}{\sqrt{\pi}} \int_x^\infty \exp(-t^2) \; \mathrm{d}t += \frac{2}{\sqrt{π}} \int_x^∞ \exp(-t^2) \; \mathrm{d}t \quad \text{for} \quad x \in \mathbb{C} \, . ``` @@ -200,11 +200,11 @@ Compute the Dawson function (scaled imaginary error function) of ``x``, defined ```math \operatorname{dawson}(x) -= \frac{\sqrt{\pi}}{2} e^{-x^2} \operatorname{erfi}(x) += \frac{\sqrt{π}}{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)`` +This is the accurate version of ``\frac{\sqrt{π}}{2} e^{-x^2} \operatorname{erfi}(x)`` for large ``x``. External links: [DLMF 7.2.5](https://dlmf.nist.gov/7.2.5), diff --git a/src/expint.jl b/src/expint.jl index 63a86b19..8400d78d 100644 --- a/src/expint.jl +++ b/src/expint.jl @@ -510,9 +510,9 @@ end Computes the exponential integral ```math -\operatorname{E}_\nu(z) = \int_1^\infty \frac{e^{-zt}}{t^\nu} \mathrm{d}t. +\operatorname{E}_ν(z) = \int_1^∞ \frac{e^{-zt}}{t^ν} \mathrm{d}t. ``` -If ``\nu`` is not specified, ``\nu=1`` is used. Arbitrary complex ``\nu`` and ``z`` are supported. +If ``ν`` is not specified, ``ν=1`` is used. Arbitrary complex ``ν`` and ``z`` are supported. External links: [DLMF 8.19](https://dlmf.nist.gov/8.19), @@ -527,10 +527,10 @@ expint(ν::Number, z::Number, niter::Int=1000) = _expint(ν, z, niter, Val{false 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. +\exp(z) \operatorname{E}_ν(z) = e^z \int_1^∞ \frac{e^{-zt}}{t^ν} \mathrm{d}t. ``` -If ``\nu`` is not specified, ``\nu = 1`` is used. Arbitrary complex -``\nu`` and ``z`` are supported. +If ``ν`` is not specified, ``ν = 1`` is used. Arbitrary complex +``ν`` and ``z`` are supported. See also: [`expint(ν, z)`](@ref SpecialFunctions.expint) """ @@ -544,7 +544,7 @@ expintx(ν::Number, z::Number, niter::Int=1000) = _expint(ν, z, niter, Val{true Computes the exponential integral function ```math -\operatorname{Ei}(x) = \int_{-\infty}^x \frac{e^t}{t} \mathrm{d}t, +\operatorname{Ei}(x) = \int_{-∞}^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. diff --git a/src/gamma.jl b/src/gamma.jl index 5004404d..c80561e8 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -225,13 +225,13 @@ end Generalized zeta function defined by ```math -\zeta(s, z) = \sum_{k=0}^\infty \frac{1}{((k+z)^2)^{s/2}}, +ζ(s, z) = \sum_{k=0}^∞ \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}``. +``\sum_{k=0}^∞ (k+z)^{-s}``. -The Riemann zeta function is recovered as ``\zeta(s) = \zeta(s,1)``. +The Riemann zeta function is recovered as ``ζ(s) = ζ(s,1)``. External links: [Riemann zeta function](https://en.wikipedia.org/wiki/Riemann_zeta_function), @@ -413,7 +413,7 @@ end Riemann zeta function ```math -\zeta(s) = \sum_{n=1}^\infty \frac{1}{n^s}\quad\text{for}\quad s\in\mathbb{C}. +ζ(s) = \sum_{n=1}^∞ \frac{1}{n^s} \quad \text{for}\quad s\in\mathbb{C}. ``` External links: [Wikipedia](https://en.wikipedia.org/wiki/Riemann_zeta_function) @@ -490,7 +490,7 @@ end Dirichlet eta function ```math -\eta(s) = \sum_{n=1}^\infty (-1)^{n-1} / n^{s}. +η(s) = \sum_{n=1}^∞ (-1)^{n-1} / n^{s}. ``` """ eta(s::Number) = _eta(float(s)) @@ -550,14 +550,10 @@ end Compute the gamma function for complex ``z``, defined by ```math -\Gamma(z) -:= -\begin{cases} - n! - & \text{for} \quad z = n+1 \;, n = 0,1,2,\dots +Γ(z) := \begin{cases} + n! & \text{for} \quad z = n+1 \;, n = 0,1,2,\dots \\ - \int_0^\infty t^{z-1} e^{-t} \, \mathrm{d}t - & \text{for} \quad \Re(z) > 0 + \int_0^∞ t^{z-1} e^{-t} \, \mathrm{d}t & \text{for} \quad \Re(z) > 0 \end{cases} ``` and by analytic continuation in the whole complex plane. @@ -586,9 +582,9 @@ External links: [Wikipedia](https://en.wikipedia.org/wiki/Gamma_function). See also: -[`loggamma(z)`](@ref SpecialFunctions.loggamma) for ``\log \Gamma(z)`` and +[`loggamma(z)`](@ref SpecialFunctions.loggamma) for ``\log Γ(z)`` and [`gamma(a,z)`](@ref SpecialFunctions.gamma(::Number,::Number)) for -the upper incomplete gamma function ``\Gamma(a,z)``. +the upper incomplete gamma function ``Γ(a,z)``. # Implementation by - `Float`: C standard math library @@ -769,7 +765,7 @@ 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) = Γ(x)Γ(y)/Γ(x+y)``. """ function beta(a::Number, b::Number) lab, sign = logabsbeta(a, b) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index 7fbc9fb0..d6858657 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -33,7 +33,7 @@ const d80 = -.652623918595309E-03 @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/Γ(a+1) - 1`` for `-0.5<=a<=1.5`: ``1/Γ(a+1) - 1``. Uses the relation `gamma(a+1) = a*gamma(a)`. """ @@ -60,7 +60,7 @@ 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/Γ(a) e^{-x} x^{a}``. Based on `DRCOMP` from the `NSWC` library. """ function rgammax(a::Float64, x::Float64) if x == 0.0 @@ -85,7 +85,7 @@ 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/Γ(x+1) = 1 + x (x-1) g(x)``, `-1 <= x <= 1`. """ function auxgam(x::Float64) @assert -1 <= x <= 1 @@ -100,7 +100,7 @@ end @doc raw""" loggamma1p(x) -Compute ``\log(\Gamma(1+x))`` for `-1 < x <= 1`. +Compute ``\log(Γ(1+x))`` for `-1 < x <= 1`. """ function loggamma1p(x::Float64) @assert -1 < x <= 1 @@ -134,7 +134,7 @@ 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{Γ(x)} - (x-0.5)*\ln{x} + x - \ln(2π) / 2``. Adapted from `stirling` in `IncgamFI`. """ function stirling_error(x::Float64) if x < floatmin(Float64)*1000.0 @@ -164,7 +164,7 @@ end ```math \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. +\dfrac{Γ(x)}{\sqrt{2 π}e^{-x + (x-0.5)\log(x)}},\quad \text{if} \quad x\leq 0. \end{cases} ``` """ @@ -181,7 +181,7 @@ end @doc raw""" lambdaeta(eta) -Compute the value of ``\lambda`` satisfying ``\eta^{2}/2 = \lambda-1-\log{\lambda}``. +Compute the value of ``λ`` satisfying ``η^{2}/2 = λ - 1 - \log λ``. """ function lambdaeta(eta::Float64) s = eta*eta*0.5 @@ -222,7 +222,7 @@ end @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} +ε (η_0,a) = ε_1(η_0,a)a / a + ε_2(η_0,a) / a^2 + ε_3(η_0,a) / a^3 ``` Refer Eqn (3.12) in the paper """ @@ -247,7 +247,7 @@ end @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} +ε(η_0,a) = ε_1(η_0,a)/a + ε_2(η_0,a)/a^2 + ε_3(η_0,a)/a^3 ``` Refer Eqn (3.12) in the paper """ @@ -296,7 +296,7 @@ end @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} +ε(η_0,a) = ε_1(η_0,a)/a + ε_2(η_0,a)/a^2 + ε_3(η_0,a)/a^3 ``` Refer Eqn (3.12) in the paper """ @@ -416,7 +416,7 @@ end 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}^∞ x^n / ((a+1)(a+2)\dots(a+n))) ``` Used when `1 <= a <= BIG` and `x <= max{a, ln 10}`. @@ -457,7 +457,7 @@ end 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} + Θ_n a_n / x^n)) ``` where `R(a,x) = rgammax(a,x)`. Used when `1 <= a <= BIG` and `x >= x0`. @@ -497,11 +497,11 @@ end 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^∞ (-x)^{n}/((a+n)n!) ``` and ``P(a,x)/x^a`` is given by: ```math -(1 - J) / \Gamma(a+1) +(1 - J) / Γ(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. @@ -544,10 +544,10 @@ end Compute ``P(a,x)`` using minimax approximations given by : ```math -1/2 * \operatorname{erfc}(\sqrt{y}) - e^{-y}/\sqrt{2\pi a} ⋅ T(a,\lambda) +1/2 * \operatorname{erfc}(\sqrt{y}) - e^{-y}/\sqrt{2π a} ⋅ T(a,λ) ``` where ```math -T(a,\lambda) = \sum_{0}^{N} c_{k}(z)a^{-k} +T(a,λ) = \sum_0^N c_{k}(z)a^{-k} ``` This is a higher accuracy approximation of Temme expansion, which deals with @@ -609,12 +609,12 @@ end @doc raw""" gamma_inc_temme(a, x, z) -Compute ``P(a,x)`` using Temme's expansion given by : +Compute ``P(a,x)`` using Temme's expansion given by: ```math -1/2 * \operatorname{erfc}(\sqrt{y}) - e^{-y}/\sqrt{2\pi a} ⋅ T(a,\lambda) +1/2 * \operatorname{erfc}(\sqrt{y}) - e^{-y}/\sqrt{2π a} ⋅ T(a,λ) ``` where ```math -T(a,\lambda) = \sum_{0}^{N} c_{k}(z)a^{-k} +T(a,λ) = \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. @@ -646,13 +646,13 @@ end 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π a} ⋅ T(a,λ) ``` where ```math -E(y) = 1/2 - (1 - y/3) ⋅ \sqrt{y/\pi} +E(y) = 1/2 - (1 - y/3) ⋅ \sqrt{y/π} ``` -Used instead of it's previous function when ``\sigma \leq e_{0}/\sqrt{a}``. +Used instead of it's previous function when ``σ ≤ e_0/\sqrt{a}``. External links: [DLMF 8.12.8](https://dlmf.nist.gov/8.12.8) """ @@ -738,10 +738,10 @@ end 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}^∞ \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 Γ(1+a))^{1/a}`` +Inverting this relation we obtain ``x = r + \sum_{k=2}^∞ c_k r^k``. """ function gamma_inc_inv_psmall(a::Float64, logr::Float64) r = exp(logr) @@ -763,10 +763,10 @@ 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 Γ(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}^∞ d_{k}/x_{0}^{k} ``` where ``b = 1-a``, ``L = \ln{x_0}``. """ @@ -800,12 +800,12 @@ end Compute `x0` - initial approximation when `a` is large. The inversion problem is rewritten as: ```math -0.5 \operatorname{erfc}(\eta \sqrt{a/2}) + R_{a}(\eta) = q +0.5 \operatorname{erfc}(η \sqrt{a/2}) + R_{a}(η) = 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: ``η(q,a) = η_0(q,a) + ε(η_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} + \cdots +ε(η_0,a) = ε_1(η_0,a)/a + ε_2(η_0,a)/a^2 + ε_3(η_0,a)/a^3 + \cdots ``` 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 @@ -831,15 +831,15 @@ end Returns a tuple ``(p, q)`` where ``p + q = 1``, and ``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} \mathrm{d}t. +P(a,x) = \frac{1}{Γ(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: ```math -Q(a,x) = \frac{1}{\Gamma (a)} \int_{x}^{\infty} e^{-t}t^{a-1} \mathrm{d}t. +Q(a,x) = \frac{1}{Γ(a)} \int_x^∞ 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)``. +``γ(a,x) = P(a,x) Γ(a)`` and the upper incomplete +gamma function is ``Γ(a,x) = Q(a,x) Γ(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. @@ -1104,13 +1104,13 @@ promotereal(x, y) = promote(x,y) Returns the upper incomplete gamma function ```math -\Gamma(a,x) = \int_x^\infty t^{a-1} e^{-t} \mathrm{d}t +Γ(a,x) = \int_x^∞ 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 ``Γ(a) = Γ(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)``. +(``γ(a,x)``) incomplete gamma functions scaled by ``Γ(a)``. External links: [DLMF 8.2.2](https://dlmf.nist.gov/8.2.2), @@ -1165,12 +1165,12 @@ end 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} \mathrm{d}t +\log Γ(a,x) = \log \int_x^∞ 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πi`` (i.e. it may employ a different branch cut). See also [`loggamma(x)`](@ref). diff --git a/src/sincosint.jl b/src/sincosint.jl index c4f00efd..cbecf599 100644 --- a/src/sincosint.jl +++ b/src/sincosint.jl @@ -265,11 +265,11 @@ Compute the cosine integral function of ``x``, defined by ```math \operatorname{Ci}(x) -:= \gamma + \log x + \int_0^x \frac{\cos (t) - 1}{t} \, \mathrm{d}t +:= γ + \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. +where ``γ`` is the Euler-Mascheroni constant. External links: [DLMF 6.2.13](https://dlmf.nist.gov/6.2.13),