diff --git a/src/Bessels.jl b/src/Bessels.jl index e11b9f5..736c5a5 100644 --- a/src/Bessels.jl +++ b/src/Bessels.jl @@ -41,6 +41,7 @@ include("bessely.jl") include("hankel.jl") include("airy.jl") include("sphericalbessel.jl") +include("modifiedsphericalbessel.jl") include("Float128/besseli.jl") include("Float128/besselj.jl") diff --git a/src/besselk.jl b/src/besselk.jl index 6ecdf42..b1e2833 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -225,9 +225,15 @@ function besselk_positive_args(nu, x::T) where T <: Union{Float32, Float64} # dispatch to avoid uniform expansion when nu = 0 iszero(nu) && return besselk0(x) + + # pre-compute the uniform asymptotic expansion cutoff: + debye_cut = besselik_debye_cutoff(nu, x) + + # check if nu is a half-integer: + (isinteger(nu-1/2) && !debye_cut) && return sphericalbesselk(nu-1/2, x)*SQRT_PID2(T)*sqrt(x) # use uniform debye expansion if x or nu is large - besselik_debye_cutoff(nu, x) && return besselk_large_orders(nu, x) + debye_cut && return besselk_large_orders(nu, x) # for integer nu use forward recurrence starting with K_0 and K_1 isinteger(nu) && return besselk_up_recurrence(x, besselk1(x), besselk0(x), 1, nu)[1] @@ -424,3 +430,4 @@ function besselk_power_series(v, x::T) where T 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 + diff --git a/src/constants.jl b/src/constants.jl index 820e8d1..8568346 100644 --- a/src/constants.jl +++ b/src/constants.jl @@ -284,3 +284,6 @@ const Q3_k1(::Type{Float64}) = ( 2.88448064302447607e1, 2.27912927104139732e0, 2.50358186953478678e-2 ) + +const SQRT_PID2(::Type{Float64}) = 1.2533141373155003 +const SQRT_PID2(::Type{Float32}) = 1.2533141f0 diff --git a/src/modifiedsphericalbessel.jl b/src/modifiedsphericalbessel.jl new file mode 100644 index 0000000..2a521af --- /dev/null +++ b/src/modifiedsphericalbessel.jl @@ -0,0 +1,40 @@ + +""" + sphericalbesselk(nu, x::T) where T <: {Float32, Float64} + +Computes `k_{ν}(x)`, the modified second-kind spherical Bessel function, and offers special branches for integer orders. +""" +sphericalbesselk(nu, x) = _sphericalbesselk(nu, float(x)) + +function _sphericalbesselk(nu, x::T) where T + isnan(x) && return NaN + if isinteger(nu) && nu < 41.5 + if x < zero(x) + return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported")) + end + # using ifelse here to hopefully cut out a branch on nu < 0 or not. The + # symmetry here is that + # k_{-n} = (...)*K_{-n + 1/2} + # = (...)*K_{|n| - 1/2} + # = (...)*K_{|n|-1 + 1/2} + # = k_{|n|-1} + _nu = ifelse(nu