Fix Issue #27

README.md

@@ -20,6 +20,12 @@ zeros if the order `nu` is one of the first few values of `0, 1, ...` and the en
 of the first values of `1, 2, 3, ...`.  See the individual function docstrings for the actual
 extents of the lookup tables.
 
+### Limitations
+The first ten zeros of any of the four functions treated here are always found
+correctly for any choice of `nu`.  For `nu ≤ 150`, any choice for `n` will produce a 
+correct result. However, for `nu > 150` and `n > 10` the results should not be trusted as some of the 
+zeros immediately following the tenth may be skipped.
+
 ### Exported Functions
 
 #### besselj_zero(nu, n)

src/FunctionZeros.jl

@@ -54,6 +54,9 @@ Asymptotic formula for the `n`th zero of the the Bessel J (Y) function of order
 """
 function bessel_zero_asymptotic(nu_in::Real, n::Integer, kind=1)
     nu = abs(nu_in)
+    if (nu ≥ 33 && n ≤ 10) || (nu ≥ 30 && n ≤ 9) || (nu ≥ 26 && n ≤ 8) || (nu ≥ 25 && n ≤ 7)
+        return bessel_zero_largenu_asymptotic(nu, n, kind, false)
+    end
     if kind == 1
         beta = MathConstants.pi * (n + nu / 2 - 1//4)
     else # kind == 2
@@ -76,6 +79,8 @@ function bessel_zero_asymptotic(nu_in::Real, n::Integer, kind=1)
     return zero_asymp
 end
 
+bessel_zero_asymptotic(nu::Integer, n::Integer, kind=1) = bessel_zero_asymptotic(float(nu), n, kind) # fix #27
+
 # Use the asymptotic values as starting values.
 # These find the correct zeros even for n = 1,2,...
 # Order 0 is 6 times slower and 50-100 times less accurate
@@ -189,8 +194,12 @@ Asymptotic formula for the `n`th zero of the the derivative of Bessel J (Y) func
 of order `nu`. `kind == 1 (2)` for Bessel function of the first (second) kind, J (Y).
 """
 function bessel_deriv_zero_asymptotic(nu_in::Real, n::Integer, kind=1)
-    # Reference: https://dlmf.nist.gov/10.21.E20
     nu = abs(nu_in)
