add wider test coverage for asym comparison

Author: heltonmc

src/besselk.jl

@@ -173,3 +173,33 @@ end
     end
     return k2
 end
+
+
+function k1(v, x::T) where T <: AbstractFloat
+    z = x / v
+    zs = sqrt(1 + z^2)
+ 
+    n = zs + log(z / (1 + zs))
+    coef = exp(-v * n) * sqrt(pi / (2*v)) / sqrt(zs)
+ 
+    p = inv(zs)
+ 
+    u0 = one(x)
+    u1 = p / 24 * (3 - 5*p^2) * -1 / v
+    u2 = p^2 / 1152 * (81 - 462*p^2 + 385*p^4) / v^2
+    u3 = p^3 / 414720 * (30375 - 369603 * p^2 + 765765*p^4 - 425425*p^6) * -1 / v^3
+    u4 = p^4 / 39813120 * (4465125 - 94121676*p^2 + 349922430*p^4 - 446185740*p^6 + 185910725*p^8) / v^4
+    u5 = p^5 / 6688604160 * (-188699385875*p^10 + 566098157625*p^8 - 614135872350*p^6 + 284499769554*p^4 - 49286948607*p^2 + 1519035525) * -1 / v^5
+    u6 = p^6 / 4815794995200 * (1023694168371875*p^12 - 3685299006138750*p^10 + 5104696716244125*p^8 - 3369032068261860*p^6 + 1050760774457901*p^4 - 127577298354750*p^2 + 2757049477875) * 1 / v^6
+    u7 = p^7 / 115579079884800 * (-221849150488590625*p^14 + 931766432052080625*p^12 - 1570320948552481125*p^10 + 1347119637570231525*p^8 - 613221795981706275*p^6 + 138799253740521843*p^4 - 12493049053044375*p^2 + 199689155040375) * -1 / v^7
+    u8 = p^8 / 22191183337881600 * (448357133137441653125*p^16 - 2152114239059719935000*p^14 + 4272845805510421639500*p^12 - 4513690624987320777000*p^10 + 2711772922412520971550*p^8 - 914113758588905038248*p^6 + 157768535329832893644*p^4 - 10960565081605263000*p^2 + 134790179652253125) * 1 / v^8
+    u9 = p^9 / 263631258054033408000 * (-64041091111686478524109375*p^18 + 345821892003106984030190625*p^16 - 790370708270219620781737500*p^14 + 992115946599792610768672500*p^12 - 741743213039573443221773250*p^10 + 334380732677827878090447630*p^8 - 87432034049652400520788332*p^6 + 11921080954211358275362500*p^4 - 659033454841709672064375*p^2 + 6427469716717690265625) * -1 / v^9
+    u10 = p^10 / 88580102706155225088000 * (290938676920391671935028890625*p^20 - 1745632061522350031610173343750*p^18 + 4513386761946134740461797128125*p^16 - 6564241639632418015173104205000*p^14 + 5876803711285273203043452095250*p^12 - 3327704366990695147540934069220*p^10 + 1177120360439828012193658602930*p^8 - 246750339886026017414509498824*p^6 + 27299183373230345667273718125*p^4 - 1230031256571145165088463750*p^2 + 9745329584487361980740625) * 1 / v^10
+    out1 = ((u10 + u9) + u8) + u7
+    out2 = ((out1 + u6 + u5) + u4) + u3
+    out = ((out2 + u2) + u1) + u0
+ 
+    return coef*out
+ end
+ 
+ 

test/besselk_test.jl

@@ -67,6 +67,18 @@ k1x_32 = besselk1x.(Float32.(x))
 @test besselk(100, 3.9) ≈ SpecialFunctions.besselk(100, 3.9)
 @test besselk(100, 234.0) ≈ SpecialFunctions.besselk(100, 234.0)
 
+# test small arguments and order
+m = 0:40; x = [1e-6; 1e-4; 1e-3; 1e-2; 0.1; 1.0:2.0:700.0]
+@test [besselk(m, x) for m in m, x in x] ≈ [SpecialFunctions.besselk(m, x) for m in m, x in x]
+
+# test medium arguments and order
+m = 30:200; x = 5.0:5.0:100.0
+@test [besselk(m, x) for m in m, x in x] ≈ [SpecialFunctions.besselk(m, x) for m in m, x in x]
+
+# test large orders
+m = 200:5:1000; x = 400.0:10.0:1200.0
+@test [besselk(m, x) for m in m, x in x] ≈ [SpecialFunctions.besselk(m, x) for m in m, x in x]
+
 @test iszero(besselk(20, 1000.0))
 @test isinf(besselk(250, 5.0))