add continued fractions

Author: heltonmc

src/besseli.jl

@@ -193,3 +193,25 @@ 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)
+    knu, knum1 = up_recurrence(x, besselk0(x), besselk1(x), nu)
+    return 1 / (x * (knum1 + knu / steed(nu, x)))
+end
+
+function steed(n, x::T) where T
+    xinv = inv(x)
+    xinv2 = 2 * xinv
+    d = x / (n + n)
+    a = d
+    h = a
+    b = muladd(2, n, 2) * xinv
+    for _ in 1:100000
+        d = inv(b + d)
+        a = muladd(b, d, -1) * a
+        h = h + a
+        b = b + xinv2
+        abs(a / h) <= eps(T) && break
+    end
+    return h
+end

src/besselk.jl

@@ -161,7 +161,7 @@ Modified Bessel function of the second kind of order nu, ``K_{nu}(x)``.
 function besselk(nu, x::T) where T <: Union{Float32, Float64, BigFloat}
     T == Float32 ? branch = 20 : branch = 50
     if nu < branch
-        return up_recurrence(x, besselk0(x), besselk1(x), nu)
+        return up_recurrence(x, besselk0(x), besselk1(x), nu)[1]
     else
         return besselk_large_orders(nu, x)
     end
@@ -175,7 +175,7 @@ Scaled modified Bessel function of the second kind of order nu, ``K_{nu}(x)*e^{x
 function besselkx(nu::Int, x::T) where T <: Union{Float32, Float64}
     T == Float32 ? branch = 20 : branch = 50
     if nu < branch
-        return up_recurrence(x, besselk0x(x), besselk1x(x), nu)
+        return up_recurrence(x, besselk0x(x), besselk1x(x), nu)[1]
     else
         return besselk_large_orders_scaled(nu, x)
     end

src/recurrence.jl

@@ -28,5 +28,5 @@ end
         k0 = k1
         k1 = k2
     end
-    return k2
+    return k2, k0
 end