From b167f1ab42e19dc3845891adcdc003cfb65656d8 Mon Sep 17 00:00:00 2001
From: Michael Helton <heltonmc@protonmail.com>
Date: Tue, 9 Aug 2022 16:32:30 -0400
Subject: [PATCH 1/6] check for underflow

---
 src/bessely.jl       | 5 +++++
 src/recurrence.jl    | 9 +++++++--
 test/bessely_test.jl | 7 +++++++
 3 files changed, 19 insertions(+), 2 deletions(-)

diff --git a/src/bessely.jl b/src/bessely.jl
index 81cf8c2..a05f554 100644
--- a/src/bessely.jl
+++ b/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))
diff --git a/src/recurrence.jl b/src/recurrence.jl
index fd808ca..c959537 100644
--- a/src/recurrence.jl
+++ b/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);
diff --git a/test/bessely_test.jl b/test/bessely_test.jl
index c2c9971..e563e44 100644
--- a/test/bessely_test.jl
+++ b/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

From fed39ca1aea4cd89bc46709f52a37102eb970794 Mon Sep 17 00:00:00 2001
From: Michael Helton <heltonmc@protonmail.com>
Date: Tue, 9 Aug 2022 16:46:51 -0400
Subject: [PATCH 2/6] add news

---
 NEWS.md        | 1 +
 src/bessely.jl | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/NEWS.md b/NEWS.md
index 467c8ab..7f1c4d7 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -18,6 +18,7 @@ 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)
 
 # Version 0.1.0
 
diff --git a/src/bessely.jl b/src/bessely.jl
index a05f554..49e32fd 100644
--- a/src/bessely.jl
+++ b/src/bessely.jl
@@ -351,7 +351,7 @@ function bessely_power_series(v, x::T) where T
     out2 = zero(T)
     a = (x/2)^v
     # check for underflow and return limit for small arguments
-    iszero(a) && return -T(Inf)
+    iszero(a) && return (-T(Inf), a)
 
     b = inv(a)
     a /= gamma(v + one(T))

From 824b0cfa0c5e56fcae86613b5507e05ceb347505 Mon Sep 17 00:00:00 2001
From: Michael Helton <heltonmc@protonmail.com>
Date: Thu, 11 Aug 2022 15:20:55 -0400
Subject: [PATCH 3/6] support better float16, integer arguments

---
 src/besseli.jl | 16 ++++++++++++----
 src/besselj.jl | 16 +++++++++++++---
 src/besselk.jl | 10 +++++++---
 src/bessely.jl | 10 +++++++---
 4 files changed, 39 insertions(+), 13 deletions(-)

diff --git a/src/besseli.jl b/src/besseli.jl
index e17e23e..ba2e04a 100644
--- a/src/besseli.jl
+++ b/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::Real, x::T) where T
+    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
     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)
 
diff --git a/src/besselj.jl b/src/besselj.jl
index c31c80b..010bb05 100644
--- a/src/besselj.jl
+++ b/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::Real, 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
     abs_nu = abs(nu)
     abs_x = abs(x)
     sg = iseven(abs_nu) ? 1 : -1
diff --git a/src/besselk.jl b/src/besselk.jl
index 0d38b50..5f163ba 100644
--- a/src/besselk.jl
+++ b/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::Real, x::T) where T
+    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
     abs_nu = abs(nu)
     abs_x = abs(x)
     sg = iseven(abs_nu) ? 1 : -1
diff --git a/src/bessely.jl b/src/bessely.jl
index 49e32fd..f9e8ec0 100644
--- a/src/bessely.jl
+++ b/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::Real, 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
     abs_nu = abs(nu)
     abs_x = abs(x)
     sg = iseven(abs_nu) ? 1 : -1

From 7ac0361c3fb06799df74383c381256b53b39b5d0 Mon Sep 17 00:00:00 2001
From: Michael Helton <heltonmc@protonmail.com>
Date: Thu, 11 Aug 2022 15:52:08 -0400
Subject: [PATCH 4/6] fix type ambiguous

---
 src/besseli.jl               | 8 ++++----
 src/besselj.jl               | 6 +++---
 src/besselk.jl               | 6 +++---
 src/bessely.jl               | 6 +++---
 test/sphericalbessel_test.jl | 4 ++--
 5 files changed, 15 insertions(+), 15 deletions(-)

diff --git a/src/besseli.jl b/src/besseli.jl
index ba2e04a..7765f8b 100644
--- a/src/besseli.jl
+++ b/src/besseli.jl
@@ -166,9 +166,9 @@ Modified Bessel function of the second kind of order nu, ``I_{nu}(x)``.
 """
 besseli(nu::Real, x::Real) = _besseli(nu, float(x))
 
-besseli(nu, x::Float16) = Float16(_besseli(nu, Float32(x)))
+_besseli(nu, x::Float16) = Float16(_besseli(nu, Float32(x)))
 
-function _besseli(nu::Real, x::T) where T
+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)
@@ -191,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
@@ -231,7 +231,7 @@ Nu must be real.
 """
 besselix(nu::Real, x::Real) = _besselix(nu, float(x))
 
-besselix(nu, x::Float16) = Float16(_besselix(nu, Float32(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)
diff --git a/src/besselj.jl b/src/besselj.jl
index 010bb05..68b76f1 100644
--- a/src/besselj.jl
+++ b/src/besselj.jl
@@ -200,9 +200,9 @@ 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)))
+_besselj(nu, x::Float16) = Float16(_besselj(nu, Float32(x)))
 
