update besselk, besseli

Author: heltonmc

src/besseli.jl

@@ -217,7 +217,7 @@ function besseli_positive_args(nu, x::T) where T <: Union{Float32, Float64}
     besseli_large_argument_cutoff(nu, x) && return besseli_large_argument(nu, x)
 
     # use uniform debye expansion if x or nu is large
-    besselik_debye_cutoff(nu, x) && return T(besseli_large_orders(nu, x))
+    besselik_debye_cutoff(nu, x) && return besseli_large_orders(nu, x)
 
     # for rest of values use the power series
     return besseli_power_series(nu, x)
@@ -236,12 +236,13 @@ _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)
+    isinf(x) && return T(Inf)
 
     # use large argument expansion if x >> nu
     besseli_large_argument_cutoff(nu, x) && return besseli_large_argument_scaled(nu, x)
 
     # use uniform debye expansion if x or nu is large
-    besselik_debye_cutoff(nu, x) && return T(besseli_large_orders_scaled(nu, x))
+    besselik_debye_cutoff(nu, x) && return besseli_large_orders_scaled(nu, x)
 
     # for rest of values use the power series
     return besseli_power_series(nu, x) * exp(-x)
@@ -259,7 +260,7 @@ end
 
 Debey's uniform asymptotic expansion for large order valid when v-> ∞ or x -> ∞
 """
-function besseli_large_orders(v, x::T) where T <: Union{Float32, Float64}
+function besseli_large_orders(v, x::T) where T
     S = promote_type(T, Float64)
     x = S(x)
     z = x / v
@@ -269,10 +270,10 @@ function besseli_large_orders(v, x::T) where T <: Union{Float32, Float64}
     p = inv(zs)
     p2  = v^2/fma(max(v,x), max(v,x), min(v,x)^2)
 
-    return coef*Uk_poly_In(p, v, p2, T)
+    return T(coef*Uk_poly_In(p, v, p2, T))
 end
 
-function besseli_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64}
+function besseli_large_orders_scaled(v, x::T) where T
     S = promote_type(T, Float64)
     x = S(x)
     z = x / v
@@ -296,42 +297,35 @@ end
 
 Debey's uniform asymptotic expansion for large order valid when v-> ∞ or x -> ∞
 """
-function besseli_large_argument(v, z::T) where T
-    MaxIter = 1000
-    a = exp(z / 2)
-    coef = a / sqrt(2 * T(pi) * z)
-    fv2 = 4 * v^2
-    term = one(T)
-    res = term
-    s = -term
-    for i in 1:MaxIter
-        i = T(i)
-        offset = muladd(2, i, -1)
-        term *= T(0.125) * muladd(offset, -offset, fv2) / (z * i)
-        res = muladd(term, s, res)
-        s = -s
-        abs(term) <= eps(T) && break
-    end
-    return res * coef * a
+function besseli_large_argument(v, x::T) where T
+    a = exp(x / 2)
+    coef = a / sqrt(2 * (π * x))
+    return T(_besseli_large_argument(v, x) * coef * a)
 end
-function besseli_large_argument_scaled(v, z::T) where T
-    MaxIter = 1000
-    coef = inv(sqrt(2 * T(pi) * z))
+
+besseli_large_argument_scaled(v, x::T) where T =  T(_besseli_large_argument(v, x) / sqrt(2 * (π * x)))
+
+function _besseli_large_argument(v, x::T) where T
+    MaxIter = 5000
+    S = promote_type(T, Float64)
+    v, x = S(v), S(x)
+
     fv2 = 4 * v^2
-    term = one(T)
+    term = one(S)
     res = term
     s = -term
     for i in 1:MaxIter
-        i = T(i)
         offset = muladd(2, i, -1)
-        term *= T(0.125) * muladd(offset, -offset, fv2) / (z * i)
+        term *= muladd(offset, -offset, fv2) / (8 * x * i)
         res = muladd(term, s, res)
         s = -s
         abs(term) <= eps(T) && break
     end
-    return res * coef
+    return res
 end
-besseli_large_argument_cutoff(nu, x) = x > maximum((30.0, nu^2 / 6))
+
+besseli_large_argument_cutoff(nu, x::Float64) = x > 23.0 && x > nu^2 / 1.8 + 23.0
+besseli_large_argument_cutoff(nu, x::Float32) = x > 18.0f0 && x > nu^2 / 19.5f0 + 18.0f0
 
 #####
 #####  Power series for I_{nu}(x)

src/besselk.jl

@@ -301,7 +301,7 @@ function besselk_large_orders_scaled(v, x::T) where T
     return T(coef*Uk_poly_Kn(p, v, p2, T))
 end
 besselik_debye_cutoff(nu, x::Float64) = nu > 25.0 || x > 35.0
-besselik_debye_cutoff(nu, x::Float32) = nu > 15.0 || x > 20.0
+besselik_debye_cutoff(nu, x::Float32) = nu > 15.0f0 || x > 20.0f0
 
 #####
 #####  Continued fraction with Wronskian for K_{nu}(x)
