check for underflow

Author: heltonmc

src/bessely.jl

@@ -336,6 +336,8 @@ end
 # this works well for small arguments x < 7.0 for rel. error ~1e-14
 # this also works well for nu > 1.35x - 4.5
 # for nu > 25 more cancellation occurs near integer values
+# There could be premature underflow when (x/2)^v == 0.
+# It might be better to use logarithms (when we get loggamma julia implementation)
 """
     bessely_power_series(nu, x::T) where T <: Float64
 
@@ -348,6 +350,9 @@ function bessely_power_series(v, x::T) where T
     out = zero(T)
     out2 = zero(T)
     a = (x/2)^v
+    # check for underflow and return limit for small arguments
+    iszero(a) && return -T(Inf)
+
     b = inv(a)
     a /= gamma(v + one(T))
     b /= gamma(-v + one(T))

src/recurrence.jl

@@ -15,14 +15,19 @@ end
 # outputs both (bessel(x, nu_end), bessel(x, nu_end+1)
 # x = 0.1; y0 = bessely0(x); y1 = bessely1(x);
 # besselj_up_recurrence(x, y1, y0, 1, 5) will give bessely(5, x)
-@inline function besselj_up_recurrence(x, jnu, jnum1, nu_start, nu_end)
+@inline function besselj_up_recurrence(x::T, jnu, jnum1, nu_start, nu_end) where T
     x2 = 2 / x
     while nu_start < nu_end + 0.5 # avoid inexact floating points when nu isa float
         jnum1, jnu = jnu, muladd(nu_start*x2, jnu, -jnum1)
         nu_start += 1
     end
-    return jnum1, jnu
+    # need to check if NaN resulted during loop from subtraction of infinities
+    # this could happen if x is very small and nu is large which eventually results in overflow (-> -Inf)
+    # we use this routine to generate bessely(nu::Int, x) in the forward direction starting from bessely0, bessely1
+    # NaN inputs need to be checked before getting to this routine 
+    return isnan(jnum1) ? (-T(Inf), -T(Inf)) : (jnum1, jnu)
 end
+
 # backward recurrence relation for besselj and bessely
 # outputs both (bessel(x, nu_end), bessel(x, nu_end-1)
 # x = 0.1; j10 = besselj(10, x); j11 = besselj(11, x);

test/bessely_test.jl

@@ -89,6 +89,13 @@ for v in nu, xx in x
     @test isapprox(Bessels.besseljy_positive_args(v, xx)[2], SpecialFunctions.bessely(v, xx), rtol=5e-12)
 end
 
+# test limits for small arguments see https://github.com/JuliaMath/Bessels.jl/issues/35
+@test bessely(185.0, 1.01) == -Inf
+@test bessely(185.5, 1.01) == -Inf
+@test bessely(2.0, 1e-300) == -Inf
+@test bessely(4.0, 1e-80) == -Inf
+@test bessely(4.5, 1e-80) == -Inf
+
 # need to test accuracy of negative orders and negative arguments and all combinations within
 # SpecialFunctions.jl doesn't provide these so will hand check against hard values
 # values taken from https://keisan.casio.com/exec/system/1180573474 which match mathematica