-function _besselj(nu::Real, x::T) where T <: Union{Float32, Float64}
+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)
@@ -228,7 +228,7 @@ function _besselj(nu::Real, x::T) where T <: Union{Float32, Float64}
     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
diff --git a/src/besselk.jl b/src/besselk.jl
index 5f163ba..0335546 100644
--- a/src/besselk.jl
+++ b/src/besselk.jl
@@ -187,9 +187,9 @@ Modified Bessel function of the second kind of order nu, ``K_{nu}(x)``.
 """
 besselk(nu::Real, x::Real) = _besselk(nu, float(x))
 
-besselk(nu, x::Float16) = Float16(_besselk(nu, Float32(x)))
+_besselk(nu, x::Float16) = Float16(_besselk(nu, Float32(x)))
 
-function _besselk(nu::Real, x::T) where T
+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)
@@ -201,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
diff --git a/src/bessely.jl b/src/bessely.jl
index f9e8ec0..1d920aa 100644
--- a/src/bessely.jl
+++ b/src/bessely.jl
@@ -237,9 +237,9 @@ nu and x must be real where nu and x can be positive or negative.
 """
 bessely(nu::Real, x::Real) = _bessely(nu, float(x))
 
-bessely(nu, x::Float16) = Float16(_bessely(nu, Float32(x)))
+_bessely(nu, x::Float16) = Float16(_bessely(nu, Float32(x)))
 
-function _bessely(nu::Real, x::T) where T <: Union{Float32, Float64}
+function _bessely(nu, x::T) where T <: Union{Float32, Float64}
     isnan(nu) || isnan(x) && return NaN
     isinteger(nu) && return _bessely(Int(nu), x)
     abs_nu = abs(nu)
@@ -264,7 +264,7 @@ function _bessely(nu::Real, x::T) where T <: Union{Float32, Float64}
         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
diff --git a/test/sphericalbessel_test.jl b/test/sphericalbessel_test.jl
index 1aaa8c0..3606e7f 100644
--- a/test/sphericalbessel_test.jl
+++ b/test/sphericalbessel_test.jl
@@ -33,8 +33,8 @@ x = 1e-15
 @test isnan(Bessels.sphericalbessely(4.0, NaN))
 
 # test Float16 types
-@test Bessels.sphericalbesselj(1.4, Float16(1.2)) isa Float16
-@test Bessels.sphericalbessely(1.4, Float16(1.2)) isa Float16
+@test Bessels.sphericalbesselj(Float16(1.4), Float16(1.2)) isa Float16
+@test Bessels.sphericalbessely(Float16(1.4), Float16(1.2)) isa Float16
 
 for x in 0.5:1.0:100.0, v in [0, 1, 5.5, 8.2, 10]
     @test isapprox(Bessels.sphericalbesselj(v, x), SpecialFunctions.sphericalbesselj(v, x), rtol=1e-12)

From c095ae28ce4f903c86e087bc288d2bf82f2cd18d Mon Sep 17 00:00:00 2001
From: Michael Helton <heltonmc@protonmail.com>
Date: Thu, 11 Aug 2022 16:00:57 -0400
Subject: [PATCH 5/6] add news update f16 test

---
 NEWS.md              | 1 +
 test/besseli_test.jl | 2 ++
 test/besselj_test.jl | 3 +++
 test/besselk_test.jl | 3 +++
 test/bessely_test.jl | 3 +++
 5 files changed, 12 insertions(+)

diff --git a/NEWS.md b/NEWS.md
index 7f1c4d7..67033da 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -19,6 +19,7 @@ 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
 
diff --git a/test/besseli_test.jl b/test/besseli_test.jl
index 35eb704..6374008 100644
--- a/test/besseli_test.jl
+++ b/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]
diff --git a/test/besselj_test.jl b/test/besselj_test.jl
index 07abd36..e679bb0 100644
--- a/test/besselj_test.jl
+++ b/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)
 
diff --git a/test/besselk_test.jl b/test/besselk_test.jl
index ab5377d..aaa8de6 100644
--- a/test/besselk_test.jl
+++ b/test/besselk_test.jl
@@ -120,6 +120,9 @@ 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 Inf
 @test iszero(besselk(2, Inf))
 
diff --git a/test/bessely_test.jl b/test/bessely_test.jl
index e563e44..a737c3c 100644
--- a/test/bessely_test.jl
+++ b/test/bessely_test.jl
@@ -89,6 +89,9 @@ 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

From df2db0c74d41c8fa57dc24012124db043fd9a033 Mon Sep 17 00:00:00 2001
From: Michael Helton <heltonmc@protonmail.com>
Date: Thu, 11 Aug 2022 16:10:50 -0400
Subject: [PATCH 6/6] fix besselk f16

---
 src/besselk.jl       | 6 +++++-
 test/besseli_test.jl | 1 +
 test/besselk_test.jl | 1 +
 3 files changed, 7 insertions(+), 1 deletion(-)

diff --git a/src/besselk.jl b/src/besselk.jl
index 0335546..ba1382f 100644
--- a/src/besselk.jl
+++ b/src/besselk.jl
@@ -243,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)
 
diff --git a/test/besseli_test.jl b/test/besseli_test.jl
index 6374008..e20dfe9 100644
--- a/test/besseli_test.jl
+++ b/test/besseli_test.jl
@@ -92,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
 
diff --git a/test/besselk_test.jl b/test/besselk_test.jl
index aaa8de6..cc8b037 100644
--- a/test/besselk_test.jl
+++ b/test/besselk_test.jl
@@ -122,6 +122,7 @@ 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))