From 5c1b7a69a128d2fac65b9a8d7949176e1c021969 Mon Sep 17 00:00:00 2001 From: Chris Geoga Date: Thu, 18 Aug 2022 09:14:09 -0400 Subject: [PATCH 1/9] Special half-integer branch for `besselk`. --- src/besselk.jl | 23 +++++++++++++++++++++++ src/constants.jl | 3 +++ test/besselk_test.jl | 6 ++++++ 3 files changed, 32 insertions(+) diff --git a/src/besselk.jl b/src/besselk.jl index 6ecdf42..e10760f 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -226,6 +226,9 @@ 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) + # check if nu is a half-integer: + isinteger(nu-1/2) && return besselk_vhalfint(nu, x) + # use uniform debye expansion if x or nu is large besselik_debye_cutoff(nu, x) && return besselk_large_orders(nu, x) @@ -424,3 +427,23 @@ 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 + + +""" + besselk_vhalfint(nu, x::T) where T <: {Float32, Float64} + +Computes `K_{ν}(x)` when `v + 1/2` is an integer using the fact that the +asymptotic expansion actually terminates and is exact for those specific `v` values. +""" +function besselk_vhalfint(v, x::T) where T + v = abs(v) + invx = inv(x) + b0 = b1 = SQRT_PID2(T)*sqrt(invx)*exp(-x) + twodx = 2*invx + _v = T(1/2) + while _v < v + b0, b1 = b1, muladd(b1, twodx*_v, b0) + _v += one(T) + end + b1 +end 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/test/besselk_test.jl b/test/besselk_test.jl index 8d82654..d305d15 100644 --- a/test/besselk_test.jl +++ b/test/besselk_test.jl @@ -67,6 +67,12 @@ k1x_32 = besselk1x.(Float32.(x)) @test besselk(100, 3.9) ≈ SpecialFunctions.besselk(100, 3.9) @test besselk(100, 234.0) ≈ SpecialFunctions.besselk(100, 234.0) +# test half-integer orders: +for v in (-3/2, -1/2, 1/2, 3/2, 5/2, 7/2, 9/2) + @test besselk(v, 7.0) ≈ SpecialFunctions.besselk(v, 7.0) + @test besselk(Float32(v), Float32(7.0)) ≈ SpecialFunctions.besselk(Float32(v), Float32(7.0)) +end + # test small arguments and order m = 0:40; x = [1e-6; 1e-4; 1e-3; 1e-2; 0.1; 1.0:2.0:500.0] for m in m, x in x From 5bd6e5931e8a96739d34261011e7de3f1471e346 Mon Sep 17 00:00:00 2001 From: Chris Geoga Date: Thu, 18 Aug 2022 11:08:45 -0400 Subject: [PATCH 2/9] Add `v` bound for half integer, and expand tests. --- src/besselk.jl | 4 +++- test/besselk_test.jl | 6 +++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/besselk.jl b/src/besselk.jl index e10760f..26da5d0 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -227,7 +227,7 @@ function besselk_positive_args(nu, x::T) where T <: Union{Float32, Float64} iszero(nu) && return besselk0(x) # check if nu is a half-integer: - isinteger(nu-1/2) && return besselk_vhalfint(nu, x) + besselk_vhalfint_check(nu, x) && return besselk_vhalfint(nu, x) # use uniform debye expansion if x or nu is large besselik_debye_cutoff(nu, x) && return besselk_large_orders(nu, x) @@ -447,3 +447,5 @@ function besselk_vhalfint(v, x::T) where T end b1 end +besselk_vhalfint_check(nu, x) = isinteger(nu-1/2) && (nu < 41.5) #@inline? + diff --git a/test/besselk_test.jl b/test/besselk_test.jl index d305d15..314b8df 100644 --- a/test/besselk_test.jl +++ b/test/besselk_test.jl @@ -68,9 +68,9 @@ k1x_32 = besselk1x.(Float32.(x)) @test besselk(100, 234.0) ≈ SpecialFunctions.besselk(100, 234.0) # test half-integer orders: -for v in (-3/2, -1/2, 1/2, 3/2, 5/2, 7/2, 9/2) - @test besselk(v, 7.0) ≈ SpecialFunctions.besselk(v, 7.0) - @test besselk(Float32(v), Float32(7.0)) ≈ SpecialFunctions.besselk(Float32(v), Float32(7.0)) +for (v,x) in Iterators.product(-35/2:1.0:81/2, range(0.0, 30.0, length=51)[2:end]) + @test besselk(v, x) ≈ SpecialFunctions.besselk(v, x) + @test besselk(Float32(v), Float32(x)) ≈ SpecialFunctions.besselk(Float32(v), Float32(x)) end # test small arguments and order From d9fc0e6b14d78a451f8584f19e297f32b2c743db Mon Sep 17 00:00:00 2001 From: Chris Geoga Date: Thu, 18 Aug 2022 14:20:52 -0400 Subject: [PATCH 3/9] Move implementation to modified spherical file, change branching a bit. --- src/Bessels.jl | 1 + src/besselk.jl | 28 +++++----------------------- src/modifiedsphericalbessel.jl | 33 +++++++++++++++++++++++++++++++++ 3 files changed, 39 insertions(+), 23 deletions(-) create mode 100644 src/modifiedsphericalbessel.jl 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 26da5d0..88048da 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -225,12 +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: - besselk_vhalfint_check(nu, x) && return besselk_vhalfint(nu, x) + (isinteger(nu-1/2) && !debye_cut) && return sphericalbesselk_int(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] @@ -428,24 +431,3 @@ 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 - -""" - besselk_vhalfint(nu, x::T) where T <: {Float32, Float64} - -Computes `K_{ν}(x)` when `v + 1/2` is an integer using the fact that the -asymptotic expansion actually terminates and is exact for those specific `v` values. -""" -function besselk_vhalfint(v, x::T) where T - v = abs(v) - invx = inv(x) - b0 = b1 = SQRT_PID2(T)*sqrt(invx)*exp(-x) - twodx = 2*invx - _v = T(1/2) - while _v < v - b0, b1 = b1, muladd(b1, twodx*_v, b0) - _v += one(T) - end - b1 -end -besselk_vhalfint_check(nu, x) = isinteger(nu-1/2) && (nu < 41.5) #@inline? - diff --git a/src/modifiedsphericalbessel.jl b/src/modifiedsphericalbessel.jl new file mode 100644 index 0000000..7910b2a --- /dev/null +++ b/src/modifiedsphericalbessel.jl @@ -0,0 +1,33 @@ + +function sphericalbesselk_int(v, x) + b0 = inv(x) + b1 = (x+one(x))/(x*x) + iszero(v) && return b0*exp(-x) + _v = one(v) + invx = inv(x) + while _v < v + _v += one(_v) + b0, b1 = b1, b0 + (2*_v - one(_v))*b1*invx + end + exp(-x)*b1 +end + +function _sphericalbesselk(nu, x::T) where T + if isinteger(nu) && nu < 41.5 + return sphericalbesselk_int(nu, x) + else + return inv(SQRT_PID2(T)*sqrt(x))*besselk(nu+1/2, x) + end +end + + +""" + 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)) + + + From 5ea2f3c4581b22efd501e9fcaddbf68775e9d393 Mon Sep 17 00:00:00 2001 From: Chris Geoga Date: Thu, 18 Aug 2022 14:28:04 -0400 Subject: [PATCH 4/9] Oof, forgot to remove _int subscript after prototyping. --- src/besselk.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/besselk.jl b/src/besselk.jl index 88048da..b1e2833 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -230,7 +230,7 @@ function besselk_positive_args(nu, x::T) where T <: Union{Float32, Float64} debye_cut = besselik_debye_cutoff(nu, x) # check if nu is a half-integer: - (isinteger(nu-1/2) && !debye_cut) && return sphericalbesselk_int(nu-1/2, x)*SQRT_PID2(T)*sqrt(x) + (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 debye_cut && return besselk_large_orders(nu, x) From 1273fcd09e52cc4862ba042b569d11ae2eedb551 Mon Sep 17 00:00:00 2001 From: Chris Geoga Date: Thu, 18 Aug 2022 14:56:35 -0400 Subject: [PATCH 5/9] Fixing docstring and adding ::Int annotation on the special method for the modified spherical bessel. --- src/modifiedsphericalbessel.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/modifiedsphericalbessel.jl b/src/modifiedsphericalbessel.jl index 7910b2a..e0a324b 100644 --- a/src/modifiedsphericalbessel.jl +++ b/src/modifiedsphericalbessel.jl @@ -1,5 +1,5 @@ -function sphericalbesselk_int(v, x) +function sphericalbesselk_int(v::Int, x) b0 = inv(x) b1 = (x+one(x))/(x*x) iszero(v) && return b0*exp(-x) @@ -24,8 +24,7 @@ end """ 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. +Computes `k_{ν}(x)`, the modified second-kind spherical Bessel function, and offers special branches for integer orders. """ sphericalbesselk(nu, x) = _sphericalbesselk(nu, float(x)) From d090b51b614884fb8a7636dbe349ebd77099ebe0 Mon Sep 17 00:00:00 2001 From: Chris Geoga Date: Thu, 18 Aug 2022 14:57:49 -0400 Subject: [PATCH 6/9] Re-order functions in file for readability. --- src/modifiedsphericalbessel.jl | 33 +++++++++++++++------------------ 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/src/modifiedsphericalbessel.jl b/src/modifiedsphericalbessel.jl index e0a324b..c31d983 100644 --- a/src/modifiedsphericalbessel.jl +++ b/src/modifiedsphericalbessel.jl @@ -1,4 +1,19 @@ +""" + 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 + if isinteger(nu) && nu < 41.5 + return sphericalbesselk_int(nu, x) + else + return inv(SQRT_PID2(T)*sqrt(x))*besselk(nu+1/2, x) + end +end + function sphericalbesselk_int(v::Int, x) b0 = inv(x) b1 = (x+one(x))/(x*x) @@ -12,21 +27,3 @@ function sphericalbesselk_int(v::Int, x) exp(-x)*b1 end -function _sphericalbesselk(nu, x::T) where T - if isinteger(nu) && nu < 41.5 - return sphericalbesselk_int(nu, x) - else - return inv(SQRT_PID2(T)*sqrt(x))*besselk(nu+1/2, x) - end -end - - -""" - 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)) - - - From 0671f11ab929305bf7144f23d643cec6a5317803 Mon Sep 17 00:00:00 2001 From: Chris Geoga Date: Thu, 18 Aug 2022 16:26:04 -0400 Subject: [PATCH 7/9] Fix negative order case for sphericalbesselk in special integer routine and add tests. --- src/modifiedsphericalbessel.jl | 9 ++++++++- test/sphericalbessel_test.jl | 6 ++++++ 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/src/modifiedsphericalbessel.jl b/src/modifiedsphericalbessel.jl index c31d983..0bd7953 100644 --- a/src/modifiedsphericalbessel.jl +++ b/src/modifiedsphericalbessel.jl @@ -8,7 +8,14 @@ sphericalbesselk(nu, x) = _sphericalbesselk(nu, float(x)) function _sphericalbesselk(nu, x::T) where T if isinteger(nu) && nu < 41.5 - return sphericalbesselk_int(nu, x) + # 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 Date: Thu, 18 Aug 2022 16:41:14 -0400 Subject: [PATCH 8/9] Faster sub for abs. --- src/modifiedsphericalbessel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/modifiedsphericalbessel.jl b/src/modifiedsphericalbessel.jl index 0bd7953..b9755e0 100644 --- a/src/modifiedsphericalbessel.jl +++ b/src/modifiedsphericalbessel.jl @@ -14,7 +14,7 @@ function _sphericalbesselk(nu, x::T) where T # = (...)*K_{|n| - 1/2} # = (...)*K_{|n|-1 + 1/2} # = k_{|n|-1} - _nu = ifelse(nu Date: Thu, 18 Aug 2022 17:13:36 -0400 Subject: [PATCH 9/9] Fix NaN propagation for x and check domain when x<0. --- src/modifiedsphericalbessel.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/modifiedsphericalbessel.jl b/src/modifiedsphericalbessel.jl index b9755e0..2a521af 100644 --- a/src/modifiedsphericalbessel.jl +++ b/src/modifiedsphericalbessel.jl @@ -7,7 +7,11 @@ Computes `k_{ν}(x)`, the modified second-kind spherical Bessel function, and of 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}