update continued fractions

Author: heltonmc

src/besseli.jl

@@ -161,11 +161,12 @@ function besselix(nu, x::T) where T <: Union{Float32, Float64}
     nu == 0 && return besseli0x(x)
     nu == 1 && return besseli1x(x)
 
-    branch = 60
-    if nu < branch
-        inup1 = besseli_large_orders_scaled(branch + 1, x)
-        inu = besseli_large_orders_scaled(branch, x)
-        return down_recurrence(x, inu, inup1, nu, branch)
+    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_small_arguments(nu, x)) * exp(-x)
+    elseif nu < 100
+        return T(_besseli_continued_fractions_scaled(nu, x))
     else
         return besseli_large_orders_scaled(nu, x)
     end
@@ -194,7 +195,6 @@ 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)
@@ -203,6 +203,14 @@ function _besseli_continued_fractions(nu, x::T) where T
     (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)
+    knu, knum1 = up_recurrence(xx, besselk0x(xx), besselk1x(xx), nu)
+    # 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)
@@ -220,7 +228,6 @@ function steed(n, x::T) where T
     end
     return h
 end
-
 function besseli_large_argument(v, z::T) where T
     MaxIter = 1000
     a = exp(z / 2)
@@ -239,6 +246,23 @@ function besseli_large_argument(v, z::T) where T
     end
     return res * coef * a
 end
+function besseli_large_argument_scaled(v, z::T) where T
+    MaxIter = 1000
+    coef = inv(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
+end
 
 function besseli_small_arguments(v, z::T) where T
     S = promote_type(T, Float64)

src/recurrence.jl

@@ -1,3 +1,5 @@
+# no longer used for besseli but could be used in future for Jn, Yn
+#=
 @inline function down_recurrence(x, in, inp1, nu, branch)
     # this prevents us from looping through large values of nu when the loop will always return zero
     (iszero(in) || iszero(inp1)) && return zero(x)
@@ -13,6 +15,7 @@
     end
     return inm1
 end
+=#
 @inline function up_recurrence(x, k0, k1, nu)
     nu == 0 && return k0
     nu == 1 && return k1

test/besseli_test.jl

@@ -87,6 +87,7 @@ t = [besseli(m, x) for m in m, x in x]
 @test t[10] isa Float64
 @test t ≈ [SpecialFunctions.besseli(m, x) for m in m, x in x]
 
-@test besselix(10, 2.0) ≈ SpecialFunctions.besselix(10, 2.0)
-@test besselix(100, 14.0) ≈ SpecialFunctions.besselix(100, 14.0)
-@test besselix(120, 504.0) ≈ SpecialFunctions.besselix(120, 504.0)
+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]
+