add more generics

Author: heltonmc

src/sphericalbessel.jl

@@ -22,37 +22,44 @@
 Spherical bessel function of the first kind of order `nu`, ``j_ν(x)``. This is the non-singular
 solution to the radial part of the Helmholz equation in spherical coordinates.
 """
-function sphericalbesselj(nu::Real, x::T) where T
+sphericalbesselj(nu::Real, x::Real) = _sphericalbesselj(nu, float(x))
+
+_sphericalbesselj(nu, x::Float16) = Float16(_sphericalbesselj(nu, Float32(x)))
+    
+function _sphericalbesselj(nu::Real, x::T) where T <: Union{Float32, Float64}
     isnan(nu) || isnan(x) && return NaN
     x < zero(T) && return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
 
-    if nu < zero(T)
-        return SQPIO2(T) * besselj(nu + 1/2, x) / sqrt(x)
-    else
-        return sphericalbesselj_positive_args(nu, x)
-    end
+    return nu < zero(T) ? sphericalbesselj_generic(nu, x) : sphericalbesselj_positive_args(nu, x)
 end
 
+sphericalbesselj_generic(nu, x::T) where T = SQPIO2(T) * besselj(nu + one(T)/2, x) / sqrt(x)
+
 #####
 ##### Positive arguments for `sphericalbesselj`
 #####
 
 function sphericalbesselj_positive_args(nu::Real, x::T) where T
-    if x^2 / (4*nu + 110) < eps(T)
-        # small arguments power series expansion
-        x2 = x^2 / 4
-        coef = evalpoly(x2, (1, -inv(3/2 + nu), -inv(5 + nu), -inv(21/2 + nu), -inv(18 + nu)))
-        a = sqrt(T(pi)/2) / (gamma(T(3)/2 + nu) * 2^(nu + T(1)/2))
-        return x^nu * a * coef
-    elseif isinteger(nu)
-        if (x >= nu && nu < 250) || (x < nu && nu < 60)
-            return sphericalbesselj_recurrence(Int(nu), x)
-        else
-            return SQPIO2(T) * besselj(nu + 1/2, x) / sqrt(x)
-        end
-    else
-        return SQPIO2(T) * besselj(nu + 1/2, x) / sqrt(x)
-    end
+    isinteger(nu) && return sphericalbesselj_positive_args(Int(nu), x)
+    return sphericalbesselj_small_args_cutoff(nu, x) ? sphericalbesselj_small_args(nu, x) : sphericalbesselj_generic(nu, x)
+end
+
+function sphericalbesselj_positive_args(nu::Integer, x::T) where T
+    sphericalbesselj_small_args_cutoff(nu, x) && return sphericalbesselj_small_args(nu, x)
+    return (x >= nu && nu < 250) || (x < nu && nu < 60) ? sphericalbesselj_recurrence(nu, x) : sphericalbesselj_generic(nu, x)
+end
+
+#####
+##### Power series expansion for small arguments
+#####
+
+sphericalbesselj_small_args_cutoff(nu, x::T) where T = x^2 / (4*nu + 110) < eps(T)
+
+function sphericalbesselj_small_args(nu, x::T) where T
+    x2 = x^2 / 4
+    coef = evalpoly(x2, (1, -inv(T(3)/2 + nu), -inv(5 + nu), -inv(T(21)/2 + nu), -inv(18 + nu)))
+    a = sqrt(T(pi)/2) / (gamma(T(3)/2 + nu) * 2^(nu + T(1)/2))
+    return x^nu * a * coef
 end
 
 #####
@@ -62,7 +69,7 @@ end
 # very accurate approach however need to consider some performance issues
 # if recurrence is stable (x>=nu) can generate very fast up to orders around 250
 # for larger orders it is more efficient to use other expansions
-# if (x<nu) we can use forward recurrence from sphericalbesselj_recurrence and
+# if (x<nu) we can use forward recurrence from sphericalbessely_recurrence and
 # then use a continued fraction approach. However, for largish orders (>60) the
 # continued fraction is slower converging and more efficient to use other methods
 function sphericalbesselj_recurrence(nu::Integer, x::T) where T
@@ -98,30 +105,32 @@ Spherical bessel function of the second kind at order `nu`, ``y_ν(x)``. This is
 solution to the radial part of the Helmholz equation in spherical coordinates. Sometimes
 known as a spherical Neumann function.
 """
-function sphericalbessely(nu::Real, x::T) where T
+sphericalbessely(nu::Real, x::Real) = _sphericalbessely(nu, float(x))
+
+_sphericalbessely(nu, x::Float16) = Float16(_sphericalbessely(nu, Float32(x)))
+
+function _sphericalbessely(nu::Real, x::T) where T <: Union{Float32, Float64}
     isnan(nu) || isnan(x) && return NaN
     x < zero(T) && return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
-    abs_nu = abs(nu)
-
-    if nu < zero(T)
-        return SQPIO2(T) * bessely(nu + 1/2, x) / sqrt(x)
-    else
-        return sphericalbessely_positive_args(nu, x)
-    end
+    return nu < zero(T) ? sphericalbessely_generic(nu, x) : sphericalbessely_positive_args(nu, x)
 end
 
+sphericalbessely_generic(nu, x::T) where T = SQPIO2(T) * bessely(nu + one(T)/2, x) / sqrt(x)
+
 #####
 ##### Positive arguments for `sphericalbesselj`
 #####
 
-function sphericalbessely_positive_args(nu::Real, x::T) where T
+sphericalbessely_positive_args(nu::Real, x) = isinteger(nu) ? sphericalbessely_positive_args(Int(nu), x) : sphericalbessely_generic(nu, x)
+
+function sphericalbessely_positive_args(nu::Integer, x::T) where T
     if besseljy_debye_cutoff(nu, x)
         # for very large orders use expansion nu >> x to avoid overflow in recurrence
         return SQPIO2(T) * besseljy_debye(nu + 1/2, x)[2] / sqrt(x)
-    elseif isinteger(nu) && nu < 250
-        return sphericalbessely_forward_recurrence(Int(nu), x)[1]
+    elseif nu < 250
+        return sphericalbessely_forward_recurrence(nu, x)[1]
     else
-        return SQPIO2(T) * bessely(nu + 1/2, x) / sqrt(x)
+        return sphericalbessely_generic(nu, x)
     end
 end
 

test/sphericalbessel_test.jl

@@ -9,7 +9,7 @@ x = 1e-15
 @test Bessels.sphericalbessely(5.5, x) ≈ SpecialFunctions.sphericalbessely(5.5, x)
 @test Bessels.sphericalbessely(10, x) ≈ SpecialFunctions.sphericalbessely(10, x)
 
-for x in 0.5:1.0:100.0, v in [0, 1, 5.5, 10]
+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)
     @test isapprox(Bessels.sphericalbessely(v, x), SpecialFunctions.sphericalbessely(v, x), rtol=1e-12)
 end