Merge pull request #53 from JuliaMath/sequence

Return sequence of Bessel functions

Author: heltonmc

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

Project.toml

@@ -1,7 +1,7 @@
 name = "Bessels"
 uuid = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
 authors = ["Michael Helton <[email protected]> and contributors"]
-version = "0.2.2"
+version = "0.2.3"
 
 [compat]
 julia = "1.6"

README.md

@@ -120,6 +120,32 @@ 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.
 
 ### Support for negative arguments
 

src/besseli.jl

@@ -164,11 +164,12 @@ end
 
 Modified Bessel function of the second kind of order nu, ``I_{nu}(x)``.
 """
-besseli(nu::Real, x::Real) = _besseli(nu, float(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, 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 +206,35 @@ 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"))
+    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)
+        tmp = @view out[1:k+1]
+        out[1:k+1] = besselk_down_recurrence!(tmp, 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}
 
@@ -232,7 +262,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)

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)))
+_besselj(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_besselj(Float32(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,41 @@ 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"))
+    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
+        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)
+            tmp = @view out[1:k+1]
+            besselj_down_recurrence!(tmp, 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
 

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)))
+_besselk(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_besselk(Float32(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,35 @@ 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"))
+    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)
+        tmp = @view out[k-1:end]
+        out[k-1:end] = besselk_up_recurrence!(tmp, 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}
 

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)))
+_bessely(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_bessely(Float32(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
 #####

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)``.

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<zero(nu), -one(nu)-nu, nu)
         return sphericalbesselk_int(Int(_nu), x)
     else
-        return inv(SQPIO2(T)*sqrt(x))*besselk(nu+1/2, x)
+        return inv(SQPIO2(T)*sqrt(x))*besselk(nu+T(1)/2, x)
     end
 end
 sphericalbesselk_cutoff(nu) = nu < 41.5
@@ -65,15 +65,15 @@ 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)))
+_sphericalbesseli(nu::Union{Int16, Float16}, x::Union{Int16, Float16}) = Float16(_sphericalbesseli(Float32(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)
+    return SQPIO2(T)*besseli(nu+T(1)/2, x) / sqrt(x)
 end
 
 function _sphericalbesseli_small_orders(nu::Integer, x::T) where T
@@ -86,7 +86,7 @@ function _sphericalbesseli_small_orders(nu::Integer, x::T) where T
     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)
+    return SQPIO2(T)*besseli(nu+T(1)/2, x) / sqrt(x)
 end
 
 function sphericalbesseli_small_args(nu, x::T) where T