MH's besselk asymptotic expansion.

Author: cgeoga

src/besselk.jl

@@ -233,7 +233,7 @@ function besselk_positive_args(nu, x::T) where T <: Union{Float32, Float64}
     (isinteger(nu-1/2) && !debye_cut) && return sphericalbesselk(nu-1/2, x)*SQRT_PID2(T)*sqrt(x)
 
     # check if the standard asymptotic expansion can be used:
-    besselk_asexp_cutoff(nu, x) && return besselk_asymptoticexp(nu, x, 10)
+    besselk_asexp_cutoff(nu, x) && return besselk_large_argument(nu, x)
 
     # use uniform debye expansion if x or nu is large
     debye_cut && return besselk_large_orders(nu, x)
@@ -442,55 +442,28 @@ Computes the asymptotic expansion of K_ν w.r.t. argument. Accurate for large x,
 """
 besselk_asexp_cutoff(nu, x::T) where T = (nu < 20.0) && (x > ASEXP_CUTOFF(T))
 
-# Compute the function for (_v, _v+1), where _v = v-floor(v), and then do
-# forward recursion. If the argument is really large, this should be faster than
-# the uniform asymptotic expansion and plenty accurate.
-function besselk_asymptoticexp(v, x::T, order) where T
-    _v  = v - floor(v)
-    (kv, kvp1) = _besselk_as_pair(_v, x, order)
-    v == _v && return kv
-    _v += one(v)
-    twodx = T(2)*inv(x)
-    while _v < v
-        (kv, kvp1) = (kvp1, muladd(kvp1, twodx*_v, kv))
-        _v += one(_v)
-    end
-    kvp1
+function besselk_large_argument(v, x::T) where T
+    a = exp(-x / 2)
+    coef = a * sqrt(pi / 2x)
+    return T(_besselk_large_argument(v, x) * coef * a)
 end
 
-# For now, no exponential improvement. It requires the exponential integral
-# function, which would either need to be lifted from SpecialFunctions.jl or
-# re-implemented. And with an order of, like, 10, this seems to be pretty
-# accurate and still faster than the uniform asymptotic expansion.
-function _besselk_as_pair(v, x::T, order) where T
-    fv = 4*v*v
-    fvp1 = 4*(v+one(v))^2
-    _z = x
-    ser_v = one(T)
-    ser_vp1 = one(T)
-    floatj = one(T)
-    ak_numv = fv   - floatj
-    ak_numvp1 = fvp1 - floatj
-    factj = one(T)
-    twofloatj = one(T)
-    eightj = T(8)
-    for _ in 1:order
-        # add to the series:
-        term_v = ak_numv/(factj*_z*eightj)
-        ser_v += term_v
-        term_vp1 = ak_numvp1/(factj*_z*eightj)
-        ser_vp1 += term_vp1
-        # update ak and _z:
-        floatj += one(T)
-        twofloatj += T(2)
-        factj  *= floatj
-        fourfloatj = twofloatj*twofloatj
-        ak_numv *= (fv   - fourfloatj)
-        ak_numvp1 *= (fvp1 - fourfloatj)
-        _z *= x
-        eightj *= T(8)
-    end
-    pre_multiply = SQRT_PID2(T)*exp(-x)/sqrt(x)
-    (pre_multiply*ser_v, pre_multiply*ser_vp1)
+besselk_large_argument_scaled(v, x::T) where T =  T(_besselk_large_argument(v, x) * sqrt(pi / 2x))
+
+function _besselk_large_argument(v, x::T) where T
+    MaxIter = 5000 
+    S = promote_type(T, Float64) 
+    v, x = S(v), S(x) 
+ 
+    fv2 = 4 * v^2 
+    term = one(S) 
+    res = term 
+    s = term 
+    for i in 1:MaxIter 
+        offset = muladd(2, i, -1) 
+        term *= muladd(offset, -offset, fv2) / (8 * x * i) 
+        res = muladd(term, s, res) 
+        abs(term) <= eps(T) && break 
+    end 
+    return res 
 end
-