update cutoffs and add tests

Author: heltonmc

src/U_polynomials.jl

@@ -26,7 +26,7 @@ function besseljy_debye(v, x)
 
     return coef_Jn * Uk_Jn, coef_Yn * Uk_Yn
 end
-hankel_debye_cutoff(nu, x) = nu < -0.0657 + 0.9976*x + Base.Math._approx_cbrt(-276.915*x)
+hankel_debye_cutoff(nu, x) = nu < 0.2 + x + Base.Math._approx_cbrt(-411*x)
 
 # valid when v < x (uniform asymptotic expansions)
 """

src/bessely.jl

@@ -139,8 +139,27 @@ function _bessely1_compute(x::Float64)
         p = p * sc[1] + w * q * sc[2]
         return p * SQ2OPI(T) / sqrt(x)
     else
-      
-        return besseljy_large_argument(1.0, x)[2]
+        xinv = inv(x)
+        x2 = xinv*xinv
+        if x < 130.0
+            p1 = (one(T), 3/16, -99/512, 6597/8192, -4057965/524288, 1113686901/8388608, -951148335159/268435456, 581513783771781/4294967296)
+            q1 = (3/8, -21/128, 1899/5120, -543483/229376, 8027901/262144, -30413055339/46137344, 9228545313147/436207616)
+            p = evalpoly(x2, p1)
+            q = evalpoly(x2, q1)
+        else
+            p2 = (one(T), 3/16, -99/512, 6597/8192)
+            q2 = (3/8, -21/128, 1899/5120, -543483/229376)
+            p = evalpoly(x2, p2)
+            q = evalpoly(x2, q2)
+        end
+
+        a = SQ2OPI(T) * sqrt(xinv) * p
+        xn = muladd(xinv, q, - 3 * PIO4(T))
+
+        # the following computes b = sin(x + xn) more accurately
+        # see src/misc.jl
+        b = sin_sum(x, xn)
+        return a * b
     end
 end
 function _bessely1_compute(x::Float32)
@@ -172,18 +191,18 @@ function _bessely(nu, x::T) where T
     nu == 0 && return bessely0(x)
     nu == 1 && return bessely1(x)
 
-    # large argument branch see src/asymptotics.jl
-    besseljy_large_argument_cutoff(nu, x) && return besseljy_large_argument(nu, x)[2]
-
     # x < ~nu branch see src/U_polynomials.jl
     besseljy_debye_cutoff(nu, x) && return besseljy_debye(nu, x)[2]
 
     # x > ~nu branch see src/U_polynomials.jl on computing Hankel function
     hankel_debye_cutoff(nu, x) && return imag(hankel_debye(nu, x))
 
-    # use forward recurrence if nu is an integer up until it becomes inefficient 
+    # large argument branch see src/asymptotics.jl
+    besseljy_large_argument_cutoff(nu, x) && return besseljy_large_argument(nu, x)[2]
+
+    # use forward recurrence if nu is an integer up until it becomes inefficient
     (isinteger(nu) && nu < 250) && return besselj_up_recurrence(x, bessely1(x), bessely0(x), 1, nu)[1]
-  
+
     # use power series for small x and for when nu > x
     bessely_series_cutoff(nu, x) && return bessely_power_series(nu, x)
 
@@ -192,7 +211,7 @@ function _bessely(nu, x::T) where T
 
     # at this point x > 19.0 (for Float64) and fairly close to nu
     # shift nu down and use the debye expansion for Hankel function (valid x > nu) then use forward recurrence
-    nu_shift = floor(nu) - ceil(Int, -4.0657 + 0.9976*x + Base.Math._approx_cbrt(-276.915*x))
+    nu_shift = floor(nu) - ceil(Int, -1.5 + x + Base.Math._approx_cbrt(-411*x))
     v2 = nu - maximum((nu_shift, modf(nu)[1] + 1))
     return besselj_up_recurrence(x, imag(hankel_debye(v2, x)), imag(hankel_debye(v2 - 1, x)), v2, nu)[1]
 end

test/bessely_test.jl

@@ -1,5 +1,5 @@
 # general array for testing input to SpecialFunctions.jl
-x = 0.01:0.01:100.0
+x = 0.01:0.01:150.0
 
 ### Tests for bessely0
 y0_SpecialFunctions = SpecialFunctions.bessely0.(big.(x))  # array to be tested against computed in BigFloats
@@ -65,6 +65,23 @@ y1_32 = bessely1.(Float32.(x))
 @test bessely1(Inf32) == zero(Float32)
 @test bessely1(Inf64) == zero(Float64)
 
-
-# briefly test the large argument is working
-@test Bessels.besseljy_large_argument(10.0, 100.0)[2] ≈ SpecialFunctions.bessely(10.0, 100.0)
+## Tests for bessely
+
+## test all numbers and orders for 0<nu<100
+x = [0.05, 0.1, 0.2, 0.4, 0.5, 0.6, 0.7, 0.75, 0.8, 0.85, 0.9, 0.92, 0.95, 0.97, 0.99, 1.0, 1.01, 1.05, 1.1, 1.2, 1.4, 2.0]
+nu = [2, 4, 6, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100]
+for v in nu, xx in x
+    xx *= v
+    sf = SpecialFunctions.bessely(BigFloat(v), BigFloat(xx))
+    @test isapprox(Bessels._bessely(v, xx), Float64(sf), rtol=7e-14)
+end
+
+# test decimal orders
+# SpecialFunctions.jl can give errors over 1e-12 so need to soften tolerance to match
+# need to switch tests over to ArbNumerics.jl for better precision tests 
+x = [0.05, 0.1, 0.2, 0.25, 0.3, 0.4, 0.5,0.55,  0.6,0.65,  0.7, 0.75, 0.8, 0.85, 0.9, 0.92, 0.95, 0.97, 0.99, 1.0, 1.01, 1.05, 1.08, 1.1, 1.2, 1.4, 1.5, 1.6, 1.8, 2.0, 2.5, 3.0, 4.0, 4.5, 4.99, 5.1]
+nu = [0.1, 0.4567, 0.8123, 1.5, 2.5, 4.1234, 6.8, 12.3, 18.9, 28.2345, 38.1235, 51.23, 72.23435, 80.5, 98.5, 104.2]
+for v in nu, xx in x
+    xx *= v
+    @test isapprox(Bessels._bessely(v, xx), SpecialFunctions.bessely(v, xx), rtol=5e-12)
+end