add tests, fix F32

Author: heltonmc

src/besseli.jl

@@ -137,37 +137,51 @@ Modified Bessel function of the first kind of order nu, ``I_{nu}(x)``.
 Nu must be real.
 """
 function besseli(nu, x::T) where T <: Union{Float32, Float64}
-    nu == 0 && return besseli0(x)
-    nu == 1 && return besseli1(x)
+    iszero(nu) && return besseli0(x)
+    isone(nu) && return besseli1(x)
 
-    if x > maximum((T(30), nu^2 / 6))
-        return T(besseli_large_argument(nu, x))
-    elseif nu > 25.0 || x > 35.0
-        return T(besseli_large_orders(nu, x))
-    else
-        return T(besseli_power_series(nu, x))
-    end
+    # use large argument expansion if x >> nu
+    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))
+
+    # for rest of values use the power series
+    return besseli_power_series(nu, x)
 end
+
 """
     besselix(nu, x::T) where T <: Union{Float32, Float64}
 
 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}
-    nu == 0 && return besseli0x(x)
-    nu == 1 && return besseli1x(x)
-
-    if x > maximum((T(30), nu^2 / 4))
-        return T(besseli_large_argument_scaled(nu, x))
-    elseif x <= 2 * sqrt(nu + 1)
-        return T(besseli_power_series(nu, x)) * exp(-x)
-    elseif nu < 100
-        return T(_besseli_continued_fractions_scaled(nu, x))
-    else
-        return besseli_large_orders_scaled(nu, x)
-    end
+    iszero(nu) && return besseli0x(x)
+    isone(nu) && return besseli1x(x)
+
+    # 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))
+
+    # for rest of values use the power series
+    return besseli_power_series(nu, x) * exp(-x)
 end
+
+#####
+#####  Debye's uniform asymptotic for I_{nu}(x)
+#####
+
+# Implements the uniform asymptotic expansion https://dlmf.nist.gov/10.41
+# In general this is valid when either x or nu is gets large
+# see the file src/U_polynomials.jl for more details
+"""
+    besseli_large_orders(nu, x::T)
+
+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}
     S = promote_type(T, Float64)
     x = S(x)
@@ -180,6 +194,7 @@ function besseli_large_orders(v, x::T) where T <: Union{Float32, Float64}
 
     return coef*Uk_poly_In(p, v, p2, T)
 end
+
 function besseli_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64}
     S = promote_type(T, Float64)
     x = S(x)
@@ -192,39 +207,18 @@ function besseli_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64}
 
     return T(coef*Uk_poly_In(p, v, p2, T))
 end
-function _besseli_continued_fractions(nu, x::T) where T
-    S = promote_type(T, Float64)
-    xx = S(x)
-    knum1, knu = besselk_up_recurrence(xx, besselk1(xx), besselk0(xx), 1, nu-1)
-    # if knu or knum1 is zero then besseli will likely overflow
-    (iszero(knu) || iszero(knum1)) && return throw(DomainError(x, "Overflow error"))
-    return 1 / (x * (knum1 + knu / steed(nu, x)))
-end
-function _besseli_continued_fractions_scaled(nu, x::T) where T
-    S = promote_type(T, Float64)
-    xx = S(x)
-    knum1, knu = besselk_up_recurrence(xx, besselk1x(xx), besselk0x(xx), 1, nu-1)
-    # if knu or knum1 is zero then besseli will likely overflow
-    (iszero(knu) || iszero(knum1)) && return throw(DomainError(x, "Overflow error"))
-    return 1 / (x * (knum1 + knu / steed(nu, x)))
-end
-function steed(n, x::T) where T
-    MaxIter = 1000
-    xinv = inv(x)
-    xinv2 = 2 * xinv
-    d = x / (n + n)
-    a = d
-    h = a
-    b = muladd(2, n, 2) * xinv
-    for _ in 1:MaxIter
-        d = inv(b + d)
-        a *= muladd(b, d, -1)
-        h = h + a
-        b = b + xinv2
-        abs(a / h) <= eps(T) && break
-    end
-    return h
-end
+
+#####
+#####  Large argument expansion (x>>nu) for I_{nu}(x)
+#####
+
+# Implements the uniform asymptotic expansion https://dlmf.nist.gov/10.40.E1
+# In general this is valid when x > nu^2
+"""
+    besseli_large_orders(nu, x::T)
+
+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)
@@ -260,7 +254,19 @@ function besseli_large_argument_scaled(v, z::T) where T
     end
     return res * coef
 end
