fix underflow for besselk

Author: heltonmc

src/besselk.jl

@@ -218,11 +218,32 @@ 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"))
-    out = Vector{T}(undef, length(nu))
-    out[1], out[2] = _besselk(nu[1], x), _besselk(nu[2], x)
-    return besselk_up_recurrence!(out, x, nu)
+    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)
+        out[k-1:end] = besselk_up_recurrence!(out[k-1:end], 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}
 

test/besselk_test.jl

@@ -129,6 +129,7 @@ end
 # test nu_range
 @test besselk(0:50, 2.0) ≈ SpecialFunctions.besselk.(0:50, 2.0) rtol=1e-13
 @test besselk(0.5:1:10.5, 12.0) ≈ SpecialFunctions.besselk.(0.5:1:10.5, 12.0) rtol=1e-13
+@test besselk(1:700, 800.0) ≈ SpecialFunctions.besselk.(1:700, 800.0) rtol=1e-13
 
 # test Float16
 @test besselk(Int16(10), Float16(1.0)) isa Float16