Besselj and Besselk for all integer and noninteger postive orders

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,114 @@ 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)
+    elseif x < 2.0 || nu > 1.6x - 1.0
+        return besselk_power_series(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
+
+# originally contributed by @cgeoga https://github.com/cgeoga/BesselK.jl/blob/main/src/besk_ser.jl
+# Equation 3.2 from Geoga, Christopher J., et al. "Fitting Mat\'ern Smoothness Parameters Using Automatic Differentiation." 
+# arXiv preprint arXiv:2201.00090 (2022).
+function besselk_power_series(v, x::T) where T
+    MaxIter = 1000
+    # precompute a handful of things:
+    xd2  = x / 2
+    xd22 = xd2 * xd2
+    half = one(T) / 2
+    # (x/2)^(±v). Writing the literal power doesn't seem to change anything here,
+    # and I think this is faster:
+    lxd2 = log(xd2)
+    xd2_v = exp(v*lxd2)
+    xd2_nv = exp(-v*lxd2)
+    # use the gamma function a couple times to start:
+    gam_v = gamma(v)
+    xp1 = abs(v) + one(T)
+    gam_nv = π / sinpi(xp1) / _gamma(xp1)
+    gam_1mv = -gam_nv*v # == gamma(one(T)-v)
+    gam_1mnv = gam_v*v   # == gamma(one(T)+v)
+    (gpv, gmv) = (gam_1mnv, gam_1mv)
+    # One final re-compression of a few things:
+    _t1 = gam_v*xd2_nv*gam_1mv
+    _t2 = gam_nv*xd2_v*gam_1mnv
+    # A couple series-specific accumulators:
+    (xd2_pow, fact_k, floatk, out) = (one(T), one(T), zero(T), zero(T))
+    for _ in 0:MaxIter
+      t1 = half*xd2_pow
+      tmp = _t1/(gmv*fact_k)
+      tmp += _t2/(gpv*fact_k)
+      term = t1*tmp
+      out += term
+      abs(term/out) < eps(T) && return out
+      # Use the trick that gamma(1+k+1+v) == gamma(1+k+v)*(1+k+v) to skip gamma calls:
+      (gpv, gmv) = (gpv*(one(T)+v+floatk), gmv*(one(T)-v+floatk)) 
+      xd2_pow *= xd22
+      fact_k *= (floatk+1)
+      floatk += one(T)
+    end
+    return out
+  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
-=#
+