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,119 @@ 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
+
+# originally contribued by @cgeoga https://github.com/cgeoga/BesselK.jl/blob/main/src/besk_ser.jl
+function _besselk_ser(v, x, maxit, tol, modify)
+  T    = promote_type(typeof(x), typeof(v))
+  out  = zero(T)
+  oneT = one(T)
+  twoT = oneT+oneT
+  # precompute a handful of things:
+  xd2      = x/twoT
+  xd22     = xd2*xd2
+  half     = oneT/twoT
+  if modify
+    e2v      = exp2(v)
+    xd2_v    = (x^(2*v))/e2v
+    xd2_nv   = e2v
+  else
+    lxd2     = log(xd2)
+    xd2_v    = exp(v*lxd2)
+    xd2_nv   = exp(-v*lxd2)
+  end
+  gam_v    = gamma(v)
+  gam_nv   = gamma(-v)
+  gam_1mv  = -gam_nv*v # == gamma(one(T)-v)
+  gam_1mnv = gam_v*v   # == gamma(one(T)+v)
+  xd2_pow  = oneT
+  fact_k   = oneT
+  floatk   = convert(T, 0)
+  (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
+  # now the loop using Oana's series expansion, with term function manually
+  # inlined for max speed:
+  for _j in 0:maxit
+    t1   = half*xd2_pow
+    tmp  = _t1/(gmv*fact_k)
+    tmp += _t2/(gpv*fact_k)
+    term = t1*tmp
+    out += term
+    ((abs(term) < tol) && _j>5) && 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*(oneT+v+floatk), gmv*(oneT-v+floatk)) 
+    xd2_pow *= xd22
+    fact_k  *= (floatk+1)
+    floatk  += 1.0
+  end
+  throw(error("$maxit iterations reached without achieving atol $tol."))
+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
-=#
+