fix some spacing, overflow, and cutoff, add docs

Author: heltonmc

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
-

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

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<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
-

test/sphericalbessel_test.jl

@@ -8,6 +8,8 @@ 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))
 
 # test zero
 @test isone(Bessels.sphericalbesselj(0, 0.0))
@@ -20,33 +22,48 @@ 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))
+
 # 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))
+
 # 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))
+
 # 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
 
 for x in 0.5:1.0:100.0, v in [0, 1, 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)
 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)
 end
 
 @test isapprox(Bessels.sphericalbessely(270, 240.0), SpecialFunctions.sphericalbessely(270, 240.0), rtol=3e-12)