From f1ff16770c6f085354345efaebafa6e94a8246a4 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 3 Oct 2022 18:26:20 -0400 Subject: [PATCH 01/11] add array support --- src/besseli.jl | 11 +++++++-- src/besselj.jl | 16 ++++++++++++-- src/besselk.jl | 11 +++++++-- src/bessely.jl | 11 +++++++-- src/recurrence.jl | 43 +++++++++++++++++++++++++++++++++--- test/besseli_test.jl | 8 ++++++- test/besselj_test.jl | 6 +++++ test/besselk_test.jl | 4 ++++ test/bessely_test.jl | 6 +++++ test/sphericalbessel_test.jl | 4 ++-- 10 files changed, 106 insertions(+), 14 deletions(-) diff --git a/src/besseli.jl b/src/besseli.jl index ce9c967..7d9e520 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -164,11 +164,11 @@ end Modified Bessel function of the second kind of order nu, ``I_{nu}(x)``. """ -besseli(nu::Real, x::Real) = _besseli(nu, float(x)) +besseli(nu, x::Real) = _besseli(nu, float(x)) _besseli(nu, x::Float16) = Float16(_besseli(nu, Float32(x))) -function _besseli(nu, x::T) where T <: Union{Float32, Float64} +function _besseli(nu::T, x::T) where T <: Union{Float32, Float64} isinteger(nu) && return _besseli(Int(nu), x) ~isfinite(x) && return x abs_nu = abs(nu) @@ -205,6 +205,13 @@ function _besseli(nu::Integer, x::T) where T <: Union{Float32, Float64} end end +function _besseli(nu::AbstractRange, x::T) where T + (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1")) + out = Vector{T}(undef, length(nu)) + out[end-1], out[end] = _besseli(nu[end-1], x), _besseli(nu[end], x) + return besselk_down_recurrence!(out, x, nu) +end + """ besseli_positive_args(nu, x::T) where T <: Union{Float32, Float64} diff --git a/src/besselj.jl b/src/besselj.jl index 4db6838..69f09b5 100644 --- a/src/besselj.jl +++ b/src/besselj.jl @@ -206,11 +206,11 @@ end Bessel function of the first kind of order nu, ``J_{nu}(x)``. nu and x must be real where nu and x can be positive or negative. """ -besselj(nu::Real, x::Real) = _besselj(nu, float(x)) +besselj(nu, x::Real) = _besselj(nu, float(x)) _besselj(nu, x::Float16) = Float16(_besselj(nu, Float32(x))) -function _besselj(nu, x::T) where T <: Union{Float32, Float64} +function _besselj(nu::T, x::T) where T <: Union{Float32, Float64} isinteger(nu) && return _besselj(Int(nu), x) abs_nu = abs(nu) abs_x = abs(x) @@ -255,6 +255,18 @@ function _besselj(nu::Integer, x::T) where T <: Union{Float32, Float64} end end +function _besselj(nu::AbstractRange, x::T) where T + (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1")) + out = Vector{T}(undef, length(nu)) + if nu[end] < x + out[1], out[2] = _besselj(nu[1], x), _besselj(nu[2], x) + return besselj_up_recurrence!(out, x, nu) + else + out[end-1], out[end] = _besselj(nu[end-1], x), _besselj(nu[end], x) + return besselj_down_recurrence!(out, x, nu) + end +end + """ besselj_positive_args(nu, x::T) where T <: Float64 diff --git a/src/besselk.jl b/src/besselk.jl index 0bfbe6b..db89517 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -187,11 +187,11 @@ end Modified Bessel function of the second kind of order nu, ``K_{nu}(x)``. """ -besselk(nu::Real, x::Real) = _besselk(nu, float(x)) +besselk(nu, x::Real) = _besselk(nu, float(x)) _besselk(nu, x::Float16) = Float16(_besselk(nu, Float32(x))) -function _besselk(nu, x::T) where T <: Union{Float32, Float64} +function _besselk(nu::T, x::T) where T <: Union{Float32, Float64} isinteger(nu) && return _besselk(Int(nu), x) abs_nu = abs(nu) abs_x = abs(x) @@ -216,6 +216,13 @@ function _besselk(nu::Integer, x::T) where T <: Union{Float32, Float64} end end +function _besselk(nu::AbstractRange, x::T) where T + (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1")) + out = Vector{T}(undef, length(nu)) + out[1], out[2] = _besselk(nu[1], x), _besselk(nu[2], x) + return besselk_up_recurrence!(out, x, nu) +end + """ besselk_positive_args(x::T) where T <: Union{Float32, Float64} diff --git a/src/bessely.jl b/src/bessely.jl index 45987eb..a460b7b 100644 --- a/src/bessely.jl +++ b/src/bessely.jl @@ -243,11 +243,11 @@ end Bessel function of the second kind of order nu, ``Y_{nu}(x)``. nu and x must be real where nu and x can be positive or negative. """ -bessely(nu::Real, x::Real) = _bessely(nu, float(x)) +bessely(nu, x::Real) = _bessely(nu, float(x)) _bessely(nu, x::Float16) = Float16(_bessely(nu, Float32(x))) -function _bessely(nu, x::T) where T <: Union{Float32, Float64} +function _bessely(nu::T, x::T) where T <: Union{Float32, Float64} isnan(nu) || isnan(x) && return NaN isinteger(nu) && return _bessely(Int(nu), x) abs_nu = abs(nu) @@ -296,6 +296,13 @@ function _bessely(nu::Integer, x::T) where T <: Union{Float32, Float64} end end +function _bessely(nu::AbstractRange, x::T) where T + (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1")) + out = Vector{T}(undef, length(nu)) + out[1], out[2] = _bessely(nu[1], x), _bessely(nu[2], x) + return besselj_up_recurrence!(out, x, nu) +end + ##### ##### `bessely` for positive arguments and orders ##### diff --git a/src/recurrence.jl b/src/recurrence.jl index c959537..cca9ff1 100644 --- a/src/recurrence.jl +++ b/src/recurrence.jl @@ -2,7 +2,7 @@ # outputs both (bessel(x, nu_end), bessel(x, nu_end+1) # x = 0.1; k0 = besselk(0,x); k1 = besselk(1,x); # besselk_up_recurrence(x, k1, k0, 1, 5) will give besselk(5, x) -@inline function besselk_up_recurrence(x, jnu, jnum1, nu_start, nu_end) +function besselk_up_recurrence(x, jnu, jnum1, nu_start, nu_end) x2 = 2 / x while nu_start < nu_end + 0.5 # avoid inexact floating points when nu isa float jnum1, jnu = jnu, muladd(nu_start*x2, jnu, jnum1) @@ -10,12 +10,21 @@ end return jnum1, jnu end +function besselk_up_recurrence!(out, x::T, nu_range) where T + x2 = 2 / x + k = 3 + for nu in nu_range[2:end-1] + out[k] = muladd(nu*x2, out[k-1], out[k-2]) + k += 1 + end + return out +end # forward recurrence relation for besselj and bessely # outputs both (bessel(x, nu_end), bessel(x, nu_end+1) # x = 0.1; y0 = bessely0(x); y1 = bessely1(x); # besselj_up_recurrence(x, y1, y0, 1, 5) will give bessely(5, x) -@inline function besselj_up_recurrence(x::T, jnu, jnum1, nu_start, nu_end) where T +function besselj_up_recurrence(x::T, jnu, jnum1, nu_start, nu_end) where T x2 = 2 / x while nu_start < nu_end + 0.5 # avoid inexact floating points when nu isa float jnum1, jnu = jnu, muladd(nu_start*x2, jnu, -jnum1) @@ -27,12 +36,21 @@ end # NaN inputs need to be checked before getting to this routine return isnan(jnum1) ? (-T(Inf), -T(Inf)) : (jnum1, jnu) end +function besselj_up_recurrence!(out, x::T, nu_range) where T + x2 = 2 / x + k = 3 + for nu in nu_range[2:end-1] + out[k] = muladd(nu*x2, out[k-1], -out[k-2]) + k += 1 + end + return out +end # backward recurrence relation for besselj and bessely # outputs both (bessel(x, nu_end), bessel(x, nu_end-1) # x = 0.1; j10 = besselj(10, x); j11 = besselj(11, x); # besselj_down_recurrence(x, j10, j11, 10, 1) will give besselj(1, x) -@inline function besselj_down_recurrence(x, jnu, jnup1, nu_start, nu_end) +function besselj_down_recurrence(x, jnu, jnup1, nu_start, nu_end) x2 = 2 / x while nu_start > nu_end - 0.5 jnup1, jnu = jnu, muladd(nu_start*x2, jnu, -jnup1) @@ -41,6 +59,16 @@ end return jnup1, jnu end +function besselj_down_recurrence!(out, x::T, nu_range) where T + x2 = 2 / x + k = length(nu_range) - 2 + for nu in nu_range[end-1:-1:2] + out[k] = muladd(nu*x2, out[k+1], -out[k+2]) + k -= 1 + end + return out +end + #= # currently not used # backward recurrence relation for besselk and besseli @@ -56,3 +84,12 @@ end return jnup1, jnu end =# +function besselk_down_recurrence!(out, x::T, nu_range) where T + x2 = 2 / x + k = length(nu_range) - 2 + for nu in nu_range[end-1:-1:2] + out[k] = muladd(nu*x2, out[k+1], out[k+2]) + k -= 1 + end + return out +end diff --git a/test/besseli_test.jl b/test/besseli_test.jl index 9f15400..b05ae79 100644 --- a/test/besseli_test.jl +++ b/test/besseli_test.jl @@ -94,7 +94,7 @@ t = [besselix(m, x) for m in m, x in x] @test t ≈ [SpecialFunctions.besselix(m, x) for m in m, x in x] @test besselix(10, Float16(1.0)) isa Float16 -## Tests for besselk +## Tests for besseli ## test all numbers and orders for 0 Date: Tue, 4 Oct 2022 10:16:39 -0400 Subject: [PATCH 02/11] fix float16 for nu arguments --- src/besseli.jl | 5 +++-- src/besselj.jl | 2 +- src/besselk.jl | 2 +- src/bessely.jl | 2 +- src/modifiedsphericalbessel.jl | 10 +++++----- src/sphericalbessel.jl | 6 +++--- test/besseli_test.jl | 4 ++-- test/besselj_test.jl | 2 +- test/besselk_test.jl | 4 ++-- test/bessely_test.jl | 2 +- test/sphericalbessel_test.jl | 4 ++-- 11 files changed, 22 insertions(+), 21 deletions(-) diff --git a/src/besseli.jl b/src/besseli.jl index 7d9e520..fed01c5 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -164,9 +164,10 @@ end Modified Bessel function of the second kind of order nu, ``I_{nu}(x)``. """ +# perhaps have two besseli(nu::Real, x::Real) and besseli(nu::AbstractRange, x::Real) besseli(nu, x::Real) = _besseli(nu, float(x)) -_besseli(nu, x::Float16) = Float16(_besseli(nu, Float32(x))) +_besseli(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_besseli(Float32(nu), Float32(x))) function _besseli(nu::T, x::T) where T <: Union{Float32, Float64} isinteger(nu) && return _besseli(Int(nu), x) @@ -239,7 +240,7 @@ Nu must be real. """ besselix(nu::Real, x::Real) = _besselix(nu, float(x)) -_besselix(nu, x::Float16) = Float16(_besselix(nu, Float32(x))) +_besselix(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_besselix(Float32(nu), Float32(x))) function _besselix(nu, x::T) where T <: Union{Float32, Float64} iszero(nu) && return besseli0x(x) diff --git a/src/besselj.jl b/src/besselj.jl index 69f09b5..34e920e 100644 --- a/src/besselj.jl +++ b/src/besselj.jl @@ -208,7 +208,7 @@ nu and x must be real where nu and x can be positive or negative. """ besselj(nu, x::Real) = _besselj(nu, float(x)) -_besselj(nu, x::Float16) = Float16(_besselj(nu, Float32(x))) +_besselj(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_besselj(Float32(nu), Float32(x))) function _besselj(nu::T, x::T) where T <: Union{Float32, Float64} isinteger(nu) && return _besselj(Int(nu), x) diff --git a/src/besselk.jl b/src/besselk.jl index db89517..fa60241 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -189,7 +189,7 @@ Modified Bessel function of the second kind of order nu, ``K_{nu}(x)``. """ besselk(nu, x::Real) = _besselk(nu, float(x)) -_besselk(nu, x::Float16) = Float16(_besselk(nu, Float32(x))) +_besselk(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_besselk(Float32(nu), Float32(x))) function _besselk(nu::T, x::T) where T <: Union{Float32, Float64} isinteger(nu) && return _besselk(Int(nu), x) diff --git a/src/bessely.jl b/src/bessely.jl index a460b7b..c392b64 100644 --- a/src/bessely.jl +++ b/src/bessely.jl @@ -245,7 +245,7 @@ nu and x must be real where nu and x can be positive or negative. """ bessely(nu, x::Real) = _bessely(nu, float(x)) -_bessely(nu, x::Float16) = Float16(_bessely(nu, Float32(x))) +_bessely(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_bessely(Float32(nu), Float32(x))) function _bessely(nu::T, x::T) where T <: Union{Float32, Float64} isnan(nu) || isnan(x) && return NaN diff --git a/src/modifiedsphericalbessel.jl b/src/modifiedsphericalbessel.jl index 86dd17e..49dfde0 100644 --- a/src/modifiedsphericalbessel.jl +++ b/src/modifiedsphericalbessel.jl @@ -19,7 +19,7 @@ Computes `k_{ν}(x)`, the modified second-kind spherical Bessel function, and of """ sphericalbesselk(nu::Real, x::Real) = _sphericalbesselk(nu, float(x)) -_sphericalbesselk(nu, x::Float16) = Float16(_sphericalbesselk(nu, Float32(x))) +_sphericalbesselk(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_sphericalbesselk(Float32(nu), Float32(x))) function _sphericalbesselk(nu, x::T) where T <: Union{Float32, Float64} if ~isfinite(x) @@ -39,7 +39,7 @@ function _sphericalbesselk(nu, x::T) where T <: Union{Float32, Float64} _nu = ifelse(nu Date: Tue, 4 Oct 2022 10:35:55 -0400 Subject: [PATCH 03/11] add hankel functions --- src/hankel.jl | 19 +++++++++++++++++++ test/hankel_test.jl | 5 +++++ 2 files changed, 24 insertions(+) diff --git a/src/hankel.jl b/src/hankel.jl index 956ce9d..6650973 100644 --- a/src/hankel.jl +++ b/src/hankel.jl @@ -194,6 +194,25 @@ function besselh(nu::Real, k::Integer, x) end end +function besselh(nu::AbstractRange, k::Integer, x::T) where T + (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1")) + if nu[end] < x + out = Vector{Complex{T}}(undef, length(nu)) + out[1], out[2] = besselh(nu[1], k, x), besselh(nu[2], k, x) + return besselj_up_recurrence!(out, x, nu) + else + Jn = besselj(nu, x) + Yn = bessely(nu, x) + if k == 1 + return complex.(Jn, Yn) + elseif k == 2 + return complex.(Jn, -Yn) + else + throw(ArgumentError("k must be 1 or 2")) + end + end +end + """ hankelh1(nu, x) Bessel function of the third kind of order `nu`, ``H^{(1)}_\\nu(x)``. diff --git a/test/hankel_test.jl b/test/hankel_test.jl index 01fde8f..9ccd0a5 100644 --- a/test/hankel_test.jl +++ b/test/hankel_test.jl @@ -26,3 +26,8 @@ v, x = 14.3, 29.4 @test isapprox(hankelh2(v, x), SpecialFunctions.hankelh2(v, x), rtol=2e-13) @test isapprox(besselh(v, 1, x), SpecialFunctions.besselh(v, 1, x), rtol=2e-13) @test isapprox(besselh(v, 2, x), SpecialFunctions.besselh(v, 2, x), rtol=2e-13) + +@test isapprox(hankelh1(1:50, 10.0), SpecialFunctions.hankelh1.(1:50, 10.0), rtol=2e-13) +@test isapprox(hankelh1(0.5:1:25.5, 15.0), SpecialFunctions.hankelh1.(0.5:1:25.5, 15.0), rtol=2e-13) +@test isapprox(hankelh1(1:50, 100.0), SpecialFunctions.hankelh1.(1:50, 100.0), rtol=2e-13) +@test isapprox(hankelh2(1:50, 10.0), SpecialFunctions.hankelh2.(1:50, 10.0), rtol=2e-13) From 135473e6e9dc9545eebda95a11ba784ea1d27fb7 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Fri, 7 Oct 2022 12:31:14 -0400 Subject: [PATCH 04/11] support underflow Jnu, Inu --- src/besseli.jl | 28 ++++++++++++++++++++++++++++ src/besselj.jl | 28 +++++++++++++++++++++++++--- test/besseli_test.jl | 2 +- test/besselj_test.jl | 6 +++--- 4 files changed, 57 insertions(+), 7 deletions(-) diff --git a/src/besseli.jl b/src/besseli.jl index fed01c5..4d1c5e4 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -213,6 +213,34 @@ function _besseli(nu::AbstractRange, x::T) where T return besselk_down_recurrence!(out, x, nu) end +function _besseli(nu::AbstractRange, x::T) where T + (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1")) + len = length(nu) + isone(len) && return [besseli(nu[1], x)] + out = zeros(T, len) + k = len + inu = zero(T) + while abs(inu) < floatmin(T) + if besseli_underflow_check(nu[k], x) + inu = zero(T) + else + inu = _besseli(nu[k], x) + end + out[k] = inu + k -= 1 + k < 1 && break + end + if k > 1 + out[k] = _besseli(nu[k], x) + out[1:k+1] = besselk_down_recurrence!(out[1:k+1], x, nu[1:k+1]) + return out + else + return out + end +end + +besseli_underflow_check(nu, x::T) where T = nu > 140 + T(1.45)*x + 53*Base.Math._approx_cbrt(x) + """ besseli_positive_args(nu, x::T) where T <: Union{Float32, Float64} diff --git a/src/besselj.jl b/src/besselj.jl index 34e920e..4d3f9ed 100644 --- a/src/besselj.jl +++ b/src/besselj.jl @@ -257,16 +257,38 @@ end function _besselj(nu::AbstractRange, x::T) where T (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1")) - out = Vector{T}(undef, length(nu)) + len = length(nu) + isone(len) && return [besselj(nu[1], x)] + + out = zeros(T, len) if nu[end] < x out[1], out[2] = _besselj(nu[1], x), _besselj(nu[2], x) return besselj_up_recurrence!(out, x, nu) else - out[end-1], out[end] = _besselj(nu[end-1], x), _besselj(nu[end], x) - return besselj_down_recurrence!(out, x, nu) + k = len + jn = zero(T) + while abs(jn) < floatmin(T) + if besselj_underflow_check(nu[k], x) + jn = zero(T) + else + jn = _besselj(nu[k], x) + end + out[k] = jn + k -= 1 + k < 1 && break + end + if k > 1 + out[k] = _besselj(nu[k], x) + out[1:k+1] = besselj_down_recurrence!(out[1:k+1], x, nu[1:k+1]) + return out + else + return out + end end end +besselj_underflow_check(nu, x::T) where T = nu > 100 + T(1.01)*x + 85*Base.Math._approx_cbrt(x) + """ besselj_positive_args(nu, x::T) where T <: Float64 diff --git a/test/besseli_test.jl b/test/besseli_test.jl index 84dc471..3f00617 100644 --- a/test/besseli_test.jl +++ b/test/besseli_test.jl @@ -107,7 +107,7 @@ for v in nu, xx in x end # test nu_range -@test besseli(0:50, 2.0) ≈ SpecialFunctions.besseli.(0:50, 2.0) rtol=1e-13 +@test besseli(0:250, 2.0) ≈ SpecialFunctions.besseli.(0:250, 2.0) rtol=1e-13 @test besseli(0.5:1:10.5, 2.0) ≈ SpecialFunctions.besseli.(0.5:1:10.5, 2.0) rtol=1e-13 ### need to fix method ambiguities for other functions ###### diff --git a/test/besselj_test.jl b/test/besselj_test.jl index 645cbb0..27d91e2 100644 --- a/test/besselj_test.jl +++ b/test/besselj_test.jl @@ -115,9 +115,9 @@ for v in nu, xx in x end # test nu_range -@test besselj(0:50, 2.0) ≈ SpecialFunctions.besselj.(0:50, 2.0) rtol=1e-11 -@test besselj(0:50, 100.0) ≈ SpecialFunctions.besselj.(0:50, 100.0) rtol=1e-11 -@test besselj(0.5:1:10.5, 2.0) ≈ SpecialFunctions.besselj.(0.5:1:10.5, 2.0) rtol=1e-11 +@test besselj(0:250, 2.0) ≈ SpecialFunctions.besselj.(0:250, 2.0) rtol=1e-11 +@test besselj(0:95, 100.0) ≈ SpecialFunctions.besselj.(0:95, 100.0) rtol=1e-11 +@test besselj(0.5:1:150.5, 2.0) ≈ SpecialFunctions.besselj.(0.5:1:150.5, 2.0) rtol=1e-11 @test besselj(0.5:1:10.5, 40.0) ≈ SpecialFunctions.besselj.(0.5:1:10.5, 40.0) rtol=1e-11 # test Float16 and Float32 From c9ff1c221f73f723a4ccc9d6f27e046c70f4d5c9 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Fri, 7 Oct 2022 12:32:39 -0400 Subject: [PATCH 05/11] delete old code --- src/besseli.jl | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/besseli.jl b/src/besseli.jl index 4d1c5e4..6113e1b 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -206,13 +206,6 @@ function _besseli(nu::Integer, x::T) where T <: Union{Float32, Float64} end end -function _besseli(nu::AbstractRange, x::T) where T - (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1")) - out = Vector{T}(undef, length(nu)) - out[end-1], out[end] = _besseli(nu[end-1], x), _besseli(nu[end], x) - return besselk_down_recurrence!(out, x, nu) -end - function _besseli(nu::AbstractRange, x::T) where T (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1")) len = length(nu) From 30526b897ee0adac94b885a78dbda8c932d2d2b0 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Fri, 7 Oct 2022 12:52:34 -0400 Subject: [PATCH 06/11] fix underflow for besselk --- src/besselk.jl | 27 ++++++++++++++++++++++++--- test/besselk_test.jl | 1 + 2 files changed, 25 insertions(+), 3 deletions(-) diff --git a/src/besselk.jl b/src/besselk.jl index fa60241..9e3596c 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -218,11 +218,32 @@ end function _besselk(nu::AbstractRange, x::T) where T (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1")) - out = Vector{T}(undef, length(nu)) - out[1], out[2] = _besselk(nu[1], x), _besselk(nu[2], x) - return besselk_up_recurrence!(out, x, nu) + len = length(nu) + isone(len) && return [besselk(nu[1], x)] + out = zeros(T, len) + k = 1 + knu = zero(T) + while abs(knu) < floatmin(T) + if besselk_underflow_check(nu[k], x) + knu = zero(T) + else + knu = _besselk(nu[k], x) + end + out[k] = knu + k += 1 + k == len && break + end + if k < len + out[k] = _besselk(nu[k], x) + out[k-1:end] = besselk_up_recurrence!(out[k-1:end], x, nu[k-1:end]) + return out + else + return out + end end +besselk_underflow_check(nu, x::T) where T = nu < T(1.45)*(x - 780) + 45*Base.Math._approx_cbrt(x - 780) + """ besselk_positive_args(x::T) where T <: Union{Float32, Float64} diff --git a/test/besselk_test.jl b/test/besselk_test.jl index e4b60bf..3ca3c9b 100644 --- a/test/besselk_test.jl +++ b/test/besselk_test.jl @@ -129,6 +129,7 @@ end # test nu_range @test besselk(0:50, 2.0) ≈ SpecialFunctions.besselk.(0:50, 2.0) rtol=1e-13 @test besselk(0.5:1:10.5, 12.0) ≈ SpecialFunctions.besselk.(0.5:1:10.5, 12.0) rtol=1e-13 +@test besselk(1:700, 800.0) ≈ SpecialFunctions.besselk.(1:700, 800.0) rtol=1e-13 # test Float16 @test besselk(Int16(10), Float16(1.0)) isa Float16 From abc57aedc92d88bee7787e844ad2793f4a4b8f00 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Fri, 7 Oct 2022 12:57:40 -0400 Subject: [PATCH 07/11] test --- test/besselk_test.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/besselk_test.jl b/test/besselk_test.jl index 3ca3c9b..bf71ed6 100644 --- a/test/besselk_test.jl +++ b/test/besselk_test.jl @@ -129,7 +129,7 @@ end # test nu_range @test besselk(0:50, 2.0) ≈ SpecialFunctions.besselk.(0:50, 2.0) rtol=1e-13 @test besselk(0.5:1:10.5, 12.0) ≈ SpecialFunctions.besselk.(0.5:1:10.5, 12.0) rtol=1e-13 -@test besselk(1:700, 800.0) ≈ SpecialFunctions.besselk.(1:700, 800.0) rtol=1e-13 +@test besselk(1:700, 800.0) ≈ SpecialFunctions.besselk.(1:700, 800.0) # test Float16 @test besselk(Int16(10), Float16(1.0)) isa Float16 From db3f78038181c68e24ca6d2ddea20a56218cb4a1 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 12 Oct 2022 11:00:33 -0400 Subject: [PATCH 08/11] fix recurrence bounds and rebase --- src/recurrence.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/recurrence.jl b/src/recurrence.jl index cca9ff1..c4cf681 100644 --- a/src/recurrence.jl +++ b/src/recurrence.jl @@ -13,7 +13,7 @@ end function besselk_up_recurrence!(out, x::T, nu_range) where T x2 = 2 / x k = 3 - for nu in nu_range[2:end-1] + for nu in nu_range[begin+1:end-1] out[k] = muladd(nu*x2, out[k-1], out[k-2]) k += 1 end @@ -39,7 +39,7 @@ end function besselj_up_recurrence!(out, x::T, nu_range) where T x2 = 2 / x k = 3 - for nu in nu_range[2:end-1] + for nu in nu_range[begin+1:end-1] out[k] = muladd(nu*x2, out[k-1], -out[k-2]) k += 1 end @@ -62,15 +62,13 @@ end function besselj_down_recurrence!(out, x::T, nu_range) where T x2 = 2 / x k = length(nu_range) - 2 - for nu in nu_range[end-1:-1:2] + for nu in nu_range[end-1:-1:begin+1] out[k] = muladd(nu*x2, out[k+1], -out[k+2]) k -= 1 end return out 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); @@ -83,11 +81,11 @@ end end return jnup1, jnu end -=# + function besselk_down_recurrence!(out, x::T, nu_range) where T x2 = 2 / x k = length(nu_range) - 2 - for nu in nu_range[end-1:-1:2] + for nu in nu_range[end-1:-1:begin+1] out[k] = muladd(nu*x2, out[k+1], out[k+2]) k -= 1 end From 5aab92fd98cfdd494246f4ac7dfa17e9e4f9ee7b Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 12 Oct 2022 11:03:52 -0400 Subject: [PATCH 09/11] update news and bump version --- NEWS.md | 5 +++++ Project.toml | 2 +- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index ee07c90..5e03686 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.3 + +### Added +- Add support for nu isa AbstractRange (i.e., `besselj(0:10, 1.0)`) to allow for fast computation of Bessel functions at many orders ([PR #53](https://github.com/JuliaMath/Bessels.jl/pull/53)). + # Version 0.2.2 ### Added diff --git a/Project.toml b/Project.toml index 264beb9..70b4730 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.2" +version = "0.2.3" [compat] julia = "1.6" From 02a9d7a0531977877a91a85fe9e10293eef875b8 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 12 Oct 2022 11:28:35 -0400 Subject: [PATCH 10/11] add an example in the readme --- README.md | 40 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/README.md b/README.md index 05e0291..48c62c4 100644 --- a/README.md +++ b/README.md @@ -120,6 +120,46 @@ julia> besselk0(x) julia> besselk1(x) 1.6750295538365835e-6 ``` +## Support for sequence of orders + +We also provide support for `besselj(nu::M, x::T)`, `bessely(nu::M, x::T)`, `besseli(nu::M, x::T)`, `besselk(nu::M, x::T)`, `besseli(nu::M, x::T)`, `besselh(nu::M, k, x::T)` when `M` is some `AbstractRange` and `T` is some float. + +```julia +julia> besselj(0:10, 1.0) +11-element Vector{Float64}: + 0.7651976865579666 + 0.44005058574493355 + 0.11490348493190049 + 0.019563353982668407 + 0.0024766389641099553 + 0.00024975773021123444 + 2.0938338002389273e-5 + 1.5023258174368085e-6 + 9.422344172604502e-8 + 5.249250179911876e-9 + 2.630615123687453e-10 +``` + +In general, this provides a fast way to generate a sequence of Bessel functions for many orders. +```julia +julia> @btime besselj(0:100, 50.0) + 443.328 ns (2 allocations: 1.75 KiB) +``` +This function will allocate so it is recommended that you calculate the Bessel functions at the top level of your function outside any hot loop. For example, + +```julia +function bessel_sequence(x, orders) + J_nu = besselj(orders, x) + out = zero(x) + for i in eachindex(J_nu) + out += sin(x*i)*J_nu[i] + end + return out +end + +julia> bessel_sequence(10.2, 1:400) +0.11404996570230919 +``` ### Support for negative arguments From 4a59336e9bf374af257625d41978f478ddfbe7ec Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 12 Oct 2022 12:27:47 -0400 Subject: [PATCH 11/11] fix allocation and readme --- README.md | 16 +--------------- src/besseli.jl | 3 ++- src/besselj.jl | 3 ++- src/besselk.jl | 3 ++- 4 files changed, 7 insertions(+), 18 deletions(-) diff --git a/README.md b/README.md index 48c62c4..21f4447 100644 --- a/README.md +++ b/README.md @@ -145,21 +145,7 @@ In general, this provides a fast way to generate a sequence of Bessel functions julia> @btime besselj(0:100, 50.0) 443.328 ns (2 allocations: 1.75 KiB) ``` -This function will allocate so it is recommended that you calculate the Bessel functions at the top level of your function outside any hot loop. For example, - -```julia -function bessel_sequence(x, orders) - J_nu = besselj(orders, x) - out = zero(x) - for i in eachindex(J_nu) - out += sin(x*i)*J_nu[i] - end - return out -end - -julia> bessel_sequence(10.2, 1:400) -0.11404996570230919 -``` +This function will allocate so it is recommended that you calculate the Bessel functions at the top level of your function outside any hot loop. ### Support for negative arguments diff --git a/src/besseli.jl b/src/besseli.jl index 877abe9..fc5d3fc 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -225,7 +225,8 @@ function _besseli(nu::AbstractRange, x::T) where T end if k > 1 out[k] = _besseli(nu[k], x) - out[1:k+1] = besselk_down_recurrence!(out[1:k+1], x, nu[1:k+1]) + tmp = @view out[1:k+1] + out[1:k+1] = besselk_down_recurrence!(tmp, x, nu[1:k+1]) return out else return out diff --git a/src/besselj.jl b/src/besselj.jl index 4d3f9ed..b2959e7 100644 --- a/src/besselj.jl +++ b/src/besselj.jl @@ -279,7 +279,8 @@ function _besselj(nu::AbstractRange, x::T) where T end if k > 1 out[k] = _besselj(nu[k], x) - out[1:k+1] = besselj_down_recurrence!(out[1:k+1], x, nu[1:k+1]) + tmp = @view out[1:k+1] + besselj_down_recurrence!(tmp, x, nu[1:k+1]) return out else return out diff --git a/src/besselk.jl b/src/besselk.jl index 6dce47e..62d6c16 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -235,7 +235,8 @@ function _besselk(nu::AbstractRange, x::T) where T end if k < len out[k] = _besselk(nu[k], x) - out[k-1:end] = besselk_up_recurrence!(out[k-1:end], x, nu[k-1:end]) + tmp = @view out[k-1:end] + out[k-1:end] = besselk_up_recurrence!(tmp, x, nu[k-1:end]) return out else return out