+besseli_large_argument_cutoff(nu, x) = x > maximum((30.0, nu^2 / 6))
+
+#####
+#####  Power series for I_{nu}(x)
+#####
+
+# Use power series form of I_v(x) which is generally accurate across all values though slower for larger x
+# https://dlmf.nist.gov/10.25.E2
+"""
+    besseli_power_series(nu, x::T) where T <: Float64
 
+Computes ``I_{nu}(x)`` using the power series for any value of nu.
+"""
 function besseli_power_series(v, x::T) where T
     MaxIter = 3000
     out = zero(T)
@@ -274,3 +280,45 @@ function besseli_power_series(v, x::T) where T
     end
     return out
 end
+
+#=
+# the following is a deprecated version of the continued fraction approach
+# using K0 and K1 as starting values then forward recurrence up till nu
+# then using the wronskian to getting I_{nu}
+# in general this method is slow and depends on starting values of K0 and K1
+# which is not very flexible for arbitary orders
+
+function _besseli_continued_fractions(nu, x::T) where T
+    S = promote_type(T, Float64)
+    xx = S(x)
+    knum1, knu = besselk_up_recurrence(xx, besselk1(xx), besselk0(xx), 1, nu-1)
+    # if knu or knum1 is zero then besseli will likely overflow
+    (iszero(knu) || iszero(knum1)) && return throw(DomainError(x, "Overflow error"))
+    return 1 / (x * (knum1 + knu / steed(nu, x)))
+end
+function _besseli_continued_fractions_scaled(nu, x::T) where T
+    S = promote_type(T, Float64)
+    xx = S(x)
+    knum1, knu = besselk_up_recurrence(xx, besselk1x(xx), besselk0x(xx), 1, nu-1)
+    # if knu or knum1 is zero then besseli will likely overflow
+    (iszero(knu) || iszero(knum1)) && return throw(DomainError(x, "Overflow error"))
+    return 1 / (x * (knum1 + knu / steed(nu, x)))
+end
+function steed(n, x::T) where T
+    MaxIter = 1000
+    xinv = inv(x)
+    xinv2 = 2 * xinv
+    d = x / (n + n)
+    a = d
+    h = a
+    b = muladd(2, n, 2) * xinv
+    for _ in 1:MaxIter
+        d = inv(b + d)
+        a *= muladd(b, d, -1)
+        h = h + a
+        b = b + xinv2
+        abs(a / h) <= eps(T) && break
+    end
+    return h
+end
+=#

src/gamma.jl

@@ -8,6 +8,8 @@ function gamma(x)
         return _gamma(x)
     end
 end
+# only have a Float64 implementations
+gamma(x::Float32) = Float32(gamma(Float64(x)))
 function _gamma(x)
     if x > 11.5
         return large_gamma(x)

test/besseli_test.jl

@@ -90,4 +90,16 @@ 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]
+
+## Tests for besselk
+
+## test all numbers and orders for 0<nu<100
+x = [0.01, 0.05, 0.1, 0.2, 0.4, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.91, 0.92, 0.93, 0.95, 0.96, 0.97, 0.98, 0.99, 0.995, 0.999, 1.0, 1.001, 1.01, 1.05, 1.1, 1.2, 1.4, 1.6, 1.8, 1.9, 2.5, 3.0, 3.5, 4.0]
+nu = [0.01,0.1, 0.5, 0.8, 1, 1.23, 2,2.56, 4,5.23, 6,9.2, 10,12.89, 15, 19.1, 20, 25, 30, 33.123, 40, 45, 50, 51.5, 55, 60, 65, 70, 72.34, 75, 80, 82.1, 85, 88.76, 90, 92.334, 95, 99.87,100, 110, 125, 145.123, 150, 160.789]
+for v in nu, xx in x
+    xx *= v
+    sf = SpecialFunctions.besseli(v, xx)
+    @test isapprox(besseli(v, xx), Float64(sf), rtol=2e-13)
+end
+