Decided to tinker a bit with the rescaled definition near the regime

change from direct series to asymptotic strategies for the xv-rescaled
methods and managed to squeeze out a couple more digits. Now the
accuracies are pretty comparable and I think reasonably satisfying.

Also gave the testing file some TLC. Now just running the tests is
actually a very convenient way to check for worst-case accuracy
regressions, especially if you comment out the allocations test, which
takes a little while.

Bumping version now due to the improved accuracy. This is probably the
end of this little tinkering spell.

Project.toml

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

src/besk.jl

@@ -104,13 +104,15 @@ end
 function adbesselkxv(v, x, maxit=100, tol=1e-12, order=5)
   (iszero(v) && iszero(x)) && return Inf
   is_ad = !(v isa AbstractFloat)
-  xcut  = is_ad ? 8.5 : 2.0
+  xcut  = is_ad ? 6.0 : 2.0
   if !isinteger(v) && (abs(x) <= 1e-8) # use Taylor at zero.
     return besselkxv_t0(v, x)
   elseif (x < xcut) && (v < 5.75) && !isnearint(v, 0.01)
     return _besselk_ser(v, x, maxit, tol, true)
   elseif (x < xcut)
     return _besselk_temme(v, x, maxit, tol, true)
+  elseif is_ad && (x > xcut) && (x < 15.0)
+    return _besselk_asv(v, x, 12, true)
   else
     return adbesselk(v, x, maxit, tol, order)*(x^v)
   end

test/runtests.jl

@@ -25,40 +25,68 @@ ad2_dbesselkxv_dv_dv(v, x) = ForwardDiff.derivative(_v->ad_dbesselkxv_dv(_v, x),
 
 # direct accuracy:
 @testset "direct eval" begin
+  println("\nDirect evaluations:")
   for (ref_fn, cand_fn, case) in ((besselk,   BesselK._besselk, :standard),
                                   (besselkxv, BesselK.adbesselkxv, :rescaled))
     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
+    ix        = findall(x-> x <= 1000.0, amos_ref)
+    thresh    = case == :standard ? 5e-11 : 2e-12
+    (maxerr, maxix) = findmax(abs, atols[ix])
+    (maxerr_v, maxerr_x) = VX[ix][maxix]
+    println("Case $case:")
+    println("worst (v,x):  ($maxerr_v, $maxerr_x)")
+    println("Ref value:     $(amos_ref[ix][maxix])")
+    println("Est value:     $(candidate[ix][maxix])")
+    println("Abs error:     $(round(maxerr, sigdigits=3))")
+    @test maxerr < thresh
   end
+  println()
 end
 
 # test derivative accuracy:
 @testset "first derivative" begin
+  println("\nFirst derivatives:")
   for (ref_fn, cand_fn, case) in ((fd_dbesselk_dv, ad_dbesselk_dv, :standard),
-                                  (fd_dbesselkxv_dv, ad_dbesselkxv_dv, :rescale))
+                                  (fd_dbesselkxv_dv, ad_dbesselkxv_dv, :rescaled))
     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]
-    thresh    = case == :standard ? 1e-8 : 1e-5
-    @test maximum(atols_f) < thresh
+    ix        = findall(x-> x <= 1000.0, amos_ref)
+    thresh    = case == :standard ? 4e-9 : 8e-7
+    (maxerr, maxix) = findmax(abs, atols[ix])
+    (maxerr_v, maxerr_x) = VX[ix][maxix]
+    println("Case $case:")
+    println("worst (v,x):  ($maxerr_v, $maxerr_x)")
+    println("Ref value:     $(amos_ref[ix][maxix])")
+    println("Est value:     $(candidate[ix][maxix])")
+    println("Abs error:     $(round(maxerr, sigdigits=3))")
+    @test maxerr < thresh
   end
+  println()
 end
 
 # test second derivative accuracy:
 @testset "second derivative" begin
+  println("\nSecond derivatives:")
   for (ref_fn, cand_fn, case) in ((fd2_dbesselk_dv_dv, ad2_dbesselk_dv_dv, :standard),
-                                  (fd2_dbesselkxv_dv_dv, ad2_dbesselkxv_dv_dv, :rescale))
+                                  (fd2_dbesselkxv_dv_dv, ad2_dbesselkxv_dv_dv, :rescaled))
     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 .<= 100.0]
-    thresh    = case == :standard ? 1e-6 : 1e-4
-    @test maximum(atols_f) < thresh
+    ix        = findall(x-> x <= 100.0, amos_ref)
+    thresh    = case == :standard ? 5e-7 : 2e-6
+    (maxerr, maxix) = findmax(abs, atols[ix])
+    (maxerr_v, maxerr_x) = VX[ix][maxix]
+    println("Case $case:")
+    println("worst (v,x):  ($maxerr_v, $maxerr_x)")
+    println("Ref value:     $(amos_ref[ix][maxix])")
+    println("Est value:     $(candidate[ix][maxix])")
+    println("Abs error:     $(round(maxerr, sigdigits=3))")
+    @test maxerr < thresh
   end
+  println()
 end
 
 # Testing the _xv versions really slows down the test script, and in general