Merge pull request #48 from cgeoga/master

Asymptotic expansion for large arg and small-ish order

Author: heltonmc

src/besselk.jl

@@ -232,6 +232,9 @@ function besselk_positive_args(nu, x::T) where T <: Union{Float32, Float64}
     # check if nu is a half-integer:
     (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_large_argument(nu, x)
+
     # use uniform debye expansion if x or nu is large
     debye_cut && return besselk_large_orders(nu, x)
 
@@ -431,3 +434,36 @@ end
 besselk_power_series_cutoff(nu, x::Float64) = x < 2.0 || nu > 1.6x - 1.0
 besselk_power_series_cutoff(nu, x::Float32) = x < 10.0f0 || nu > 1.65f0*x - 8.0f0
 
+
+"""
+  besselk_asymptoticexp(v, x, order)
+
+Computes the asymptotic expansion of K_ν w.r.t. argument. Accurate for large x, and faster than uniform asymptotic expansion for small to small-ish orders. The default order of the expansion in `Bessels.besselk` is 10.
+"""
+besselk_asexp_cutoff(nu, x::T) where T = (nu < 20.0) && (x > ASEXP_CUTOFF(T))
+
+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
+
+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

src/constants.jl

@@ -287,3 +287,6 @@ const Q3_k1(::Type{Float64}) = (
 
 const SQRT_PID2(::Type{Float64}) = 1.2533141373155003
 const SQRT_PID2(::Type{Float32}) = 1.2533141f0
+
+const ASEXP_CUTOFF(::Type{Float64}) = 35.0
+const ASEXP_CUTOFF(::Type{Float32}) = 35.0f0 # temporary