continued fraction

Author: heltonmc

src/U_polynomials.jl

@@ -70,10 +70,10 @@ function Uk_poly_Hankel(p, v, p2, x::T) where T <: Float64
 end
 Uk_poly_Hankel(p, v, p2, x) = Uk_poly_Jn(p, v, p2, x::BigFloat)
 
-Uk_poly_In(p, v, p2, ::Type{Float32}) = Uk_poly3(p, v, p2)[1]
-Uk_poly_In(p, v, p2, ::Type{Float64}) = Uk_poly5(p, v, p2)[1]
-Uk_poly_Kn(p, v, p2, ::Type{Float32}) = Uk_poly3(p, v, p2)[2]
-Uk_poly_Kn(p, v, p2, ::Type{Float64}) = Uk_poly5(p, v, p2)[2]
+Uk_poly_In(p, v, p2, ::Type{Float32}) = Uk_poly5(p, v, p2)[1]
+Uk_poly_In(p, v, p2, ::Type{Float64}) = Uk_poly10(p, v, p2)[1]
+Uk_poly_Kn(p, v, p2, ::Type{Float32}) = Uk_poly5(p, v, p2)[2]
+Uk_poly_Kn(p, v, p2, ::Type{Float64}) = Uk_poly10(p, v, p2)[2]
 
 @inline function split_evalpoly(x, P)
     # polynomial P must have an even number of terms

src/besseli.jl

@@ -139,18 +139,15 @@ Nu must be real.
 function besseli(nu, x::T) where T <: Union{Float32, Float64}
     nu == 0 && return besseli0(x)
     nu == 1 && return besseli1(x)
-    
-    if x > maximum((T(30), nu^2 / 4))
+
+    if x > maximum((T(30), nu^2 / 6))
         return T(besseli_large_argument(nu, x))
-    elseif x <= 2 * sqrt(nu + 1)
-        return T(besseli_small_arguments(nu, x))
-    elseif nu < 100
-        return T(_besseli_continued_fractions(nu, x))
-    else
+    elseif nu > 25.0 || x > 35.0
         return T(besseli_large_orders(nu, x))
+    else
+        return T(besseli_power_series(nu, x))
     end
 end
-
 """
     besselix(nu, x::T) where T <: Union{Float32, Float64}
 
@@ -264,25 +261,16 @@ function besseli_large_argument_scaled(v, z::T) where T
     return res * coef
 end
 
-function besseli_small_arguments(v, z::T) where T
-    S = promote_type(T, Float64)
-    x = S(z)
-    if v < 20
-        coef = (x / 2)^v / factorial(v)
-    else
-        vinv = inv(v)
-        coef = sqrt(vinv / (2 * π)) * MathConstants.e^(v * (log(x / (2 * v)) + 1)) 
-        coef *= evalpoly(vinv, (1, -1/12, 1/288,  139/51840, -571/2488320, -163879/209018880, 5246819/75246796800, 534703531/902961561600))
-    end
-
-    MaxIter = 1000
-    out = one(S)
-    zz = x^2 / 4
-    a = one(S)
-    for k in 1:MaxIter
-        a *= zz / (k * (k + v))
+function besseli_power_series(v, x::T) where T
+    MaxIter = 3000
+    out = zero(T)
+    xs = (x/2)^v
+    a = xs / gamma(v + one(T))
+    t2 = (x/2)^2
+    for i in 0:MaxIter
         out += a
-        a <= eps(T) && break
+        abs(a) < eps(T) * abs(out) && break
+        a *= inv((v + i + one(T)) * (i + one(T))) * t2
     end
-    return coef * out
+    return out
 end

src/besselk.jl

@@ -153,6 +153,7 @@ In other words, for large values of x and small to moderate values of nu,
 this routine will underflow to zero.
 For small to medium values of x and large values of nu this will overflow and return Inf.
 =#
+#=
 """
     besselk(x::T) where T <: Union{Float32, Float64}
 
@@ -166,7 +167,7 @@ function besselk(nu, x::T) where T <: Union{Float32, Float64, BigFloat}
         return besselk_large_orders(nu, x)
     end
 end
-
+=#
 """
     besselk(x::T) where T <: Union{Float32, Float64}
 
@@ -203,4 +204,70 @@ function besselk_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64,
     p2  = v^2/fma(max(v,x), max(v,x), min(v,x)^2)
 
     return T(coef*Uk_poly_Kn(p, v, p2, T))
-end
+end
+
+function besselk(nu, x::T) where T <: Union{Float32, Float64, BigFloat}
+    (isinteger(nu) && nu < 250) && return besselk_up_recurrence(x, besselk1(x), besselk0(x), 1, nu)[1]
+
+    if nu > 25.0 || x > 35.0
+        return besselk_large_orders(nu, x)
+    else
+        return besselk_continued_fraction(nu, x)
+    end
+end
+
+# could also use the continued fraction for inu/inmu1
+# but it seems like adapting the besseli_power series
+# to give both nu and nu-1 is faster
+
+#inum1 = besseli_power_series(nu-1, x)
+#H_inu = steed(nu, x)
+#inu = besseli_power_series(nu, x)#inum1 * H_inu
+function besselk_continued_fraction(nu, x)
+    inu, inum1 = besseli_power_series_inu_inum1(nu, x)
+    H_knu = besselk_ratio_knu_knup1(nu-1, x)
+    return 1 / (x * (inum1 + inu / H_knu))
+end
+
+function besseli_power_series_inu_inum1(v, x::T) where T
+    MaxIter = 3000
+    out = zero(T)
+    out2 = zero(T)
+    x2 = x / 2
+    xs = (x2)^v
+    gmx = xs / gamma(v)
+    a = gmx / v
+    b = gmx / x2
+    t2 = x2*x2
+    for i in 0:MaxIter
+        out += a
+        out2 += b
+        abs(a) < eps(T) * abs(out) && break
+        a *= inv((v + i + one(T)) * (i + one(T))) * t2
+        b *= inv((v + i) * (i + one(T))) * t2
+    end
+    return out, out2
+end
+
+
+# slightly modified version of https://github.com/heltonmc/Bessels.jl/issues/17#issuecomment-1195726642 from @cgeoga
+function besselk_ratio_knu_knup1(v, x::T) where T
+    MaxIter = 1000
+    # do the modified Lentz method:
+    (hn, Dn, Cn) = (1e-50, zero(v), 1e-50)
+
+    jf   = one(T)
+    vv   = v*v
+    for _ in 1:MaxIter
+        an  = (vv - ((2*jf - 1)^2) * T(0.25))
+        bn  = 2 * (x + jf)
+        Cn  = an / Cn + bn      
+        Dn  = inv(muladd(an, Dn, bn))
+        del = Dn * Cn
+        hn *= del
+        abs(del - 1) < eps(T) && break
+        jf += one(T)
+    end
+    xinv = inv(x)
+    return xinv * (v + x + 1/2) + xinv * hn
+end

src/recurrence.jl

@@ -36,7 +36,7 @@ end
     return jnup1, jnu
 end
 
-#=
+
 # currently not used
 # backward recurrence relation for besselk and besseli
 # outputs both (bessel(x, nu_end), bessel(x, nu_end-1)
@@ -50,4 +50,4 @@ end
     end
     return jnup1, jnu
 end
-=#
+