+    if (nu ≥ 33 && n ≤ 10) || (nu ≥ 30 && n ≤ 9) || (nu ≥ 26 && n ≤ 8) || (nu ≥ 25 && n ≤ 7)
+        return bessel_zero_largenu_asymptotic(nu, n, kind, true)
+    end
+
+    # Reference: https://dlmf.nist.gov/10.21.E20
     if kind == 1
         beta = MathConstants.pi * (n + nu / 2 - 3//4)
     else # kind == 2
@@ -214,6 +223,8 @@ function bessel_deriv_zero_asymptotic(nu_in::Real, n::Integer, kind=1)
     return zero_asymp
 end
 
+bessel_deriv_zero_asymptotic(nu::Integer, n::Integer, kind=1) = bessel_deriv_zero_asymptotic(float(nu), n, kind) # fix #27
+
 """
     _besselj_deriv_zero(nu, n)
 
@@ -293,4 +304,96 @@ function bessely_deriv_zero(nu::Union{Integer,Float64}, n::Integer)
     end
 end
 
+"""
+    airy_zeros()
+
+Return the first few negative zeros of the functions `airyai`, `airybi`, `airyaiprime`, and `airybiprime`
+
+## Return Value
+- `(; ai, bi, aiprime, biprime)`: A named tuple containing the first few negative zeros of the functions
+  `airyai`, `airybi`, `airyaiprime`, and `airybiprime`, respectively, as defined in the `SpecialFunctions`
+   package. Each field in the named tuple consists of a tuple of `n` increasingly negative values. Here `n`
+   is a small integer, currently 10. 
+"""
+@inline function airy_zeros()
+    ai = (-2.338107410459767, -4.087949444130973, -5.520559828095556, -6.7867080900717625, 
+          -7.94413358712085, -9.02265085334098, -10.040174341558084, -11.008524303733264,
+          -11.936015563236262, -12.828776752865759)
+    bi = (-1.173713222709127, -3.2710933028363516, -4.8307378416620095, -6.169852128310234,
+          -7.3767620793677535, -8.491948846509374, -9.538194379346248, -10.529913506705357,
+          -11.47695355127878, -12.38641713858274)
+    aiprime = (-1.0187929716474717, -3.248197582179841, -4.820099211178737, -6.163307355639495,
+               -7.372177255047777, -8.488486734019723, -9.535449052433547, -10.527660396957408,
+               -11.475056633480246, -12.384788371845747)
+    biprime = (-2.2944396826141227, -4.073155089071816, -5.5123957296635835, -6.781294445990291,
+               -7.940178689168584, -9.019583358794248, -10.037696334908555, -11.00646266771229,
+               -11.934261645014844, -12.82725830917722)
+    #return (; ai, bi, aiprime, biprime) # Not compatible with Julia 1.0
+    return (ai=ai, bi=bi, aiprime=aiprime, biprime=biprime) # Compatible with Julia 1.0
+end
+
+"""
+    bessel_zero_largenu_asymptotic(nu::Real, m::Integer, kind::Integer, deriv::Bool)
+
+Compute an asymptotic approximation for a zero of the J or Y Bessel function (or their derivative) suitable for large order.
+
+## Input Arguments
+- `nu`: The (nonnegative) order of the Bessel function.
+- `m`: A small positive integer enumerating which zero is sought. Can not exceed `length(airy_zeros().ai)`.
+  The asymptotic formulas used herein weaken as `m` increases.
+- `kind`: Kind of Bessel function whose zero is sought. Either 1 for the J Bessel function or 2 for the Y Bessel function.
+- `deriv`: If false, returns the zero of the function.  If true, returns the zero of the function's derivative.
+
+## Return Value
+`z`: A `Float64` approximation of the desired zero.
+
+## Reference
+- https://dlmf.nist.gov/10.21.vii       
+"""
+function bessel_zero_largenu_asymptotic(nu::Real, m::Integer, kind::Integer, deriv::Bool)
+    abfactor = inv(cbrt(2.0))
+
+    if kind == 1
+        alpha = -abfactor * (deriv ? airy_zeros().aiprime[m] : airy_zeros().ai[m])
+    elseif kind == 2
+        alpha = -abfactor * (deriv ? airy_zeros().biprime[m] : airy_zeros().bi[m])
+    else
+        throw(ArgumentError("kind must be 1 or 2 but is $kind"))
+    end
+
+    alphap2 = alpha * alpha
+    alphap3 = alpha * alphap2
+    alphap4 = alpha * alphap3
+    alphap5 = alpha * alphap4
+
+    alpha0 = 1.0
+    alpha1 = alpha
+
+    if deriv
+        alpha2 = 3 * alphap2 / 10 - inv(10 * alpha)
+        alpha3 = -(alphap3 / 350 + inv(25) + inv(200 * alphap3))
+        alpha4 = -479 * alphap4 / 630_000 + 509 * alpha / 31_500 + inv(1500 * alphap2) - inv(2_000 * alphap5)
+        alpha5 = 0.0
+    else
+        alpha2 = 3 * alphap2 / 10
+        alpha3 = -alphap3 / 350 + inv(70)
+        alpha4 = -(479 * alphap4 / 63_000 + alpha / 3150)
+        alpha5 = 20_231 * alphap5 / 8_085_000 - 551 * alphap2 / 161_700
+    end
+
+    x = inv(cbrt(nu)^2)
+    xpk = 1.0
+    zsum = lastterm = alpha0
+    for alphak in (alpha1, alpha2, alpha3, alpha4, alpha5)
+        xpk *= x
+        nextterm = alphak * xpk
+        abs(nextterm) > abs(lastterm) && break # Asymptotic series starting to diverge
+        zsum += nextterm
+        lastterm = nextterm
+    end
+
+    zsum *= nu
+    return zsum
+end
+
 end # module FunctionZeros

test/runtests.jl

@@ -1,4 +1,6 @@
 using FunctionZeros
+using FunctionZeros: bessel_zero_asymptotic, bessel_deriv_zero_asymptotic
+
 using Test
 using SpecialFunctions: besselj, bessely
 
@@ -120,20 +122,48 @@ end
 @testset "high nu" begin
     @test isapprox(besselj_zero(6, 1), 9.936109524217688)
     @test isapprox(besselj_zero(7, 2), 14.821268727013171)
+    @test isapprox(besselj_zero(50, 1), 57.116899160119196)
+    @test isapprox(besselj_zero(60, 1), 67.52878576502944)
+    @test isapprox(besselj_zero(60, 2), 73.50669452996178)
+    @test isapprox(bessely_zero(20, 1), 22.625159280072324)
+end
 
-    if Int == Int64
-        @test isapprox(besselj_zero(50, 1), 57.116899160119196)
-    else
-        @test_broken isapprox(besselj_zero(50, 1), 57.116899160119196)
+@testset "higher nu" begin
+    nu = 93
+    zs = (101.63574127054777, 108.39596476159262, 114.11931644545753, 119.31845259184884, 124.18623309985139, 128.82065510579764)
+    for n in eachindex(zs)
+        @test isapprox(besselj_zero(nu, n), zs[n])
     end
-    # FIXME: This gives the incorrect result
-    # @test_broken isapprox(besselj_zero(60, 1), xxx)
-    if Int == Int64
-        @test isapprox(besselj_zero(60, 2), 73.50669452996178)
-    else
-        @test_broken isapprox(besselj_zero(60, 2), 73.50669452996178)
+    zs = (97.27824281868035, 105.20857143856234, 111.34228417034764, 116.76901331276218, 121.78633625489611, 126.5283633227293)
+    for n in eachindex(zs)
+        @test isapprox(bessely_zero(nu, n), zs[n])
+    end
+    zs = (96.67901927598446, 105.11090346045167, 111.2934441467005, 116.73708027868332, 121.76274824610373, 126.50968709216339)
+    for n in eachindex(zs)
+        @test isapprox(besselj_deriv_zero(nu, n), zs[n])
+    end
+    zs = (101.4575922200051, 108.33036194055933, 114.08063131302099, 119.29130692486132, 124.16538547400332, 128.8037392041164)
+    for n in eachindex(zs)
+        @test isapprox(bessely_deriv_zero(nu, n), zs[n])
+    end
+
+    nu = 3000.5
+    zs = (3027.337764568362, 3047.5168782184123, 3064.097399434406)
+    for n in eachindex(zs)
+        @test isapprox(besselj_zero(nu, n), zs[n])
+    end
+    zs = (3013.9544634926597, 3038.086941162325, 3056.1069360002584)
+    for n in eachindex(zs)
+        @test isapprox(bessely_zero(nu, n), zs[n])
+    end
+    zs = (3012.1679250547336, 3037.8201731857043, 3055.981972098408)
+    for n in eachindex(zs)
+        @test isapprox(besselj_deriv_zero(nu, n), zs[n])
+    end
+    zs = (3026.831390317444, 3047.3437713977214, 3064.0011561225597)
+    for n in eachindex(zs)
+        @test isapprox(bessely_deriv_zero(nu, n), zs[n])
     end
-    @test isapprox(bessely_zero(20, 1), 22.625159280072324)
 end
 
 @testset "asymptotic" begin
@@ -377,4 +407,14 @@ let nu = min(FunctionZeros.nupre_max, 1)
     @test abs(z - zbig) < 5 * eps()
     @test abs(bessely(nu-1, zbig) - bessely(nu+1, zbig)) / 2 < 2 * eps(BigFloat)
 end
+
+@testset "Issue 27" begin
+ @test besselj_deriv_zero(37, 1) ≈ besselj_deriv_zero(37.0, 1) ≈ 39.71488992674072
+ @test bessel_deriv_zero_asymptotic(37, 2, 1) ≈ 46.17514728525214
+ @test bessel_deriv_zero_asymptotic(BigInt(37), 2, 1) ≈ big"46.17514728525213690979739144407799904081911499170621436441096165051529368770411"
+ @test bessel_zero_asymptotic(87, 2, 1) ≈ 102.08833844612316
+ @test bessel_zero_asymptotic(BigInt(87), 2, 1) ≈ big"102.0883384461231427738288799405520802571778132437833489160205247116488836660219"
+end
+
+
 end # let