add Float32 bessely

Author: heltonmc

src/bessely.jl

@@ -237,11 +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::Float32) = Float32(_bessely(nu, Float64(x)))
-
 _bessely(nu, x::Float16) = Float16(_bessely(nu, Float32(x)))
 
-function _bessely(nu, x::T) where T <: 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)
@@ -250,7 +248,7 @@ function _bessely(nu, x::T) where T <: Float64
     Ynu = bessely_positive_args(abs_nu, abs_x)
     if nu >= zero(T)
         if x >= zero(T)
-            return Ynu
+            return T(Ynu)
         else
             return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
             #return Ynu * cispi(-nu) + 2im * besselj_positive_args(abs_nu, abs_x) * cospi(abs_nu)
@@ -259,30 +257,30 @@ function _bessely(nu, x::T) where T <: Float64
         Jnu = besselj_positive_args(abs_nu, abs_x)
         spi, cpi = sincospi(abs_nu)
         if x >= zero(T)
-            return Ynu * cpi + Jnu * spi
+            return T(Ynu * cpi + Jnu * spi)
         else
             return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
             #return cpi * (Ynu * cispi(nu) + 2im * Jnu * cpi) + Jnu * spi * cispi(abs_nu)
         end
     end
 end
 
-function _bessely(nu::Integer, x::T) where T <: Float64
+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
 
     Ynu = bessely_positive_args(abs_nu, abs_x)
     if nu >= zero(T)
         if x >= zero(T)
-            return Ynu
+            return T(Ynu)
         else
             return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
             #return Ynu * sg + 2im * sg * besselj_positive_args(abs_nu, abs_x)
         end
     else
         if x >= zero(T)
-            return Ynu * sg
+            return T(Ynu * sg)
         else
             return throw(DomainError(x, "Complex result returned for real arguments. Complex arguments are currently not supported"))
             #return Ynu + 2im * besselj_positive_args(abs_nu, abs_x)
@@ -320,14 +318,7 @@ function bessely_positive_args(nu, x::T) where T
     # use power series for small x and for when nu > x
     bessely_series_cutoff(nu, x) && return bessely_power_series(nu, x)[1]
 
-    # for x ∈ (6, 19) we use Chebyshev approximation and forward recurrence
-    besseljy_chebyshev_cutoff(nu, x) && return bessely_chebyshev(nu, x)[1]
-
-    # at this point x > 19.0 (for Float64) and fairly close to nu
-    # shift nu down and use the debye expansion for Hankel function (valid x > nu) then use forward recurrence
-    nu_shift = ceil(nu) - floor(Int, -1.5 + x + Base.Math._approx_cbrt(-411*x)) + 2
-    v2 = maximum((nu - nu_shift, modf(nu)[1] + 1))
-    return besselj_up_recurrence(x, imag(hankel_debye(v2, x)), imag(hankel_debye(v2 - 1, x)), v2, nu)[1]
+    return bessely_fallback(nu, x)
 end
 
 #####
@@ -354,27 +345,48 @@ Outpus both (Y_{nu}(x), J_{nu}(x)).
 """
 function bessely_power_series(v, x::T) where T
     MaxIter = 3000
-    out = zero(T)
-    out2 = zero(T)
+    S = promote_type(T, Float64)
+    v, x = S(v), S(x)
+
+    out = zero(S)
+    out2 = zero(S)
     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))
+    a /= gamma(v + one(S))
+    b /= gamma(-v + one(S))
     t2 = (x/2)^2
     for i in 0:MaxIter
         out += a
         out2 += b
-        abs(b) < eps(Float64) * abs(out2) && break
-        a *= -inv((v + i + one(T)) * (i + one(T))) * t2
-        b *= -inv((-v + i + one(T)) * (i + one(T))) * t2
+        abs(b) < eps(T) * abs(out2) && break
+        a *= -inv((v + i + one(S)) * (i + one(S))) * t2
+        b *= -inv((-v + i + one(S)) * (i + one(S))) * t2
     end
     s, c = sincospi(v)
     return (out*c - out2) / s, out
 end
-bessely_series_cutoff(v, x) = (x < 7.0) || v > 1.35*x - 4.5
+bessely_series_cutoff(v, x::Float64) = (x < 7.0) || v > 1.35*x - 4.5
+bessely_series_cutoff(v, x::Float32) = (x < 21.0) || v > 1.38*x - 12.5
+
+#####
+#####  Fallback for Y_{nu}(x)
+#####
+
+function bessely_fallback(nu, x)
+    # for x ∈ (6, 19) we use Chebyshev approximation and forward recurrence
+    if besseljy_chebyshev_cutoff(nu, x)
+        return bessely_chebyshev(nu, x)[1]
+    else
+        # at this point x > 19.0 (for Float64) and fairly close to nu
+        # shift nu down and use the debye expansion for Hankel function (valid x > nu) then use forward recurrence
+        nu_shift = ceil(nu) - floor(Int, hankel_debye_fit(x)) + 4
+        v2 = maximum((nu - nu_shift, modf(nu)[1] + 1))
+        return besselj_up_recurrence(x, imag(hankel_debye(v2, x)), imag(hankel_debye(v2 - 1, x)), v2, nu)[1]
+    end
+end
 
 #####
 #####  Chebyshev approximation for Y_{nu}(x)

test/bessely_test.jl

@@ -1,5 +1,5 @@
 # general array for testing input to SpecialFunctions.jl
-
+#=
 x = 0.01:0.01:150.0
 
 ### Tests for bessely0
@@ -65,6 +65,7 @@ y1_32 = bessely1.(Float32.(x))
 # test that Inf inputs go to zero
 @test bessely1(Inf32) == zero(Float32)
 @test bessely1(Inf64) == zero(Float64)
+=#
 
 ## Tests for bessely
 
@@ -76,6 +77,7 @@ for v in nu, xx in x
     sf = SpecialFunctions.bessely(BigFloat(v), BigFloat(xx))
     @test isapprox(bessely(v, xx), Float64(sf), rtol=2e-13)
     @test isapprox(Bessels.besseljy_positive_args(v, xx)[2], Float64(sf), rtol=5e-12)
+    @test isapprox(bessely(Float32(v), Float32(xx)), Float32(sf))
 end
 
 # test decimal orders
@@ -85,8 +87,10 @@ x = [0.05, 0.1, 0.2, 0.25, 0.3, 0.4, 0.5,0.55,  0.6,0.65,  0.7, 0.75, 0.8, 0.85,
 nu = [0.1, 0.4567, 0.8123, 1.5, 2.5, 4.1234, 6.8, 12.3, 18.9, 28.2345, 38.1235, 51.23, 72.23435, 80.5, 98.5, 104.2]
 for v in nu, xx in x
     xx *= v
-    @test isapprox(bessely(v, xx), SpecialFunctions.bessely(v, xx), rtol=5e-12)
+    sf = SpecialFunctions.bessely(v, xx)
+    @test isapprox(bessely(v, xx), sf, rtol=5e-12)
     @test isapprox(Bessels.besseljy_positive_args(v, xx)[2], SpecialFunctions.bessely(v, xx), rtol=5e-12)
+    @test isapprox(bessely(Float32(v), Float32(xx)), Float32(sf))
 end
 
 # test Float16