support underflow Jnu, Inu

Author: heltonmc

src/besseli.jl

@@ -213,6 +213,34 @@ function _besseli(nu::AbstractRange, x::T) where T
     return besselk_down_recurrence!(out, x, nu) 
 end
 
+function _besseli(nu::AbstractRange, x::T) where T
+    (nu[1] >= 0 && step(nu) == 1) || throw(ArgumentError("nu must be >= 0 with step(nu)=1"))
+    len = length(nu)
+    isone(len) && return [besseli(nu[1], x)]
+    out = zeros(T, len)
+    k = len
+    inu = zero(T)
+    while abs(inu) < floatmin(T)
+        if besseli_underflow_check(nu[k], x)
+            inu = zero(T)
+        else
+            inu = _besseli(nu[k], x)
+        end
+        out[k] = inu
+        k -= 1
+        k < 1 && break
+    end
+    if k > 1
+        out[k] = _besseli(nu[k], x)
+        out[1:k+1] = besselk_down_recurrence!(out[1:k+1], x, nu[1:k+1])
+        return out
+    else
+        return out
+    end
+end
+
+besseli_underflow_check(nu, x::T) where T = nu > 140 + T(1.45)*x + 53*Base.Math._approx_cbrt(x)
+
 """
     besseli_positive_args(nu, x::T) where T <: Union{Float32, Float64}
 

src/besselj.jl

@@ -257,16 +257,38 @@ end
 
 function _besselj(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))
+    len = length(nu)
+    isone(len) && return [besselj(nu[1], x)]
+
+    out = zeros(T, len)
     if nu[end] < x
         out[1], out[2] = _besselj(nu[1], x), _besselj(nu[2], x)
         return besselj_up_recurrence!(out, x, nu)
     else
-        out[end-1], out[end] = _besselj(nu[end-1], x), _besselj(nu[end], x)
-        return besselj_down_recurrence!(out, x, nu)
+        k = len
+        jn = zero(T)
+        while abs(jn) < floatmin(T)
+            if besselj_underflow_check(nu[k], x)
+                jn = zero(T)
+            else
+                jn = _besselj(nu[k], x)
+            end
+            out[k] = jn
+            k -= 1
+            k < 1 && break
+        end
+        if k > 1
+            out[k] = _besselj(nu[k], x)
+            out[1:k+1] = besselj_down_recurrence!(out[1:k+1], x, nu[1:k+1])
+            return out
+        else
+            return out
+        end
     end
 end
 
+besselj_underflow_check(nu, x::T) where T = nu > 100 + T(1.01)*x + 85*Base.Math._approx_cbrt(x)
+
 """
     besselj_positive_args(nu, x::T) where T <: Float64
 

test/besseli_test.jl

@@ -107,7 +107,7 @@ for v in nu, xx in x
 end
 
 # test nu_range
-@test besseli(0:50, 2.0) ≈ SpecialFunctions.besseli.(0:50, 2.0) rtol=1e-13
+@test besseli(0:250, 2.0) ≈ SpecialFunctions.besseli.(0:250, 2.0) rtol=1e-13
 @test besseli(0.5:1:10.5, 2.0) ≈ SpecialFunctions.besseli.(0.5:1:10.5, 2.0) rtol=1e-13
 
 ### need to fix method ambiguities for other functions ###### 

test/besselj_test.jl

@@ -115,9 +115,9 @@ for v in nu, xx in x
 end
 
 # test nu_range
-@test besselj(0:50, 2.0) ≈ SpecialFunctions.besselj.(0:50, 2.0) rtol=1e-11
-@test besselj(0:50, 100.0) ≈ SpecialFunctions.besselj.(0:50, 100.0) rtol=1e-11
-@test besselj(0.5:1:10.5, 2.0) ≈ SpecialFunctions.besselj.(0.5:1:10.5, 2.0) rtol=1e-11
+@test besselj(0:250, 2.0) ≈ SpecialFunctions.besselj.(0:250, 2.0) rtol=1e-11
+@test besselj(0:95, 100.0) ≈ SpecialFunctions.besselj.(0:95, 100.0) rtol=1e-11
+@test besselj(0.5:1:150.5, 2.0) ≈ SpecialFunctions.besselj.(0.5:1:150.5, 2.0) rtol=1e-11
 @test besselj(0.5:1:10.5, 40.0) ≈ SpecialFunctions.besselj.(0.5:1:10.5, 40.0) rtol=1e-11
 
 # test Float16 and Float32