From 6d0a561f57af12983d31709d7973cff6526f2d83 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Fri, 19 Aug 2022 12:29:34 -0400 Subject: [PATCH 01/11] fix some spacing, overflow, and cutoff, add docs --- src/besselk.jl | 10 +++----- src/constants.jl | 3 --- src/modifiedsphericalbessel.jl | 46 +++++++++++++++++++++++----------- test/sphericalbessel_test.jl | 17 +++++++++++++ 4 files changed, 52 insertions(+), 24 deletions(-) diff --git a/src/besselk.jl b/src/besselk.jl index b1e2833..10fa654 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -226,14 +226,11 @@ 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) + # check if nu is a half-integer + (isinteger(nu-1/2) && sphericalbesselk_cutoff(nu)) && return sphericalbesselk_int(Int(nu-1/2), x)*SQPIO2(T)*sqrt(x) # use uniform debye expansion if x or nu is large - debye_cut && return besselk_large_orders(nu, x) + besselik_debye_cutoff(nu, x) && 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] @@ -430,4 +427,3 @@ 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 8568346..820e8d1 100644 --- a/src/constants.jl +++ b/src/constants.jl @@ -284,6 +284,3 @@ 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 index 2a521af..4f144cf 100644 --- a/src/modifiedsphericalbessel.jl +++ b/src/modifiedsphericalbessel.jl @@ -1,40 +1,58 @@ - +# Modified Spherical Bessel functions +# +# sphericalbesselk(nu, x) +# +# A numerical routine to compute the modified spherical bessel functions of the second kind. +# For moderate sized integer orders, forward recurrence is used starting from explicit formulas for k0(x) [1] and k1(x) [2]. +# Large orders are determined from the uniform asymptotic expansions (see src/besselk.jl for details) +# For non-integer orders, we directly call the besselk routine using the relation k_{n}(x) = sqrt(pi/(2x))*besselk(n+1/2, x) [3]. +# +# [1] http://dlmf.nist.gov/10.49.E12 +# [2] http://dlmf.nist.gov/10.49.E13 +# [3] http://dlmf.nist.gov/10.47.E9 +# """ 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)) +sphericalbesselk(nu::Real, x::Real) = _sphericalbesselk(nu, float(x)) + +_sphericalbesselk(nu, x::Float16) = Float16(_sphericalbesselk(nu, Float32(x))) -function _sphericalbesselk(nu, x::T) where T - isnan(x) && return NaN - if isinteger(nu) && nu < 41.5 +function _sphericalbesselk(nu, x::T) where T <: Union{Float32, Float64} + if ~isfinite(x) + isnan(x) && return x + isinf(x) && return zero(x) + end + if isinteger(nu) && sphericalbesselk_cutoff(nu) 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 + # using ifelse here to 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} + # = k_{|n|-1} _nu = ifelse(nu Date: Fri, 19 Aug 2022 13:53:13 -0400 Subject: [PATCH 02/11] update news --- NEWS.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index cb83b31..fdbd68a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -11,7 +11,8 @@ For bug fixes, performance enhancements, or fixes to unexported functions we wil # Unreleased ### Added - - Add more optimized methods for Float32 calculations that are faster([PR #43](https://github.com/JuliaMath/Bessels.jl/pull/43)) + - Add more optimized methods for Float32 calculations that are faster ([PR #43](https://github.com/JuliaMath/Bessels.jl/pull/43)) + - Add methods for computing modified spherical bessel function of second ([PR #46](https://github.com/JuliaMath/Bessels.jl/pull/46) and ([PR #47](https://github.com/JuliaMath/Bessels.jl/pull/47))) currently unexported closes ([Issue #25](https://github.com/JuliaMath/Bessels.jl/issues/25)) ### Fixed - Reduce compile time and time to first call of besselj and bessely ([PR #42](https://github.com/JuliaMath/Bessels.jl/pull/42)) From 43b5e200ac78b64a09041302110c678a4d055b92 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Fri, 19 Aug 2022 14:03:49 -0400 Subject: [PATCH 03/11] update reference --- src/modifiedsphericalbessel.jl | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/modifiedsphericalbessel.jl b/src/modifiedsphericalbessel.jl index 4f144cf..4723d65 100644 --- a/src/modifiedsphericalbessel.jl +++ b/src/modifiedsphericalbessel.jl @@ -3,13 +3,11 @@ # sphericalbesselk(nu, x) # # A numerical routine to compute the modified spherical bessel functions of the second kind. -# For moderate sized integer orders, forward recurrence is used starting from explicit formulas for k0(x) [1] and k1(x) [2]. +# For moderate sized integer orders, forward recurrence is used starting from explicit formulas for k0(x) = exp(-x) / x and k1(x) = k0(x) * (x+1) / x [1]. # Large orders are determined from the uniform asymptotic expansions (see src/besselk.jl for details) -# For non-integer orders, we directly call the besselk routine using the relation k_{n}(x) = sqrt(pi/(2x))*besselk(n+1/2, x) [3]. +# For non-integer orders, we directly call the besselk routine using the relation k_{n}(x) = sqrt(pi/(2x))*besselk(n+1/2, x) [1]. # -# [1] http://dlmf.nist.gov/10.49.E12 -# [2] http://dlmf.nist.gov/10.49.E13 -# [3] http://dlmf.nist.gov/10.47.E9 +# [1] https://mathworld.wolfram.com/ModifiedSphericalBesselFunctionoftheSecondKind.html # """ sphericalbesselk(nu, x::T) where T <: {Float32, Float64} From a3c767fb32ba63a72f8ddda9f79ea99dd4f02e9d Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Fri, 19 Aug 2022 16:05:10 -0400 Subject: [PATCH 04/11] add sphericalbesseli --- src/besseli.jl | 9 +++++- src/modifiedsphericalbessel.jl | 54 ++++++++++++++++++++++++++++++---- test/sphericalbessel_test.jl | 39 +++++++++++++++++++++++- 3 files changed, 94 insertions(+), 8 deletions(-) diff --git a/src/besseli.jl b/src/besseli.jl index 6dc87ba..2f146b5 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -170,6 +170,10 @@ _besseli(nu, x::Float16) = Float16(_besseli(nu, Float32(x))) function _besseli(nu, x::T) where T <: Union{Float32, Float64} isinteger(nu) && return _besseli(Int(nu), x) + if ~isfinite(x) + isnan(x) && return x + isinf(x) && return x + end abs_nu = abs(nu) abs_x = abs(x) @@ -192,6 +196,10 @@ function _besseli(nu, x::T) where T <: Union{Float32, Float64} end end function _besseli(nu::Integer, x::T) where T <: Union{Float32, Float64} + if ~isfinite(x) + isnan(x) && return x + isinf(x) && return x + end abs_nu = abs(nu) abs_x = abs(x) sg = iseven(abs_nu) ? 1 : -1 @@ -211,7 +219,6 @@ Modified Bessel function of the first kind of order nu, ``I_{nu}(x)`` for positi function besseli_positive_args(nu, x::T) where T <: Union{Float32, Float64} iszero(nu) && return besseli0(x) isone(nu) && return besseli1(x) - isinf(x) && return T(Inf) # use large argument expansion if x >> nu besseli_large_argument_cutoff(nu, x) && return besseli_large_argument(nu, x) diff --git a/src/modifiedsphericalbessel.jl b/src/modifiedsphericalbessel.jl index 4723d65..86dd17e 100644 --- a/src/modifiedsphericalbessel.jl +++ b/src/modifiedsphericalbessel.jl @@ -1,13 +1,16 @@ # Modified Spherical Bessel functions # -# sphericalbesselk(nu, x) +# sphericalbesseli(nu, x), sphericalbesselk(nu, x) # -# A numerical routine to compute the modified spherical bessel functions of the second kind. -# For moderate sized integer orders, forward recurrence is used starting from explicit formulas for k0(x) = exp(-x) / x and k1(x) = k0(x) * (x+1) / x [1]. +# A numerical routine to compute the modified spherical bessel functions of the first and second kind. +# The modified spherical bessel function of the first kind is computed using the power series for small arguments, +# explicit formulas for (nu=0,1,2), and using its relation to besseli for other arguments [1]. +# The modified bessel function of the second kind is computed for small to moderate integer orders using forward recurrence starting from explicit formulas for k0(x) = exp(-x) / x and k1(x) = k0(x) * (x+1) / x [2]. # Large orders are determined from the uniform asymptotic expansions (see src/besselk.jl for details) -# For non-integer orders, we directly call the besselk routine using the relation k_{n}(x) = sqrt(pi/(2x))*besselk(n+1/2, x) [1]. -# -# [1] https://mathworld.wolfram.com/ModifiedSphericalBesselFunctionoftheSecondKind.html +# For non-integer orders, we directly call the besselk routine using the relation k_{n}(x) = sqrt(pi/(2x))*besselk(n+1/2, x) [2]. +# +# [1] https://mathworld.wolfram.com/ModifiedBesselFunctionoftheFirstKind.html +# [2] https://mathworld.wolfram.com/ModifiedSphericalBesselFunctionoftheSecondKind.html # """ sphericalbesselk(nu, x::T) where T <: {Float32, Float64} @@ -54,3 +57,42 @@ function sphericalbesselk_int(v::Int, x) end b1 end + +""" + sphericalbesseli(nu, x::T) where T <: {Float32, Float64} + +Computes `i_{ν}(x)`, the modified first-kind spherical Bessel function. +""" +sphericalbesseli(nu::Real, x::Real) = _sphericalbesseli(nu, float(x)) + +_sphericalbesseli(nu, x::Float16) = Float16(_sphericalbesseli(nu, Float32(x))) + +function _sphericalbesseli(nu, x::T) where T <: Union{Float32, Float64} + isinf(x) && return x + x < zero(x) && throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported")) + + sphericalbesselj_small_args_cutoff(nu, x::T) && return sphericalbesseli_small_args(nu, x) + isinteger(nu) && return _sphericalbesseli_small_orders(Int(nu), x) + return SQPIO2(T)*besseli(nu+1/2, x) / sqrt(x) +end + +function _sphericalbesseli_small_orders(nu::Integer, x::T) where T + # prone to cancellation in the subtraction + # best to expand and group + nu_abs = abs(nu) + x2 = x*x + sinhx = sinh(x) + coshx = cosh(x) + nu_abs == 0 && return sinhx / x + nu_abs == 1 && return (x*coshx - sinhx) / x2 + nu_abs == 2 && return (x2*sinhx + 3*(sinhx - x*coshx)) / (x2*x) + return SQPIO2(T)*besseli(nu+1/2, x) / sqrt(x) +end + +function sphericalbesseli_small_args(nu, x::T) where T + iszero(x) && return iszero(nu) ? one(T) : x + x2 = x^2 / 4 + coef = evalpoly(x2, (1, inv(T(3)/2 + nu), inv(5 + nu), inv(T(21)/2 + nu), inv(18 + nu))) + a = SQPIO2(T) / (gamma(T(3)/2 + nu) * 2^(nu + T(1)/2)) + return x^nu * a * coef +end diff --git a/test/sphericalbessel_test.jl b/test/sphericalbessel_test.jl index 6188ca4..342afcc 100644 --- a/test/sphericalbessel_test.jl +++ b/test/sphericalbessel_test.jl @@ -11,6 +11,14 @@ x = 1e-15 @test Bessels.sphericalbesselk(5.5, x) ≈ SpecialFunctions.besselk(5.5 + 1/2, x) * sqrt( 2 / (x*pi)) @test Bessels.sphericalbesselk(10, x) ≈ SpecialFunctions.besselk(10 + 1/2, x) * sqrt( 2 / (x*pi)) +x = 1e-20 +@test Bessels.sphericalbesseli(0, x) ≈ SpecialFunctions.besseli(0 + 1/2, x) * sqrt( pi / (x*2)) +@test Bessels.sphericalbesseli(1, x) ≈ SpecialFunctions.besseli(1 + 1/2, x) * sqrt( pi / (x*2)) +@test Bessels.sphericalbesseli(2, x) ≈ SpecialFunctions.besseli(2 + 1/2, x) * sqrt( pi / (x*2)) +@test Bessels.sphericalbesseli(3, x) ≈ SpecialFunctions.besseli(3 + 1/2, x) * sqrt( pi / (x*2)) +@test Bessels.sphericalbesseli(4, x) ≈ SpecialFunctions.besseli(4 + 1/2, x) * sqrt( pi / (x*2)) +@test Bessels.sphericalbesseli(6.5, x) ≈ SpecialFunctions.besseli(6.5 + 1/2, x) * sqrt( pi / (x*2)) + # test zero @test isone(Bessels.sphericalbesselj(0, 0.0)) @test iszero(Bessels.sphericalbesselj(3, 0.0)) @@ -26,6 +34,14 @@ x = 1e-15 @test isinf(Bessels.sphericalbesselk(4, 0.0)) @test isinf(Bessels.sphericalbesselk(10.2, 0.0)) +x = 0.0 +@test isone(Bessels.sphericalbesseli(0, x)) +@test iszero(Bessels.sphericalbesseli(1, x)) +@test iszero(Bessels.sphericalbesseli(2, x)) +@test iszero(Bessels.sphericalbesseli(3, x)) +@test iszero(Bessels.sphericalbesseli(4, x)) +@test iszero(Bessels.sphericalbesseli(6.4, x)) + # test Inf @test iszero(Bessels.sphericalbesselj(1, Inf)) @test iszero(Bessels.sphericalbesselj(10.2, Inf)) @@ -36,6 +52,14 @@ x = 1e-15 @test iszero(Bessels.sphericalbesselk(4, Inf)) @test iszero(Bessels.sphericalbesselk(10.2, Inf)) +x = Inf +@test isinf(Bessels.sphericalbesseli(0, x)) +@test isinf(Bessels.sphericalbesseli(1, x)) +@test isinf(Bessels.sphericalbesseli(2, x)) +@test isinf(Bessels.sphericalbesseli(3, x)) +@test isinf(Bessels.sphericalbesseli(4, x)) +@test isinf(Bessels.sphericalbesseli(6.4, x)) + # test NaN @test isnan(Bessels.sphericalbesselj(1.4, NaN)) @test isnan(Bessels.sphericalbesselj(4.0, NaN)) @@ -45,6 +69,14 @@ x = 1e-15 @test isnan(Bessels.sphericalbesselk(1.4, NaN)) @test isnan(Bessels.sphericalbesselk(4.0, NaN)) +x = NaN +@test isnan(Bessels.sphericalbesseli(0, x)) +@test isnan(Bessels.sphericalbesseli(1, x)) +@test isnan(Bessels.sphericalbesseli(2, x)) +@test isnan(Bessels.sphericalbesseli(3, x)) +@test isnan(Bessels.sphericalbesseli(4, x)) +@test isnan(Bessels.sphericalbesseli(6.4, x)) + # test Float16, Float32 types @test Bessels.sphericalbesselj(Float16(1.4), Float16(1.2)) isa Float16 @test Bessels.sphericalbessely(Float16(1.4), Float16(1.2)) isa Float16 @@ -54,16 +86,21 @@ x = 1e-15 @test Bessels.sphericalbesselk(Float16(1.4), Float16(1.2)) isa Float16 @test Bessels.sphericalbesselk(1.0f0, 1.2f0) isa Float32 -for x in 0.5:1.0:100.0, v in [0, 1, 5.5, 8.2, 10] +@test Bessels.sphericalbesseli(Float16(1.4), Float16(1.2)) isa Float16 +@test Bessels.sphericalbesseli(1.0f0, 1.2f0) isa Float32 + +for x in 0.5:1.5:100.0, v in [0, 1, 2, 3, 4, 5.5, 8.2, 10] @test isapprox(Bessels.sphericalbesselj(v, x), SpecialFunctions.sphericalbesselj(v, x), rtol=1e-12) @test isapprox(Bessels.sphericalbessely(v, x), SpecialFunctions.sphericalbessely(v, x), rtol=1e-12) @test isapprox(Bessels.sphericalbesselk(v, x), SpecialFunctions.besselk(v+1/2, x) * sqrt( 2 / (x*pi)), rtol=1e-12) + @test isapprox(Bessels.sphericalbesseli(v, x), SpecialFunctions.besseli(v+1/2, x) * sqrt( pi / (x*2)), rtol=1e-12) end for x in 5.5:4.0:160.0, v in [20, 25.0, 32.4, 40.0, 45.12, 50.0, 55.2, 60.124, 70.23, 75.0, 80.0, 92.3, 100.0, 120.0] @test isapprox(Bessels.sphericalbesselj(v, x), SpecialFunctions.sphericalbesselj(v, x), rtol=3e-12) @test isapprox(Bessels.sphericalbessely(v, x), SpecialFunctions.sphericalbessely(v, x), rtol=3e-12) @test isapprox(Bessels.sphericalbesselk(v, x), SpecialFunctions.besselk(v+1/2, x) * sqrt( 2 / (x*pi)), rtol=1e-12) + @test isapprox(Bessels.sphericalbesseli(v, x), SpecialFunctions.besseli(v+1/2, x) * sqrt( pi / (x*2)), rtol=1e-12) end @test isapprox(Bessels.sphericalbessely(270, 240.0), SpecialFunctions.sphericalbessely(270, 240.0), rtol=3e-12) From 4f41688389420e553ea3ec5cb215ca43f1e44b40 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Fri, 19 Aug 2022 16:21:38 -0400 Subject: [PATCH 05/11] return x ~isfinite --- src/besseli.jl | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/src/besseli.jl b/src/besseli.jl index 2f146b5..ce9c967 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -170,10 +170,7 @@ _besseli(nu, x::Float16) = Float16(_besseli(nu, Float32(x))) function _besseli(nu, x::T) where T <: Union{Float32, Float64} isinteger(nu) && return _besseli(Int(nu), x) - if ~isfinite(x) - isnan(x) && return x - isinf(x) && return x - end + ~isfinite(x) && return x abs_nu = abs(nu) abs_x = abs(x) @@ -196,10 +193,7 @@ function _besseli(nu, x::T) where T <: Union{Float32, Float64} end end function _besseli(nu::Integer, x::T) where T <: Union{Float32, Float64} - if ~isfinite(x) - isnan(x) && return x - isinf(x) && return x - end + ~isfinite(x) && return x abs_nu = abs(nu) abs_x = abs(x) sg = iseven(abs_nu) ? 1 : -1 From 2f0c45eed7dffceb0c177e3740c72c7bbc067f2b Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 19 Sep 2022 12:40:40 -0400 Subject: [PATCH 06/11] resolve some conflicts --- src/besselk.jl | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/src/besselk.jl b/src/besselk.jl index 10fa654..c90d669 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -229,6 +229,9 @@ function besselk_positive_args(nu, x::T) where T <: Union{Float32, Float64} # check if nu is a half-integer (isinteger(nu-1/2) && sphericalbesselk_cutoff(nu)) && return sphericalbesselk_int(Int(nu-1/2), x)*SQPIO2(T)*sqrt(x) + # check if the standard asymptotic expansion can be used + besseli_large_argument_cutoff(nu, x) && return besselk_large_argument(nu, x) + # use uniform debye expansion if x or nu is large besselik_debye_cutoff(nu, x) && return besselk_large_orders(nu, x) @@ -427,3 +430,37 @@ 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 + +##### +##### Large argument expansion for K_{nu}(x) +##### + +# Computes the asymptotic expansion of K_ν w.r.t. argument. +# Accurate for large x, and faster than uniform asymptotic expansion for small to small-ish orders +# See http://dlmf.nist.gov/10.40.E2 + +function besselk_large_argument(v, x::T) where T + a = exp(-x / 2) + coef = a * sqrt(pi / 2x) + return T(_besselk_large_argument(v, x) * coef * a) +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 + MaxIter = 5000 + S = promote_type(T, Float64) + v, x = S(v), S(x) + + fv2 = 4 * v^2 + term = one(S) + res = term + s = term + for i in 1:MaxIter + offset = muladd(2, i, -1) + term *= muladd(offset, -offset, fv2) / (8 * x * i) + res = muladd(term, s, res) + abs(term) <= eps(T) && break + end + return res +end From 3ccf5df8feb853928c097d162be02b9855b8853e Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 19 Sep 2022 12:51:51 -0400 Subject: [PATCH 07/11] fix merge --- src/besselk.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/besselk.jl b/src/besselk.jl index 04d990d..ba8e524 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -232,8 +232,8 @@ function besselk_positive_args(nu, x::T) where T <: Union{Float32, Float64} # check if the standard asymptotic expansion can be used besseli_large_argument_cutoff(nu, x) && return besselk_large_argument(nu, x) - # check if the standard asymptotic expansion can be used: - besselk_asexp_cutoff(nu, x) && return besselk_large_argument(nu, x) + # check if the standard asymptotic expansion can be used + besseli_large_argument_cutoff(nu, x) && return besselk_large_argument(nu, x) # use uniform debye expansion if x or nu is large besselik_debye_cutoff(nu, x) && return besselk_large_orders(nu, x) From 08b1c53f230b972755fa9a82f0e4ecb7373930ed Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 19 Sep 2022 13:00:24 -0400 Subject: [PATCH 08/11] update docs and add news --- NEWS.md | 3 ++- src/besselk.jl | 2 ++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index fdbd68a..564f3ef 100644 --- a/NEWS.md +++ b/NEWS.md @@ -12,7 +12,8 @@ For bug fixes, performance enhancements, or fixes to unexported functions we wil ### Added - Add more optimized methods for Float32 calculations that are faster ([PR #43](https://github.com/JuliaMath/Bessels.jl/pull/43)) - - Add methods for computing modified spherical bessel function of second ([PR #46](https://github.com/JuliaMath/Bessels.jl/pull/46) and ([PR #47](https://github.com/JuliaMath/Bessels.jl/pull/47))) currently unexported closes ([Issue #25](https://github.com/JuliaMath/Bessels.jl/issues/25)) + - Add methods for computing modified spherical bessel function of second ([PR #46](https://github.com/JuliaMath/Bessels.jl/pull/46) contributed by @cgeoga and first ([PR #47](https://github.com/JuliaMath/Bessels.jl/pull/47))) kind. These functions are currently not exported. Closes ([Issue #25](https://github.com/JuliaMath/Bessels.jl/issues/25)) + - Asymptotic expansion for x >> nu was added ([PR #48](https://github.com/JuliaMath/Bessels.jl/pull/48)) that decreases computation time for large arguments. Contributed by @cgeoga ### Fixed - Reduce compile time and time to first call of besselj and bessely ([PR #42](https://github.com/JuliaMath/Bessels.jl/pull/42)) diff --git a/src/besselk.jl b/src/besselk.jl index ba8e524..57e2c80 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -160,6 +160,7 @@ end # # The identities are computed by calling the `besseli_positive_args(nu, x)` function which computes K_{ν}(x) # for positive arguments and orders. For large orders, Debye's uniform asymptotic expansions are used. +# For large arguments x >> nu, large argument expansion is used [9]. # For small value and when nu > ~x the power series is used. The rest of the values are computed using slightly different methods. # The power series for besseli is modified to give both I_{v} and I_{v-1} where the ratio K_{v+1} / K_{v} is computed using continued fractions [8]. # The wronskian connection formula is then used to compute K_v. @@ -178,6 +179,7 @@ end # arXiv preprint arXiv:2201.00090 (2022). # [8] Cuyt, A. A., Petersen, V., Verdonk, B., Waadeland, H., & Jones, W. B. (2008). # Handbook of continued fractions for special functions. Springer Science & Business Media. +# [9] http://dlmf.nist.gov/10.40.E2 # """ From f38aef6461f4779d3912b89e8229590863b679cd Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 19 Sep 2022 13:01:31 -0400 Subject: [PATCH 09/11] fix repeat line --- src/besselk.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/besselk.jl b/src/besselk.jl index 57e2c80..27e666a 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -234,9 +234,6 @@ function besselk_positive_args(nu, x::T) where T <: Union{Float32, Float64} # check if the standard asymptotic expansion can be used besseli_large_argument_cutoff(nu, x) && return besselk_large_argument(nu, x) - # check if the standard asymptotic expansion can be used - besseli_large_argument_cutoff(nu, x) && return besselk_large_argument(nu, x) - # use uniform debye expansion if x or nu is large besselik_debye_cutoff(nu, x) && return besselk_large_orders(nu, x) From 20de7e2499d47c89bae331aa0c00ba6f04104676 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 19 Sep 2022 13:05:41 -0400 Subject: [PATCH 10/11] use expansion in scaled --- src/besselk.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/besselk.jl b/src/besselk.jl index 27e666a..0bfbe6b 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -259,6 +259,9 @@ function _besselkx(nu, x::T) where T <: Union{Float32, Float64} # dispatch to avoid uniform expansion when nu = 0 iszero(nu) && return besselk0x(x) + # check if the standard asymptotic expansion can be used + besseli_large_argument_cutoff(nu, x) && return besselk_large_argument_scaled(nu, x) + # use uniform debye expansion if x or nu is large besselik_debye_cutoff(nu, x) && return besselk_large_orders_scaled(nu, x) From 26ae05c1134ccd42b2191dddcde44ada5575d232 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 19 Sep 2022 13:11:53 -0400 Subject: [PATCH 11/11] bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 3709c06..bf23bbf 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.0" +version = "0.2.1" [compat] julia = "1.6"