Special half-integer branch for besselk.

Author: cgeoga

src/besselk.jl

@@ -226,6 +226,9 @@ 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)
 
+    # check if nu is a half-integer:
+    isinteger(nu-1/2) && return besselk_vhalfint(nu, x)
+
     # use uniform debye expansion if x or nu is large
     besselik_debye_cutoff(nu, x) && return besselk_large_orders(nu, x)
 
@@ -424,3 +427,23 @@ 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
+
+
+"""
+    besselk_vhalfint(nu, x::T) where T <: {Float32, Float64}
+
+Computes `K_{ν}(x)` when `v + 1/2` is an integer using the fact that the
+asymptotic expansion actually terminates and is exact for those specific `v` values.
+"""
+function besselk_vhalfint(v, x::T) where T
+    v = abs(v)
+    invx = inv(x)
+    b0 = b1 = SQRT_PID2(T)*sqrt(invx)*exp(-x) 
+    twodx = 2*invx
+    _v  = T(1/2)
+    while _v < v
+        b0, b1 = b1, muladd(b1, twodx*_v, b0)
+        _v += one(T)
+    end
+    b1
+end

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

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 in (-3/2, -1/2, 1/2, 3/2, 5/2, 7/2, 9/2)
+  @test besselk(v, 7.0) ≈ SpecialFunctions.besselk(v, 7.0)
+  @test besselk(Float32(v), Float32(7.0)) ≈ SpecialFunctions.besselk(Float32(v), Float32(7.0))
+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