stabilize Float32 functions

Author: oscardssmith

src/BesselFunctions/U_polynomials.jl

@@ -40,7 +40,7 @@ function besseljy_debye(v, x::T) where T
     p2  = v^2 / vmx
 
     Uk_Jn, Uk_Yn = Uk_poly_Jn(p, v, p2, T(x))
-    return coef_Jn * Uk_Jn, coef_Yn * Uk_Yn
+    return T(coef_Jn * Uk_Jn), T(coef_Yn * Uk_Yn)
 end
 
 # Cutoffs for besseljy_debye expansions
@@ -102,7 +102,7 @@ function hankel_debye(v, x::T) where T
     p2  = v^2 / vmx
 
     _, Uk_Yn = Uk_poly_Hankel(p*im, v, -p2, T(x))
-    return coef_Yn * Uk_Yn
+    return T(coef_Yn * Uk_Yn)
 end
 
 # Cutoffs for hankel_debye expansions

src/BesselFunctions/asymptotics.jl

@@ -51,7 +51,7 @@ function besseljy_large_argument(v, x::T) where T
     s3 = CMS * sα
     s4 = CPS * cα
 
-    return SQ2O2(S) * (s1 + s2) * b, SQ2O2(S) * (s3 - s4) * b
+    return T(SQ2O2(S) * (s1 + s2) * b), T(SQ2O2(S) * (s3 - s4) * b)
 end
 
 # Float64
@@ -101,11 +101,11 @@ end
 function _α_αp_asymptotic(v, x::Float32)
     v, x = Float64(v), Float64(x)
     if x > 4*v
-        return _α_αp_poly_5(v, x)
+        return Float32.(_α_αp_poly_5(v, x))
     elseif x > 1.8*v
-        return _α_αp_poly_10(v, x)
+        return Float32.(_α_αp_poly_10(v, x))
     else
-        return _α_αp_poly_30(v, x)
+        return Float32.(_α_αp_poly_30(v, x))
     end
 end
 

src/BesselFunctions/bessely.jl

@@ -420,7 +420,7 @@ function bessely_power_series(v, x::T) where T
     out2 = zero(S)
     a = (x/2)^v
     # check for underflow and return limit for small arguments
-    iszero(a) && return (-T(Inf), a)
+    iszero(a) && return (-T(Inf), T(a))
 
     b = inv(a)
     a /= gamma(v + one(S))
@@ -434,7 +434,7 @@ function bessely_power_series(v, x::T) where T
         b *= -inv((-v + i + one(S)) * (i + one(S))) * t2
     end
     s, c = sincospi(v)
-    return (out*c - out2) / s, out
+    return T((out*c - out2) / s), T(out)
 end
 bessely_series_cutoff(v, x::Float64) = (x < 7.0) || v > 1.35*x - 4.5
 bessely_series_cutoff(v, x::Float32) = (x < 21.0f0) || v > 1.38f0*x - 12.5f0
@@ -484,7 +484,7 @@ function bessely_chebyshev_low_orders(v, x)
     x1 = 2*(x - 6)/13 - 1
     v1 = v - 1
     v2 = v
-    a = clenshaw_chebyshev.(x1, bessely_cheb_weights)
+    a = clenshaw_chebyshev.(x1, map(Base.Fix1(map, typeof(x1)),bessely_cheb_weights))
     return clenshaw_chebyshev(v1, a), clenshaw_chebyshev(v2, a)
 end
 

src/BesselFunctions/hankel.jl

@@ -133,7 +133,7 @@ function besseljy_positive_args(nu::Real, x::T) where T
 
     # 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
+    nu_shift = ceil(nu) - floor(Int, -3//2 + x + Base.Math._approx_cbrt(-411*x)) + 2
     v2 = maximum((nu - nu_shift, modf(nu)[1] + 1))
 
     Hnu = hankel_debye(v2, x)