@@ -389,35 +389,38 @@ No checks are performed on nu so this is not accurate when nu is an integer.
 """
 function besselk_power_series(v, x::T) where T
     MaxIter = 1000
+    S = promote_type(T, Float64)
+    v, x = S(v), S(x)
+
     z  = x / 2
     zz = z * z
-    
     logz = log(z)
     xd2_v = exp(v*logz)
     xd2_nv = inv(xd2_v)
 
     # use the reflection identify to calculate gamma(-v)
     # use relation gamma(v)*v = gamma(v+1) to avoid two gamma calls
     gam_v = gamma(v)
-    xp1 = abs(v) + one(T)
+    xp1 = abs(v) + one(S)
     gam_nv = π / (sinpi(xp1) * gam_v * v)
     gam_1mv = -gam_nv * v
     gam_1mnv = gam_v * v
 
     _t1 = gam_v * xd2_nv * gam_1mv
     _t2 = gam_nv * xd2_v * gam_1mnv
-    (xd2_pow, fact_k, out) = (one(T), one(T), zero(T))
+    (xd2_pow, fact_k, out) = (one(S), one(S), zero(S))
     for k in 0:MaxIter
         t1 = xd2_pow * T(0.5)
         tmp = muladd(_t1, gam_1mnv, _t2 * gam_1mv)
         tmp *= inv(gam_1mv * gam_1mnv * fact_k)
         term = t1 * tmp
         out += term
         abs(term / out) < eps(T) && break
-        (gam_1mnv, gam_1mv) = (gam_1mnv*(one(T) + v + k), gam_1mv*(one(T) - v + k)) 
+        (gam_1mnv, gam_1mv) = (gam_1mnv*(one(S) + v + k), gam_1mv*(one(S) - v + k)) 
         xd2_pow *= zz
-        fact_k *= k + one(T)
+        fact_k *= k + one(S)
     end
-    return out
+    return T(out)
 end
-besselk_power_series_cutoff(nu, x) = x < 2.0 || nu > 1.6x - 1.0
+besselk_power_series_cutoff(nu, x::Float64) = x < 2.0 || nu > 1.6x - 1.0
+besselk_power_series_cutoff(nu, x::Float32) = x < 10.0f0 || nu > 1.65f0*x - 8.0f0

src/bessely.jl

@@ -34,7 +34,11 @@
 
 Bessel function of the second kind of order zero, ``Y_0(x)``.
 """
-function bessely0(x::T) where T <: Union{Float32, Float64}
+bessely0(x::Real) = _bessely0(float(x))
+
+_bessely0(x::Float16) = Float16(_bessely0(Float32(x)))
+
+function _bessely0(x::T) where T <: Union{Float32, Float64}
     if x <= zero(x)
         if iszero(x)
             return T(-Inf)
@@ -110,7 +114,11 @@ end
 
 Bessel function of the second kind of order one, ``Y_1(x)``.
 """
-function bessely1(x::T) where T <: Union{Float32, Float64}
+bessely1(x::Real) = _bessely1(float(x))
+
+_bessely1(x::Float16) = Float16(_bessely1(Float32(x)))
+
+function _bessely1(x::T) where T <: Union{Float32, Float64}
     if x <= zero(x)
         if iszero(x)
             return T(-Inf)
@@ -318,6 +326,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]
 
+    # shift nu down and use upward recurrence starting from either chebyshev approx or hankel expansion
     return bessely_fallback(nu, x)
 end
 
@@ -369,7 +378,7 @@ function bessely_power_series(v, x::T) where T
     return (out*c - out2) / s, 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.0) || v > 1.38*x - 12.5
+bessely_series_cutoff(v, x::Float32) = (x < 21.0f0) || v > 1.38f0*x - 12.5f0
 
 #####
 #####  Fallback for Y_{nu}(x)
@@ -382,7 +391,7 @@ function bessely_fallback(nu, x)
     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
+        nu_shift = ceil(Int, 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

test/besseli_test.jl

@@ -103,6 +103,7 @@ for v in nu, xx in x
     xx *= v
     sf = SpecialFunctions.besseli(v, xx)
     @test isapprox(besseli(v, xx), Float64(sf), rtol=2e-13)
+    @test isapprox(besseli(Float32(v), Float32(xx)), Float32(sf))
 end
 
 # test Inf

test/besselk_test.jl

@@ -108,7 +108,6 @@ t = [besselk(m, x) for m in m, x in x]
 @test besselkx(3.2, 1.1) ≈ SpecialFunctions.besselk(3.2, 1.1)*exp(1.1)
 @test besselkx(8.2, 9.1) ≈ SpecialFunctions.besselk(8.2, 9.1)*exp(9.1)
 
-
 ## Tests for besselk
 
 ## test all numbers and orders for 0<nu<100
@@ -118,6 +117,7 @@ for v in nu, xx in x
     xx *= v
     sf = SpecialFunctions.besselk(v, xx)
     @test isapprox(besselk(v, xx), Float64(sf), rtol=2e-13)
+    @test isapprox(besselk(Float32(v), Float32(xx)), Float32(sf))
 end
 
 # test Float16