diff --git a/NEWS.md b/NEWS.md index bd7f17f..ee07c90 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,6 +10,11 @@ For bug fixes, performance enhancements, or fixes to unexported functions we wil # Unreleased +# Version 0.2.2 + +### Added + - Support for all airy functions and derivatives in entire complex plane ([PR #51](https://github.com/JuliaMath/Bessels.jl/pull/51)). These are specialized routines for airy functions instead of relying on connection to besselk. These are a couple digits more accurate and faster than previous versions. + # Version 0.2.1 ### Added diff --git a/Project.toml b/Project.toml index bf23bbf..264beb9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Bessels" uuid = "0e736298-9ec6-45e8-9647-e4fc86a2fe38" authors = ["Michael Helton and contributors"] -version = "0.2.1" +version = "0.2.2" [compat] julia = "1.6" diff --git a/src/Bessels.jl b/src/Bessels.jl index 736c5a5..63e7536 100644 --- a/src/Bessels.jl +++ b/src/Bessels.jl @@ -34,6 +34,8 @@ export airyaiprime export airybi export airybiprime +const ComplexOrReal{T} = Union{T,Complex{T}} + include("besseli.jl") include("besselj.jl") include("besselk.jl") diff --git a/src/airy.jl b/src/airy.jl index 1ed969f..1c8f920 100644 --- a/src/airy.jl +++ b/src/airy.jl @@ -1,36 +1,65 @@ # Airy functions # -# airyai(x), airybi(nu, x) -# airyaiprime(x), airybiprime(nu, x) +# airyai(z), airybi(nu, z) +# airyaiprime(z), airybiprime(nu, z) # -# A numerical routine to compute the airy functions and their derivatives. -# These routines use their relations to other special functions using https://dlmf.nist.gov/9.6. -# Specifically see (NIST 9.6.E1 - 9.6.E9) for computation from the defined bessel functions. -# For negative arguments these definitions are prone to some cancellation leading to higher errors. -# In the future, these could be replaced with more custom routines as they depend on a single variable. +# A numerical routine to compute the airy functions and their derivatives in the entire complex plane. +# These routines are based on the methods reported in [1] which use a combination of the power series +# for small arguments and a large argument expansion for (x > ~10). The primary difference between [1] +# and what is used here is that the regions where the power series and large argument expansions +# do not provide good results they are filled by relation to other special functions (besselk and besseli) +# using https://dlmf.nist.gov/9.6 (NIST 9.6.E1 - 9.6.E9). In this case the power series of besseli is used and then besselk +# is calculated using the continued fraction approach. This method is described in more detail in src/besselk.jl. +# However, care must be taken when computing besseli because when the imaginary component is much larger than the real part +# cancellation will occur. This can be overcome by shifting the order of besseli to be much larger and then using the power series +# and downward recurrence to get besseli(1/3, x). Another difficult region is when -10 zero(T) - return isinf(z) ? zero(T) : sqrt(x / 3) * besselk(one(T)/3, z) / π - elseif x < zero(T) - Jv, Yv = besseljy_positive_args(one(T)/3, z) - Jmv = (Jv - sqrt(T(3)) * Yv) / 2 - return isinf(z) ? throw(DomainError(x, "airyai(x) is not defined.")) : sqrt(x_abs) * (Jmv + Jv) / 3 - elseif iszero(x) - return T(0.3550280538878172) + # use relation to besselk (http://dlmf.nist.gov/9.6.E1) + zz = 2 * z * sqrt(z) / 3 + return sqrt(z / 3) * besselk_continued_fraction_shift(one(T)/3, zz) / T(π) + else + # z is close to the negative real axis + # for imag(z) == 0 use reflection to compute in terms of bessel functions of first kind (http://dlmf.nist.gov/9.6.E5) + # use computation for real numbers then convert to input type for stability + # for imag(z) != 0 use rotation identity (http://dlmf.nist.gov/9.2.E14) + if iszero(y) + xabs = abs(x) + xx = 2 * xabs * sqrt(xabs) / 3 + Jv, Yv = besseljy_positive_args(one(T)/3, xx) + Jmv = (Jv - sqrt(T(3)) * Yv) / 2 + return convert(eltype(z), sqrt(xabs) * (Jmv + Jv) / 3) + else + return cispi(one(T)/3) * _airyai(-z*cispi(one(T)/3)) + cispi(-one(T)/3) * _airyai(-z*cispi(-one(T)/3)) + end end end @@ -38,23 +67,41 @@ end airyaiprime(x) Derivative of the Airy function of the first kind ``\\operatorname{Ai}'(x)``. """ -airyaiprime(x::Real) = _airyaiprime(float(x)) +airyaiprime(z::Number) = _airyaiprime(float(z)) _airyaiprime(x::Float16) = Float16(_airyaiprime(Float32(x))) +_airyaiprime(z::ComplexF16) = ComplexF16(_airyaiprime(ComplexF32(z))) -function _airyaiprime(x::T) where T <: Union{Float32, Float64} - isnan(x) && return x - x_abs = abs(x) - z = 2 * x_abs^(T(3)/2) / 3 +function _airyaiprime(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} + if ~isfinite(z) + if abs(angle(z)) < 2*T(π)/3 + return -exp(-z) + else + return 1 / z + end + end + x, y = real(z), imag(z) + airy_large_argument_cutoff(z) && return airyaiprime_large_argument(z) + airyai_power_series_cutoff(x, y) && return airyaiprime_power_series(z) if x > zero(T) - return isinf(z) ? zero(T) : -x * besselk(T(2)/3, z) / (sqrt(T(3)) * π) - elseif x < zero(T) - Jv, Yv = besseljy_positive_args(T(2)/3, z) - Jmv = -(Jv + sqrt(T(3))*Yv) / 2 - return isinf(z) ? throw(DomainError(x, "airyaiprime(x) is not defined.")) : x_abs * (Jv - Jmv) / 3 - elseif iszero(x) - return T(-0.2588194037928068) + # use relation to besselk (http://dlmf.nist.gov/9.6.E2) + zz = 2 * z * sqrt(z) / 3 + return -z * besselk_continued_fraction_shift(T(2)/3, zz) / (T(π) * sqrt(T(3))) + else + # z is close to the negative real axis + # for imag(z) == 0 use reflection to compute in terms of bessel functions of first kind (http://dlmf.nist.gov/9.6.E5) + # use computation for real numbers then convert to input type for stability + # for imag(z) != 0 use rotation identity (http://dlmf.nist.gov/9.2.E14) + if iszero(y) + xabs = abs(x) + xx = 2 * xabs * sqrt(xabs) / 3 + Jv, Yv = besseljy_positive_args(T(2)/3, xx) + Jmv = -(Jv + sqrt(T(3))*Yv) / 2 + return convert(eltype(z), xabs * (Jv - Jmv) / 3) + else + return -(cispi(T(2)/3) * _airyaiprime(-z * cispi(one(T)/3)) + cispi(-T(2)/3) * _airyaiprime(-z * cispi(-one(T)/3))) + end end end @@ -62,23 +109,43 @@ end airybi(x) Airy function of the second kind ``\\operatorname{Bi}(x)``. """ -airybi(x::Real) = _airybi(float(x)) +airybi(z::Number) = _airybi(float(z)) _airybi(x::Float16) = Float16(_airybi(Float32(x))) +_airybi(z::ComplexF16) = ComplexF16(_airybi(ComplexF32(z))) -function _airybi(x::T) where T <: Union{Float32, Float64} - isnan(x) && return x - x_abs = abs(x) - z = 2 * x_abs^(T(3)/2) / 3 +function _airybi(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} + if ~isfinite(z) + if abs(angle(z)) < 2π/3 + return exp(z) + else + return 1 / z + end + end + x, y = real(z), imag(z) + airy_large_argument_cutoff(z) && return airybi_large_argument(z) + airybi_power_series_cutoff(x, y) && return airybi_power_series(z) if x > zero(T) - return isinf(z) ? z : sqrt(x / 3) * (besseli(-one(T)/3, z) + besseli(one(T)/3, z)) - elseif x < zero(T) - Jv, Yv = besseljy_positive_args(one(T)/3, z) - Jmv = (Jv - sqrt(T(3)) * Yv) / 2 - return isinf(z) ? throw(DomainError(x, "airybi(x) is not defined.")) : sqrt(x_abs/3) * (Jmv - Jv) - elseif iszero(x) - return T(0.6149266274460007) + zz = 2 * z * sqrt(z) / 3 + shift = 20 + order = one(T)/3 + inu, inum1 = besseli_power_series_inu_inum1(order + shift, zz) + inu, inum1 = besselk_down_recurrence(zz, inum1, inu, order + shift - 1, order) + + inu2, inum2 = besseli_power_series_inu_inum1(-order + shift, zz) + inu2, inum2 = besselk_down_recurrence(zz, inum2, inu2, -order + shift - 1, -order) + return sqrt(z/3) * (inu + inu2) + else + if iszero(y) + xabs = abs(x) + xx = 2 * xabs * sqrt(xabs) / 3 + Jv, Yv = besseljy_positive_args(one(T)/3, xx) + Jmv = (Jv - sqrt(T(3)) * Yv) / 2 + return convert(eltype(z), sqrt(xabs/3) * (Jmv - Jv)) + else + return cispi(one(T)/3) * _airybi(-z * cispi(one(T)/3)) + cispi(-one(T)/3) * _airybi(-z*cispi(-one(T)/3)) + end end end @@ -86,22 +153,357 @@ end airybiprime(x) Derivative of the Airy function of the second kind ``\\operatorname{Bi}'(x)``. """ -airybiprime(x::Real) = _airybiprime(float(x)) + +airybiprime(z::Number) = _airybiprime(float(z)) _airybiprime(x::Float16) = Float16(_airybiprime(Float32(x))) +_airybiprime(z::ComplexF16) = ComplexF16(_airybiprime(ComplexF32(z))) -function _airybiprime(x::T) where T <: Union{Float32, Float64} - isnan(x) && return x - x_abs = abs(x) - z = 2 * x_abs^(T(3)/2) / 3 +function _airybiprime(z::ComplexOrReal{T}) where T <: Union{Float32, Float64} + if ~isfinite(z) + if abs(angle(z)) < 2*T(π)/3 + return exp(z) + else + return -1 / z + end + end + x, y = real(z), imag(z) + airy_large_argument_cutoff(z) && return airybiprime_large_argument(z) + airybi_power_series_cutoff(x, y) && return airybiprime_power_series(z) if x > zero(T) - return isinf(z) ? z : x * (besseli(T(2)/3, z) + besseli(-T(2)/3, z)) / sqrt(T(3)) - elseif x < zero(T) - Jv, Yv = besseljy_positive_args(T(2)/3, z) - Jmv = -(Jv + sqrt(T(3))*Yv) / 2 - return isinf(z) ? throw(DomainError(x, "airybiprime(x) is not defined.")) : x_abs * (Jv + Jmv) / sqrt(T(3)) - elseif iszero(x) - return T(0.4482883573538264) + zz = 2 * z * sqrt(z) / 3 + shift = 20 + order = T(2)/3 + inu, inum1 = besseli_power_series_inu_inum1(order + shift, zz) + inu, inum1 = besselk_down_recurrence(zz, inum1, inu, order + shift - 1, order) + + inu2, inum2 = besseli_power_series_inu_inum1(-order + shift, zz) + inu2, inum2 = besselk_down_recurrence(zz, inum2, inu2, -order + shift - 1, -order) + return z / sqrt(3) * (inu + inu2) + else + if iszero(y) + xabs = abs(x) + xx = 2 * xabs * sqrt(xabs) / 3 + Jv, Yv = besseljy_positive_args(T(2)/3, xx) + Jmv = -(Jv + sqrt(T(3))*Yv) / 2 + return convert(eltype(z), xabs * (Jv + Jmv) / sqrt(T(3))) + else + return -(cispi(T(2)/3) * _airybiprime(-z*cispi(one(T)/3)) + cispi(-T(2)/3) * _airybiprime(-z*cispi(-one(T)/3))) + end end end + +##### +##### Power series for airyai(x) +##### + +# cutoffs for power series valid for both airyai and airyaiprime +airyai_power_series_cutoff(x::T, y::T) where T <: Float64 = x < 2 && abs(y) > -1.4*(x + 5.5) +airyai_power_series_cutoff(x::T, y::T) where T <: Float32 = x < 4.5f0 && abs(y) > -1.4f0*(x + 9.5f0) + +function airyai_power_series(x::ComplexOrReal{T}; tol=eps(T)) where T + S = eltype(x) + iszero(x) && return S(0.3550280538878172) + MaxIter = 3000 + ai1 = zero(S) + ai2 = zero(S) + x2 = x*x + x3 = x2*x + t = one(S) / GAMMA_TWO_THIRDS(T) + t2 = 3*x / GAMMA_ONE_THIRD(T) + + for i in 0:MaxIter + ai1 += t + ai2 += t2 + abs(t) < tol * abs(ai1) && break + t *= x3 * inv(9*(i + one(T))*(i + T(2)/3)) + t2 *= x3 * inv(9*(i + one(T))*(i + T(4)/3)) + end + return (ai1*3^(-T(2)/3) - ai2*3^(-T(4)/3)) +end +airyai_power_series(x::Float32) = Float32(airyai_power_series(Float64(x), tol=eps(Float32))) +airyai_power_series(x::ComplexF32) = ComplexF32(airyai_power_series(ComplexF64(x), tol=eps(Float32))) + +##### +##### Power series for airyaiprime(x) +##### + +function airyaiprime_power_series(x::ComplexOrReal{T}; tol=eps(T)) where T + S = eltype(x) + iszero(x) && return S(-0.2588194037928068) + MaxIter = 3000 + ai1 = zero(S) + ai2 = zero(S) + x2 = x*x + x3 = x2*x + t = one(S) / GAMMA_ONE_THIRD(T) + t2 = 3*x2 / (2*GAMMA_TWO_THIRDS(T)) + + for i in 0:MaxIter + ai1 += t + ai2 += t2 + abs(t) < tol * abs(ai1) && break + t *= x3 * inv(9*(i + one(T))*(i + T(1)/3)) + t2 *= x3 * inv(9*(i + one(T))*(i + T(5)/3)) + end + return -ai1*3^(-T(1)/3) + ai2*3^(-T(5)/3) +end +airyaiprime_power_series(x::Float32) = Float32(airyaiprime_power_series(Float64(x), tol=eps(Float32))) +airyaiprime_power_series(x::ComplexF32) = ComplexF32(airyaiprime_power_series(ComplexF64(x), tol=eps(Float32))) + +##### +##### Power series for airybi(x) +##### + +# cutoffs for power series valid for both airybi and airybiprime +# has a more complicated validity as it works well close to positive real line and for small negative arguments also works for angle(z) ~ 2pi/3 +# the statements are somewhat complicated but we want to hit this branch when we can as the other algorithms are 10x slower +# the Float32 cutoff can be simplified because it overlaps with the large argument expansion so there isn't a need for more complicated statements +airybi_power_series_cutoff(x::T, y::T) where T <: Float64 = (abs(y) > -1.4*(x + 5.5) && abs(y) < -2.2*(x - 4)) || (x > zero(T) && abs(y) < 3) +airybi_power_series_cutoff(x::T, y::T) where T <: Float32 = abs(complex(x, y)) < 9 + +function airybi_power_series(x::ComplexOrReal{T}; tol=eps(T)) where T + S = eltype(x) + iszero(x) && return S(0.6149266274460007) + MaxIter = 3000 + ai1 = zero(S) + ai2 = zero(S) + x2 = x*x + x3 = x2*x + t = one(S) / GAMMA_TWO_THIRDS(T) + t2 = 3*x / GAMMA_ONE_THIRD(T) + + for i in 0:MaxIter + ai1 += t + ai2 += t2 + abs(t) < tol * abs(ai1) && break + t *= x3 * inv(9*(i + one(T))*(i + T(2)/3)) + t2 *= x3 * inv(9*(i + one(T))*(i + T(4)/3)) + end + return (ai1*3^(-T(1)/6) + ai2*3^(-T(5)/6)) +end +airybi_power_series(x::Float32) = Float32(airybi_power_series(Float64(x), tol=eps(Float32))) +airybi_power_series(x::ComplexF32) = ComplexF32(airybi_power_series(ComplexF64(x), tol=eps(Float32))) + +##### +##### Power series for airybiprime(x) +##### + +function airybiprime_power_series(x::ComplexOrReal{T}; tol=eps(T)) where T + S = eltype(x) + iszero(x) && return S(0.4482883573538264) + MaxIter = 3000 + ai1 = zero(S) + ai2 = zero(S) + x2 = x*x + x3 = x2*x + t = one(S) / GAMMA_ONE_THIRD(T) + t2 = 3*x2 / (2*GAMMA_TWO_THIRDS(T)) + + for i in 0:MaxIter + ai1 += t + ai2 += t2 + abs(t) < tol * abs(ai1) && break + t *= x3 * inv(9*(i + one(T))*(i + T(1)/3)) + t2 *= x3 * inv(9*(i + one(T))*(i + T(5)/3)) + end + return (ai1*3^(T(1)/6) + ai2*3^(-T(7)/6)) +end +airybiprime_power_series(x::Float32) = Float32(airybiprime_power_series(Float64(x), tol=eps(Float32))) +airybiprime_power_series(x::ComplexF32) = ComplexF32(airybiprime_power_series(ComplexF64(x), tol=eps(Float32))) + +##### +##### Large argument expansion for airy functions +##### +airy_large_argument_cutoff(z::ComplexOrReal{Float64}) = abs(z) > 8 +airy_large_argument_cutoff(z::ComplexOrReal{Float32}) = abs(z) > 4 + +function airyai_large_argument(x::Real) + x < zero(x) && return real(airyai_large_argument(complex(x))) + return airy_large_arg_a(abs(x)) +end + +function airyai_large_argument(z::Complex{T}) where T + x, y = real(z), imag(z) + a = airy_large_arg_a(z) + if x < zero(T) && abs(y) < 5.5 + b = airy_large_arg_b(z) + y >= zero(T) ? (return a + im*b) : (return a - im*b) + end + return a +end + +function airyaiprime_large_argument(x::Real) + x < zero(x) && return real(airyaiprime_large_argument(complex(x))) + return airy_large_arg_c(abs(x)) +end + +function airyaiprime_large_argument(z::Complex{T}) where T + x, y = real(z), imag(z) + c = airy_large_arg_c(z) + if x < zero(T) && abs(y) < 5.5 + d = airy_large_arg_d(z) + y >= zero(T) ? (return c + im*d) : (return c - im*d) + end + return c +end + +function airybi_large_argument(x::Real) + if x < zero(x) + return 2*real(airy_large_arg_b(complex(x))) + else + return 2*(airy_large_arg_b(x)) + end +end + +function airybi_large_argument(z::Complex{T}) where T + x, y = real(z), imag(z) + b = airy_large_arg_b(z) + abs(y) <= 1.7*(x - 6) && return 2*b + + check_conj = false + if y < zero(T) + z = conj(z) + b = conj(b) + y = abs(y) + check_conj = true + end + + a = airy_large_arg_a(z) + if x < zero(T) && y < 5 + out = b + im*a + check_conj && (out = conj(out)) + return out + else + out = 2*b + im*a + check_conj && (out = conj(out)) + return out + end +end + +function airybiprime_large_argument(x::Real) + if x < zero(x) + return 2*real(airy_large_arg_d(complex(x))) + else + return 2*(airy_large_arg_d(x)) + end +end + +function airybiprime_large_argument(z::Complex{T}) where T + x, y = real(z), imag(z) + d = airy_large_arg_d(z) + abs(y) <= 1.7*(x - 6) && return 2*d + + check_conj = false + if y < zero(T) + z = conj(z) + d = conj(d) + y = abs(y) + check_conj = true + end + + c = airy_large_arg_c(z) + if x < zero(T) && y < 5 + out = d + im*c + check_conj && (out = conj(out)) + return out + else + out = 2*d + im*c + check_conj && (out = conj(out)) + return out + end +end + +# see equations 24 and relations using eq 25 and 26 in [1] +function airy_large_arg_a(x::ComplexOrReal{T}; tol=eps(T)*40) where T + S = eltype(x) + MaxIter = 3000 + xsqr = sqrt(x) + + out = zero(S) + t = GAMMA_ONE_SIXTH(T) * GAMMA_FIVE_SIXTHS(T) / 4 + a = 4*xsqr*x + for i in 0:MaxIter + out += t + abs(t) < tol*abs(out) && break + t *= -3*(i + one(T)/6) * (i + T(5)/6) / (a*(i + one(T))) + end + return out * exp(-a / 6) / (sqrt(T(π)^3) * sqrt(xsqr)) +end + +function airy_large_arg_b(x::ComplexOrReal{T}; tol=eps(T)*40) where T + S = eltype(x) + MaxIter = 3000 + xsqr = sqrt(x) + + out = zero(S) + t = GAMMA_ONE_SIXTH(T) * GAMMA_FIVE_SIXTHS(T) / 4 + a = 4*xsqr*x + for i in 0:MaxIter + out += t + abs(t) < tol*abs(out) && break + t *= 3*(i + one(T)/6) * (i + T(5)/6) / (a*(i + one(T))) + end + return out * exp(a / 6) / (sqrt(T(π)^3) * sqrt(xsqr)) +end + +function airy_large_arg_c(x::ComplexOrReal{T}; tol=eps(T)*40) where T + S = eltype(x) + MaxIter = 3000 + xsqr = sqrt(x) + + out = zero(S) + # use identities of gamma to relate to defined constants + # t = gamma(-one(T) / 6) * gamma(T(7) / 6) / 4 + t = -GAMMA_FIVE_SIXTHS(T) * GAMMA_ONE_SIXTH(T) / 4 + a = 4*xsqr*x + for i in 0:MaxIter + out += t + abs(t) < tol*abs(out) && break + t *= -3*(i - one(T)/6) * (i + T(7)/6) / (a*(i + one(T))) + end + return out * exp(-a / 6) * sqrt(xsqr) / sqrt(T(π)^3) +end + +function airy_large_arg_d(x::ComplexOrReal{T}; tol=eps(T)*40) where T + S = eltype(x) + MaxIter = 3000 + xsqr = sqrt(x) + + out = zero(S) + # use identities of gamma to relate to defined constants + # t = gamma(-one(T) / 6) * gamma(T(7) / 6) / 4 + t = -GAMMA_FIVE_SIXTHS(T) * GAMMA_ONE_SIXTH(T) / 4 + a = 4*xsqr*x + for i in 0:MaxIter + out += t + abs(t) < tol*abs(out) && break + t *= 3*(i - one(T)/6) * (i + T(7)/6) / (a*(i + one(T))) + end + return -out * exp(a / 6) * sqrt(xsqr) / sqrt(T(π)^3) +end + +# negative arguments of airy functions oscillate around zero so as x -> -Inf it is more prone to cancellation +# to give best accuracy it is best to promote to Float64 numbers until the Float32 tolerance +airy_large_arg_a(x::Float32) = (airy_large_arg_a(Float64(x), tol=eps(Float32))) +airy_large_arg_a(x::ComplexF32) = (airy_large_arg_a(ComplexF64(x), tol=eps(Float32))) + +airy_large_arg_b(x::Float32) = Float32(airy_large_arg_b(Float64(x), tol=eps(Float32))) +airy_large_arg_b(x::ComplexF32) = ComplexF32(airy_large_arg_b(ComplexF64(x), tol=eps(Float32))) + +airy_large_arg_c(x::Float32) = Float32(airy_large_arg_c(Float64(x), tol=eps(Float32))) +airy_large_arg_c(x::ComplexF32) = ComplexF32(airy_large_arg_c(ComplexF64(x), tol=eps(Float32))) + +airy_large_arg_d(x::Float32) = Float32(airy_large_arg_d(Float64(x), tol=eps(Float32))) +airy_large_arg_d(x::ComplexF32) = ComplexF32(airy_large_arg_d(ComplexF64(x), tol=eps(Float32))) +# calculates besselk from the power series of besseli using the continued fraction and wronskian +# this shift the order higher first to avoid cancellation in the power series of besseli along the imaginary axis +# for real arguments this is not needed because besseli can be computed stably over the entire real axis +function besselk_continued_fraction_shift(nu, x) + shift = 20 + inu, inum1 = besseli_power_series_inu_inum1(nu + shift, x) + inu, inum1 = besselk_down_recurrence(x, inum1, inu, nu + shift - 1, nu) + H_knu = besselk_ratio_knu_knup1(nu-1, x) + return 1 / (x * (inum1 + inu / H_knu)) +end diff --git a/src/besseli.jl b/src/besseli.jl index ce9c967..e4afe97 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -339,9 +339,10 @@ besseli_large_argument_cutoff(nu, x::Float32) = x > 18.0f0 && x > nu^2 / 19.5f0 Computes ``I_{nu}(x)`` using the power series for any value of nu. """ -function besseli_power_series(v, x::T) where T +function besseli_power_series(v, x::ComplexOrReal{T}) where T MaxIter = 3000 - out = zero(T) + S = eltype(x) + out = zero(S) xs = (x/2)^v a = xs / gamma(v + one(T)) t2 = (x/2)^2 diff --git a/src/besselk.jl b/src/besselk.jl index 0bfbe6b..eca9deb 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -339,10 +339,11 @@ end # a modified version of the I_{nu} power series to compute both I_{nu} and I_{nu-1} # use this along with the continued fractions for besselk -function besseli_power_series_inu_inum1(v, x::T) where T +function besseli_power_series_inu_inum1(v, x::ComplexOrReal{T}) where T MaxIter = 3000 - out = zero(T) - out2 = zero(T) + S = eltype(x) + out = zero(S) + out2 = zero(S) x2 = x / 2 xs = x2^v gmx = xs / gamma(v) @@ -361,11 +362,14 @@ end # computes K_{nu+1}/K_{nu} using continued fractions and the modified Lentz method # generally slow to converge for small x -function besselk_ratio_knu_knup1(v, x::T) where T +besselk_ratio_knu_knup1(v, x::Float32) = Float32(besselk_ratio_knu_knup1(v, Float64(x))) +besselk_ratio_knu_knup1(v, x::ComplexF32) = ComplexF32(besselk_ratio_knu_knup1(v, ComplexF64(x))) +function besselk_ratio_knu_knup1(v, x::ComplexOrReal{T}) where T MaxIter = 1000 - (hn, Dn, Cn) = (1e-50, zero(v), 1e-50) + S = eltype(x) + (hn, Dn, Cn) = (S(1e-50), zero(S), S(1e-50)) - jf = one(T) + jf = one(S) vv = v * v for _ in 1:MaxIter an = (vv - ((2*jf - 1)^2) * T(0.25)) @@ -398,9 +402,12 @@ Computes ``K_{nu}(x)`` using the power series when nu is not an integer. In general, this is most accurate for small arguments and when nu > x. No checks are performed on nu so this is not accurate when nu is an integer. """ -function besselk_power_series(v, x::T) where T +besselk_power_series(v, x::Float32) = Float32(besselk_power_series(v, Float64(x))) +besselk_power_series(v, x::ComplexF32) = ComplexF32(besselk_power_series(v, ComplexF64(x))) + +function besselk_power_series(v, x::ComplexOrReal{T}) where T MaxIter = 1000 - S = promote_type(T, Float64) + S = eltype(x) v, x = S(v), S(x) z = x / 2 @@ -431,7 +438,7 @@ function besselk_power_series(v, x::T) where T xd2_pow *= zz fact_k *= k + one(S) end - return T(out) + return out end besselk_power_series_cutoff(nu, x::Float64) = x < 2.0 || nu > 1.6x - 1.0 besselk_power_series_cutoff(nu, x::Float32) = x < 10.0f0 || nu > 1.65f0*x - 8.0f0 @@ -452,9 +459,11 @@ end besselk_large_argument_scaled(v, x::T) where T = T(_besselk_large_argument(v, x) * sqrt(pi / 2x)) -function _besselk_large_argument(v, x::T) where T +_besselk_large_argument(v, x::Float32) = Float32(_besselk_large_argument(v, Float64(x))) +_besselk_large_argument(v, x::ComplexF32) = ComplexF32(_besselk_large_argument(v, ComplexF64(x))) +function _besselk_large_argument(v, x::ComplexOrReal{T}) where T MaxIter = 5000 - S = promote_type(T, Float64) + S = eltype(x) v, x = S(v), S(x) fv2 = 4 * v^2 diff --git a/src/gamma.jl b/src/gamma.jl index 488ad00..25aa086 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -1,16 +1,18 @@ # Adapted from Cephes Mathematical Library (MIT license https://en.smath.com/view/CephesMathLibrary/license) by Stephen L. Moshier -function gamma(x) +gamma(z::Number) = _gamma(float(z)) +_gamma(x::Float32) = Float32(_gamma(Float64(x))) + +function _gamma(x::Float64) if x < 0 isinteger(x) && throw(DomainError(x, "NaN result for non-NaN input.")) xp1 = abs(x) + 1.0 - return π / sinpi(xp1) / _gamma(xp1) + return π / sinpi(xp1) / _gammax(xp1) else - return _gamma(x) + return _gammax(x) end end # only have a Float64 implementations -gamma(x::Float32) = Float32(gamma(Float64(x))) -function _gamma(x) +function _gammax(x) if x > 11.5 return large_gamma(x) elseif x <= 11.5 diff --git a/src/math_constants.jl b/src/math_constants.jl index f98e221..30d06ea 100644 --- a/src/math_constants.jl +++ b/src/math_constants.jl @@ -22,3 +22,13 @@ const TWOOPI(::Type{Float32}) = 0.636619772367581343075535f0 const THPIO4(::Type{Float32}) = 2.35619449019234492885f0 const SQ2O2(::Type{Float32}) = 0.7071067811865476f0 const SQPIO2(::Type{Float32}) = 1.25331413731550025f0 + +const GAMMA_TWO_THIRDS(::Type{Float64}) = 1.3541179394264005 +const GAMMA_ONE_THIRD(::Type{Float64}) = 2.6789385347077475 +const GAMMA_ONE_SIXTH(::Type{Float64}) = 5.566316001780235 +const GAMMA_FIVE_SIXTHS(::Type{Float64}) = 1.128787029908126 + +const GAMMA_TWO_THIRDS(::Type{Float32}) = 1.3541179394264005f0 +const GAMMA_ONE_THIRD(::Type{Float32}) = 2.6789385347077475f0 +const GAMMA_ONE_SIXTH(::Type{Float32}) = 5.566316001780235f0 +const GAMMA_FIVE_SIXTHS(::Type{Float32}) = 1.128787029908126f0 diff --git a/src/recurrence.jl b/src/recurrence.jl index c959537..4321cf8 100644 --- a/src/recurrence.jl +++ b/src/recurrence.jl @@ -41,8 +41,6 @@ end return jnup1, jnu end -#= -# currently not used # backward recurrence relation for besselk and besseli # outputs both (bessel(x, nu_end), bessel(x, nu_end-1) # x = 0.1; k0 = besseli(10,x); k1 = besseli(11,x); @@ -55,4 +53,3 @@ end end return jnup1, jnu end -=# diff --git a/test/airy_test.jl b/test/airy_test.jl index aa09707..e63f0b4 100644 --- a/test/airy_test.jl +++ b/test/airy_test.jl @@ -3,27 +3,65 @@ # this is amplified for negative arguments # the implementations provided by SpecialFunctions.jl suffer from same inaccuracies for x in [0.0; 1e-17:0.1:100.0] - @test isapprox(airyai(x), SpecialFunctions.airyai(x), rtol=1e-12) - @test isapprox(airyai(-x), SpecialFunctions.airyai(-x), rtol=1e-9) + @test isapprox(airyai(x), SpecialFunctions.airyai(x), rtol=2e-13) + @test isapprox(airyai(-x), SpecialFunctions.airyai(-x), rtol=3e-12) - @test isapprox(airyaiprime(x), SpecialFunctions.airyaiprime(x), rtol=1e-12) - @test isapprox(airyaiprime(-x), SpecialFunctions.airyaiprime(-x), rtol=1e-9) + @test isapprox(airyaiprime(x), SpecialFunctions.airyaiprime(x), rtol=2e-13) + @test isapprox(airyaiprime(-x), SpecialFunctions.airyaiprime(-x), rtol=5e-12) - @test isapprox(airybi(x), SpecialFunctions.airybi(x), rtol=1e-12) - @test isapprox(airybi(-x), SpecialFunctions.airybi(-x), rtol=5e-8) + @test isapprox(airybi(x), SpecialFunctions.airybi(x), rtol=2e-13) + @test isapprox(airybi(-x), SpecialFunctions.airybi(-x), rtol=5e-12) - @test isapprox(airybiprime(x), SpecialFunctions.airybiprime(x), rtol=1e-12) - @test isapprox(airybiprime(-x), SpecialFunctions.airybiprime(-x), rtol=1e-9) + @test isapprox(airybiprime(x), SpecialFunctions.airybiprime(x), rtol=2e-13) + @test isapprox(airybiprime(-x), SpecialFunctions.airybiprime(-x), rtol=5e-12) +end + +# Float32 +for x in [0.0; 0.5:0.5:30.0] + @test isapprox(airyai(x), SpecialFunctions.airyai(x), rtol=2e-13) + @test isapprox(airyai(-x), SpecialFunctions.airyai(-x), rtol=3e-12) + + @test isapprox(airyaiprime(x), SpecialFunctions.airyaiprime(x), rtol=2e-13) + @test isapprox(airyaiprime(-x), SpecialFunctions.airyaiprime(-x), rtol=5e-12) + + @test isapprox(airybi(x), SpecialFunctions.airybi(x), rtol=2e-13) + @test isapprox(airybi(-x), SpecialFunctions.airybi(-x), rtol=5e-12) + + @test isapprox(airybiprime(x), SpecialFunctions.airybiprime(x), rtol=2e-13) + @test isapprox(airybiprime(-x), SpecialFunctions.airybiprime(-x), rtol=5e-12) +end + +for x in [0.0, 0.01, 0.5, 1.0, 2.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 50.0], a in 0:pi/12:2pi + z = x*exp(im*a) + @test isapprox(airyai(z), SpecialFunctions.airyai(z), rtol=5e-10) + @test isapprox(airyaiprime(z), SpecialFunctions.airyaiprime(z), rtol=5e-10) + @test isapprox(airybi(z), SpecialFunctions.airybi(z), rtol=1e-11) + @test isapprox(airybiprime(z), SpecialFunctions.airybiprime(z), rtol=1e-11) end # test Inf + @test iszero(airyai(Inf)) @test iszero(airyaiprime(Inf)) @test isinf(airybi(Inf)) @test isinf(airybiprime(Inf)) +@test airyai(Inf + 0.0im) === exp(-(Inf + 0.0im)) +@test airyaiprime(Inf + 0.0im) === -exp(-(Inf + 0.0im)) +@test airyai(-Inf + 0.0im) === 1 / (-Inf + 0.0im) +@test airyaiprime(-Inf + 0.0im) === 1 / (-Inf + 0.0im) +@test airybi(Inf + Inf*im) === exp((Inf + Inf*im)) +@test airybi(-Inf + 10.0*im) === 1 / (-Inf + 10.0*im) +@test airybiprime(Inf + 0.0*im) === exp((Inf + 0.0*im)) +@test airybiprime(-Inf + 0.0*im) === -1 / (-Inf + 0.0*im) + + # test Float16 types @test airyai(Float16(1.2)) isa Float16 +@test airyai(ComplexF16(1.2)) isa ComplexF16 @test airyaiprime(Float16(1.9)) isa Float16 +@test airyaiprime(ComplexF16(1.2)) isa ComplexF16 @test airybi(Float16(1.2)) isa Float16 +@test airybi(ComplexF16(1.2)) isa ComplexF16 @test airybiprime(Float16(1.9)) isa Float16 +@test airybiprime(ComplexF16(1.2)) isa ComplexF16