Merge pull request #46 from cgeoga/master

Special half-integer branch for `besselk`

Author: heltonmc

src/Bessels.jl

@@ -41,6 +41,7 @@ include("bessely.jl")
 include("hankel.jl")
 include("airy.jl")
 include("sphericalbessel.jl")
+include("modifiedsphericalbessel.jl")
 
 include("Float128/besseli.jl")
 include("Float128/besselj.jl")

src/besselk.jl

@@ -225,9 +225,15 @@ 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)
 
     # use uniform debye expansion if x or nu is large
-    besselik_debye_cutoff(nu, x) && return besselk_large_orders(nu, x)
+    debye_cut && 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]
@@ -424,3 +430,4 @@ 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,3 +284,6 @@ 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

@@ -0,0 +1,40 @@
+
+"""
+    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))
+
+function _sphericalbesselk(nu, x::T) where T
+    isnan(x) && return NaN
+    if isinteger(nu) && nu < 41.5
+        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 
+        # k_{-n} = (...)*K_{-n     + 1/2}
+        #        = (...)*K_{|n|    - 1/2}
+        #        = (...)*K_{|n|-1  + 1/2}
+        #        = 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)
+    end
+end
+
+function sphericalbesselk_int(v::Int, x)
+    b0 = inv(x)
+    b1 = (x+one(x))/(x*x)
+    iszero(v) && return b0*exp(-x)
+    _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
+end
+

test/besselk_test.jl

@@ -67,6 +67,12 @@ k1x_32 = besselk1x.(Float32.(x))
 @test besselk(100, 3.9) ≈ SpecialFunctions.besselk(100, 3.9)
 @test besselk(100, 234.0) ≈ SpecialFunctions.besselk(100, 234.0)
 
+# test half-integer orders:
+for (v,x) in Iterators.product(-35/2:1.0:81/2, range(0.0, 30.0, length=51)[2:end])
+  @test besselk(v, x) ≈ SpecialFunctions.besselk(v, x)
+  @test besselk(Float32(v), Float32(x)) ≈ SpecialFunctions.besselk(Float32(v), Float32(x))
+end
+
 # test small arguments and order
 m = 0:40; x = [1e-6; 1e-4; 1e-3; 1e-2; 0.1; 1.0:2.0:500.0]
 for m in m, x in x

test/sphericalbessel_test.jl

@@ -58,3 +58,9 @@ v, x = -4.0, 5.6
 v, x = -6.8, 15.6
 @test isapprox(Bessels.sphericalbesselj(v, x), 0.04386355397884301866595, rtol=3e-12)
 @test isapprox(Bessels.sphericalbessely(v, x), 0.05061013363904335437354, rtol=3e-12)
+
+# test for negative order of spherical modified besselk in the special integer
+# routine:
+for v in 1:10
+  @test Bessels.sphericalbesselk(-v, 1.1) ≈ Bessels.sphericalbesselk(v-1, 1.1)
+end