Merge pull request #40 from JuliaMath/bessely_nan_error

Fix #35 and correct for NaN return for small arguments in bessely

Author: heltonmc

NEWS.md

@@ -18,6 +18,8 @@ For bug fixes, performance enhancements, or fixes to unexported functions we wil
 
 ### Fixed
  - fix cutoff in `bessely` to not return error for integer orders and small arguments (#33)
+ - fix NaN return for small arguments (issue [#35]) in bessely (#40)
+ - allow calling with integer argument and add float16 support (#40)
 
 # Version 0.1.0
 

src/besseli.jl

@@ -164,8 +164,12 @@ end
 
 Modified Bessel function of the second kind of order nu, ``I_{nu}(x)``.
 """
-function besseli(nu::Real, x::T) where T
-    isinteger(nu) && return besseli(Int(nu), x)
+besseli(nu::Real, x::Real) = _besseli(nu, float(x))
+
+_besseli(nu, x::Float16) = Float16(_besseli(nu, Float32(x)))
+
+function _besseli(nu, x::T) where T <: Union{Float32, Float64}
+    isinteger(nu) && return _besseli(Int(nu), x)
     abs_nu = abs(nu)
     abs_x = abs(x)
 
@@ -187,7 +191,7 @@ function besseli(nu::Real, x::T) where T
         end
     end
 end
-function besseli(nu::Integer, x::T) where T
+function _besseli(nu::Integer, x::T) where T <: Union{Float32, Float64}
     abs_nu = abs(nu)
     abs_x = abs(x)
     sg = iseven(abs_nu) ? 1 : -1
@@ -225,7 +229,11 @@ end
 Scaled modified Bessel function of the first kind of order nu, ``I_{nu}(x)*e^{-x}``.
 Nu must be real.
 """
-function besselix(nu, x::T) where T <: Union{Float32, Float64}
+besselix(nu::Real, x::Real) = _besselix(nu, float(x))
+
+_besselix(nu, x::Float16) = Float16(_besselix(nu, Float32(x)))
+
+function _besselix(nu, x::T) where T <: Union{Float32, Float64}
     iszero(nu) && return besseli0x(x)
     isone(nu) && return besseli1x(x)
 

src/besselj.jl

@@ -192,8 +192,18 @@ end
 ##### Generic routine for `besselj`
 #####
 
-function besselj(nu::Real, x::T) where T
-    isinteger(nu) && return besselj(Int(nu), x)
+"""
+    besselj(nu, x::T) where T <: Union{Float32, Float64}
+
+Bessel function of the first kind of order nu, ``J_{nu}(x)``.
+nu and x must be real where nu and x can be positive or negative.
+"""
+besselj(nu::Real, x::Real) = _besselj(nu, float(x))
+
+_besselj(nu, x::Float16) = Float16(_besselj(nu, Float32(x)))
+
+function _besselj(nu, x::T) where T <: Union{Float32, Float64}
+    isinteger(nu) && return _besselj(Int(nu), x)
     abs_nu = abs(nu)
     abs_x = abs(x)
 
@@ -218,7 +228,7 @@ function besselj(nu::Real, x::T) where T
     end
 end
 
-function besselj(nu::Integer, x::T) where T
+function _besselj(nu::Integer, x::T) where T <: Union{Float32, Float64}
     abs_nu = abs(nu)
     abs_x = abs(x)
     sg = iseven(abs_nu) ? 1 : -1

src/besselk.jl

@@ -185,8 +185,12 @@ end
 
 Modified Bessel function of the second kind of order nu, ``K_{nu}(x)``.
 """
-function besselk(nu::Real, x::T) where T
-    isinteger(nu) && return besselk(Int(nu), x)
+besselk(nu::Real, x::Real) = _besselk(nu, float(x))
+
+_besselk(nu, x::Float16) = Float16(_besselk(nu, Float32(x)))
+
+function _besselk(nu, x::T) where T <: Union{Float32, Float64}
+    isinteger(nu) && return _besselk(Int(nu), x)
     abs_nu = abs(nu)
     abs_x = abs(x)
 
@@ -197,7 +201,7 @@ function besselk(nu::Real, x::T) where T
         #return cispi(-abs_nu)*besselk_positive_args(abs_nu, abs_x) - besseli_positive_args(abs_nu, abs_x) * im * π
     end
 end
-function besselk(nu::Integer, x::T) where T
+function _besselk(nu::Integer, x::T) where T <: Union{Float32, Float64}
     abs_nu = abs(nu)
     abs_x = abs(x)
     sg = iseven(abs_nu) ? 1 : -1
@@ -239,7 +243,11 @@ end
 
 Scaled modified Bessel function of the second kind of order nu, ``K_{nu}(x)*e^{x}``.
 """
-function besselkx(nu, x::T) where T <: Union{Float32, Float64}
+besselkx(nu::Real, x::Real) = _besselkx(nu, float(x))
+
+_besselkx(nu, x::Float16) = Float16(_besselkx(nu, Float32(x)))
+
+function _besselkx(nu, x::T) where T <: Union{Float32, Float64}
     # dispatch to avoid uniform expansion when nu = 0 
     iszero(nu) && return besselk0x(x)
 

src/bessely.jl

@@ -235,9 +235,13 @@ end
 Bessel function of the second kind of order nu, ``Y_{nu}(x)``.
 nu and x must be real where nu and x can be positive or negative.
 """
-function bessely(nu::Real, x::T) where T
+bessely(nu::Real, x::Real) = _bessely(nu, float(x))
+
+_bessely(nu, x::Float16) = Float16(_bessely(nu, Float32(x)))
+
+function _bessely(nu, x::T) where T <: Union{Float32, Float64}
     isnan(nu) || isnan(x) && return NaN
-    isinteger(nu) && return bessely(Int(nu), x)
+    isinteger(nu) && return _bessely(Int(nu), x)
     abs_nu = abs(nu)
     abs_x = abs(x)
 
@@ -260,7 +264,7 @@ function bessely(nu::Real, x::T) where T
         end
     end
 end
-function bessely(nu::Integer, x::T) where T
+function _bessely(nu::Integer, x::T) where T <: Union{Float32, Float64}
     abs_nu = abs(nu)
     abs_x = abs(x)
     sg = iseven(abs_nu) ? 1 : -1
@@ -336,6 +340,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 +354,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), a)
+
     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/besseli_test.jl

@@ -80,6 +80,8 @@ end
 #Float 64
 m = 0:1:200; x = 0.1:0.5:150.0
 @test besseli(10, 1.0) isa Float64
+@test besseli(10, Float16(1.0)) isa Float16
+
 @test besseli(2, 80.0) isa Float64
 @test besseli(112, 80.0) isa Float64
 t = [besseli(m, x) for m in m, x in x]
@@ -90,6 +92,7 @@ t = [besseli(m, x) for m in m, x in x]
 t = [besselix(m, x) for m in m, x in x]
 @test t[10] isa Float64
 @test t ≈ [SpecialFunctions.besselix(m, x) for m in m, x in x]
+@test besselix(10, Float16(1.0)) isa Float16
 
 ## Tests for besselk
 

test/besselj_test.jl

@@ -110,6 +110,9 @@ for v in nu, xx in x
     @test isapprox(Bessels.besseljy_positive_args(v, xx)[1], Float64(sf), rtol=5e-11)
 end
 
+# test Float16
+@test besselj(10, Float16(1.0)) isa Float16
+
 ## test large arguments
 @test isapprox(besselj(10.0, 150.0), SpecialFunctions.besselj(10.0, 150.0), rtol=1e-12)
 

test/besselk_test.jl

@@ -120,6 +120,10 @@ for v in nu, xx in x
     @test isapprox(besselk(v, xx), Float64(sf), rtol=2e-13)
 end
 
+# test Float16
+@test besselk(10, Float16(1.0)) isa Float16
+@test besselkx(10, Float16(1.0)) isa Float16
+
 # test Inf
 @test iszero(besselk(2, Inf))
 

test/bessely_test.jl

@@ -89,6 +89,16 @@ 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 Float16
+@test bessely(10, Float16(1.0)) isa Float16
+
+# 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