Merge pull request #47 from JuliaMath/cgeoga/master

Update sphericalbesselk and add sphericalbesseli

Author: heltonmc

NEWS.md

@@ -11,7 +11,9 @@ 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) 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))

Project.toml

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

src/besseli.jl

@@ -170,6 +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)
+    ~isfinite(x) && return x
     abs_nu = abs(nu)
     abs_x = abs(x)
 
@@ -192,6 +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}
+    ~isfinite(x) && return x
     abs_nu = abs(nu)
     abs_x = abs(x)
     sg = iseven(abs_nu) ? 1 : -1
@@ -211,7 +213,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)

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
 #
 
 """
@@ -226,17 +228,14 @@ 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) && sphericalbesselk_cutoff(nu)) && return sphericalbesselk_int(Int(nu-1/2), x)*SQPIO2(T)*sqrt(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 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
-    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]
@@ -260,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)
 
@@ -434,13 +436,13 @@ 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)
+#####
 
-"""
-  besselk_asymptoticexp(v, x, order)
-
-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. The default order of the expansion in `Bessels.besselk` is 10.
-"""
-besselk_asexp_cutoff(nu, x::T) where T = (nu < 20.0) && (x > ASEXP_CUTOFF(T))
+# 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)

src/constants.jl

@@ -284,9 +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
-
-const ASEXP_CUTOFF(::Type{Float64}) = 35.0
-const ASEXP_CUTOFF(::Type{Float32}) = 35.0f0 # temporary

src/modifiedsphericalbessel.jl

@@ -1,40 +1,98 @@
-
+#                           Modified Spherical Bessel functions
+#
+#                                 sphericalbesseli(nu, x), sphericalbesselk(nu, x)
+#
+#    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) [2].
+#   
+# [1] https://mathworld.wolfram.com/ModifiedBesselFunctionoftheFirstKind.html 
+# [2] https://mathworld.wolfram.com/ModifiedSphericalBesselFunctionoftheSecondKind.html
+#
 """
     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<zero(nu), -one(nu)-nu, nu)
         return sphericalbesselk_int(Int(_nu), x)
     else
-        return inv(SQRT_PID2(T)*sqrt(x))*besselk(nu+1/2, x)
+        return inv(SQPIO2(T)*sqrt(x))*besselk(nu+1/2, x)
     end
 end
+sphericalbesselk_cutoff(nu) = nu < 41.5
 
 function sphericalbesselk_int(v::Int, x)
-    b0 = inv(x)
-    b1 = (x+one(x))/(x*x)
-    iszero(v) && return b0*exp(-x)
+    xinv = inv(x)
+    b0 = exp(-x) * xinv
+    b1 = b0 * (x + one(x)) * xinv
+    iszero(v) && return b0
     _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
+    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

test/sphericalbessel_test.jl

@@ -8,6 +8,16 @@ x = 1e-15
 @test Bessels.sphericalbessely(1, x) ≈ SpecialFunctions.sphericalbessely(1, x)
 @test Bessels.sphericalbessely(5.5, x) ≈ SpecialFunctions.sphericalbessely(5.5, x)
 @test Bessels.sphericalbessely(10, x) ≈ SpecialFunctions.sphericalbessely(10, x)
+@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))
@@ -20,33 +30,77 @@ x = 1e-15
 @test Bessels.sphericalbessely(10, 0.0) == -Inf
 @test Bessels.sphericalbessely(290, 0.0) == -Inf
 
+@test isinf(Bessels.sphericalbesselk(0, 0.0))
+@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))
 @test iszero(Bessels.sphericalbessely(3, Inf))
 @test iszero(Bessels.sphericalbessely(4.5, Inf))
 
+@test iszero(Bessels.sphericalbesselk(0, Inf))
+@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))
 @test isnan(Bessels.sphericalbessely(1.4, NaN))
 @test isnan(Bessels.sphericalbessely(4.0, NaN))
 
+@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
 @test Bessels.sphericalbesselj(1.4f0, 1.2f0) isa Float32
 @test Bessels.sphericalbessely(1.4f0, 1.2f0) isa Float32
 
+@test Bessels.sphericalbesselk(Float16(1.4), Float16(1.2)) isa Float16
+@test Bessels.sphericalbesselk(1.0f0, 1.2f0) isa Float32
+
+@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.0:100.0, v in [0, 1, 5.5, 8.2, 10]
+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)