Catch a few very subtle edge cases in the xv modifications, and added

them to the tests. Bumped all the deps.

Project.toml

@@ -1,7 +1,7 @@
 name = "BesselK"
 uuid = "432ab697-7a72-484f-bc4a-bc531f5c819b"
 authors = ["Chris Geoga <christopher.geoga@rutgers.edu>"]
-version = "0.1.1"
+version = "0.2.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/besk.jl

@@ -32,7 +32,8 @@ function _besselk(v, x, maxit=100, tol=1e-12, order=6)
     if is_ad
       return _besselk_ser(v, x, maxit, tol, false) 
     else
-      return _besselk_as(v, x, Int(ceil(v)), is_ad) 
+      # Probably don't need the is_ad correction here.
+      return _besselk_as(v, x, Int(ceil(v)), false) 
     end
   end
   #
@@ -58,7 +59,7 @@ function _besselk(v, x, maxit=100, tol=1e-12, order=6)
     if abs(v) > 1.5
       return _besselk_asv(v,x,8) 
     else
-      return _besselk_as(v,x,order,is_ad) 
+      return _besselk_as(v,x,order) 
     end
   end
 end
@@ -94,13 +95,18 @@ end
   evalpoly(x, cof)
 end
 
+# Note that the special series cutofs are low compared to the above function. In
+# general, the adbesselk* functions that are exported really should be pretty
+# near machine precision here or should just fall back to AMOS when possible. It
+# turns out that the *xv modifications in the code are really only helpful when
+# the argument is pretty small.
 function adbesselkxv(v, x, maxit=100, tol=1e-12, order=5)
   (iszero(v) && iszero(x)) && return Inf
   if !isinteger(v) && (abs(x) <= 1e-8) # use Taylor at zero.
     return besselkxv_t0(v, x)
-  elseif (abs(x) < 8.5) && (v < 3.0) && !isnearint(v, 0.01)
+  elseif (abs(x) < 2.0) && (v < 3.0) && !isnearint(v, 0.01)
     return _besselk_ser(v, x, maxit, tol, true)
-  elseif (abs(x) < 8.5)
+  elseif (abs(x) < 2.0)
     return _besselk_temme(v, x, maxit, tol, true)
   else
     return adbesselk(v, x, maxit, tol, order)*(x^v)

src/besk_ser.jl

@@ -30,13 +30,13 @@ function _besselk_ser(v, x, maxit, tol, modify)
   _t2 = gam_nv*xd2_v*gam_1mnv
   # now the loop using Oana's series expansion, with term function manually
   # inlined for max speed:
-  for _ in 0:maxit
+  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 && return out
+    ((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

test/runtests.jl

@@ -9,6 +9,8 @@ const REF_FD2 = central_fdm(10,2)
 
 atolfun(tru, est) = isnan(est) ? NaN : (isinf(tru) ? 0.0 : abs(tru-est))
 
+besselkxv(v,x) = besselk(v,x)*(x^v)
+
 fd_dbesselk_dv(v, x) = REF_FD1(_v->besselk(_v, x), v)
 ad_dbesselk_dv(v, x) = ForwardDiff.derivative(_v->BesselK.adbesselk(_v, x), v)
 
@@ -17,11 +19,14 @@ ad2_dbesselk_dv_dv(v, x) = ForwardDiff.derivative(_v->ad_dbesselk_dv(_v, x), v)
 
 # direct accuracy:
 @testset "direct eval" begin
-  amos_ref  = map(vx->besselk(vx[1], vx[2]), VX)
-  candidate = map(vx->BesselK._besselk(vx[1], vx[2]), VX)
-  atols     = map(a_c->atolfun(a_c[1], a_c[2]), zip(amos_ref, candidate))
-  atols_f   = atols[amos_ref .<= 1000.0]
-  @test maximum(atols_f) < 1e-10
+  for (ref_fn, cand_fn) in ((besselk,   BesselK._besselk),
+                            (besselkxv, BesselK.adbesselkxv))
+    amos_ref  = map(vx->ref_fn(vx[1], vx[2]), VX)
+    candidate = map(vx->cand_fn(vx[1], vx[2]), VX)
+    atols     = map(a_c->atolfun(a_c[1], a_c[2]), zip(amos_ref, candidate))
+    atols_f   = atols[amos_ref .<= 1000.0]
+    @test maximum(atols_f) < 1e-10
+  end
 end
 
 # test derivative accuracy: