From 494816f24451817aee6124dc8d3980e0dd1f5351 Mon Sep 17 00:00:00 2001 From: Peter Simon Date: Fri, 17 Oct 2025 15:04:11 -0700 Subject: [PATCH 1/8] Fix Issue #27 --- src/FunctionZeros.jl | 4 ++++ test/runtests.jl | 12 ++++++++++++ 2 files changed, 16 insertions(+) diff --git a/src/FunctionZeros.jl b/src/FunctionZeros.jl index 4112e84..d523d20 100644 --- a/src/FunctionZeros.jl +++ b/src/FunctionZeros.jl @@ -76,6 +76,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 @@ -214,6 +216,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) diff --git a/test/runtests.jl b/test/runtests.jl index 4f14889..e5979dd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,6 @@ using FunctionZeros +using FunctionZeros: bessel_zero_asymptotic, bessel_deriv_zero_asymptotic + using Test using SpecialFunctions: besselj, bessely @@ -377,4 +379,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) ≈ 47.22868085729076 + @test bessel_deriv_zero_asymptotic(BigInt(37), 2, 1) ≈ big"47.22868085729076149030198319523818371811907166937005677029525603864958297422779" + @test bessel_zero_asymptotic(87, 2, 1) ≈ 105.69324719238529 + @test bessel_zero_asymptotic(BigInt(87), 2, 1) ≈ big"105.693247192385300157734443161801368795079219923515831608871687185822839565836" +end + + end # let From 4410a5e9f4fee091fd1dac761b6784344945efcc Mon Sep 17 00:00:00 2001 From: Peter Simon Date: Fri, 17 Oct 2025 15:37:55 -0700 Subject: [PATCH 2/8] Convert @test_broken to @test --- test/runtests.jl | 16 +++------------- 1 file changed, 3 insertions(+), 13 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index e5979dd..6141844 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -122,19 +122,9 @@ end @testset "high nu" begin @test isapprox(besselj_zero(6, 1), 9.936109524217688) @test isapprox(besselj_zero(7, 2), 14.821268727013171) - - if Int == Int64 - @test isapprox(besselj_zero(50, 1), 57.116899160119196) - else - @test_broken isapprox(besselj_zero(50, 1), 57.116899160119196) - 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) - end + @test isapprox(besselj_zero(50, 1), 57.116899160119196) + @test isapprox(besselj_zero(60, 1), 73.5066945299618) + @test isapprox(besselj_zero(60, 2), 73.50669452996178) @test isapprox(bessely_zero(20, 1), 22.625159280072324) end From fa52c859645d60fd663315a04c979fb8699f37b8 Mon Sep 17 00:00:00 2001 From: Peter Simon Date: Sun, 19 Oct 2025 15:47:25 -0700 Subject: [PATCH 3/8] Add Limitations section to README.md --- README.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/README.md b/README.md index 498b1d9..e1268e1 100644 --- a/README.md +++ b/README.md @@ -20,6 +20,11 @@ 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 three zeros (`n` ∈ {1,2,3}) of any of the four functions treated here are always found +correctly for any choice of `nu`. Similarly, when `nu ≤ 93`, any choice for `n` will produce a +correct result. However, for `nu > 93` and `n > 3` some of the zeros may be skipped or found in the wrong order. + ### Exported Functions #### besselj_zero(nu, n) From 890a138db748986c4ae52384cd79babe8dcdf771 Mon Sep 17 00:00:00 2001 From: Peter Simon Date: Sun, 19 Oct 2025 15:48:12 -0700 Subject: [PATCH 4/8] Add additional tests for large nu values --- test/runtests.jl | 48 +++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 43 insertions(+), 5 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 6141844..61066e2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -123,11 +123,49 @@ end @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), 73.5066945299618) + @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 +@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 + 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 +end + @testset "asymptotic" begin FunctionZeros.bessel_zero_asymptotic(4, 2, 1) == FunctionZeros.besselj_zero_asymptotic(4, 2) FunctionZeros.bessel_zero_asymptotic(4, 2, 2) == FunctionZeros.bessely_zero_asymptotic(4, 2) @@ -372,10 +410,10 @@ 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) ≈ 47.22868085729076 - @test bessel_deriv_zero_asymptotic(BigInt(37), 2, 1) ≈ big"47.22868085729076149030198319523818371811907166937005677029525603864958297422779" - @test bessel_zero_asymptotic(87, 2, 1) ≈ 105.69324719238529 - @test bessel_zero_asymptotic(BigInt(87), 2, 1) ≈ big"105.693247192385300157734443161801368795079219923515831608871687185822839565836" + @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 From 2c26ca78589cdb23a29a6726045764a75c9fb2ef Mon Sep 17 00:00:00 2001 From: Peter Simon Date: Sun, 19 Oct 2025 16:01:38 -0700 Subject: [PATCH 5/8] Implement asymptotic formulas valid for large nu and small n --- src/FunctionZeros.jl | 92 +++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 91 insertions(+), 1 deletion(-) diff --git a/src/FunctionZeros.jl b/src/FunctionZeros.jl index d523d20..00efc28 100644 --- a/src/FunctionZeros.jl +++ b/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 > 25 && n ≤ length(airy_zeros().ai) + return bessel_zero_largenu_asymptotic(nu, n, kind, false) + end if kind == 1 beta = MathConstants.pi * (n + nu / 2 - 1//4) else # kind == 2 @@ -191,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 > 25 && n ≤ length(airy_zeros().ai) + 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 @@ -297,4 +304,87 @@ 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, say 2 or 3. +""" +@inline function airy_zeros() + ai = (-2.338107410459767, -4.087949444130973, -5.520559828095556) + bi = (-1.173713222709127, -3.2710933028363516, -4.8307378416620095) + aiprime = (-1.0187929716474717, -3.248197582179841, -4.820099211178737) + biprime = (-2.2944396826141227, -4.073155089071816, -5.5123957296635835) + return (; ai, bi, aiprime, biprime) +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 From 2ddbd98e34b8d6eeef1145892c4d1ddb4aaec43c Mon Sep 17 00:00:00 2001 From: Peter Simon Date: Sun, 19 Oct 2025 16:22:19 -0700 Subject: [PATCH 6/8] Fix named tuple constructor for Julia 1.0 compatibility --- src/FunctionZeros.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/FunctionZeros.jl b/src/FunctionZeros.jl index 00efc28..a0cfe5e 100644 --- a/src/FunctionZeros.jl +++ b/src/FunctionZeros.jl @@ -320,7 +320,8 @@ Return the first few negative zeros of the functions `airyai`, `airybi`, `airyai bi = (-1.173713222709127, -3.2710933028363516, -4.8307378416620095) aiprime = (-1.0187929716474717, -3.248197582179841, -4.820099211178737) biprime = (-2.2944396826141227, -4.073155089071816, -5.5123957296635835) - return (; ai, bi, aiprime, biprime) + #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 """ From c0d7ba354412a8726f0f9339ba41bb846be04dd1 Mon Sep 17 00:00:00 2001 From: Peter Simon Date: Mon, 20 Oct 2025 10:58:45 -0700 Subject: [PATCH 7/8] Increase size of table for airy function zeros to 10 --- README.md | 7 ++++--- src/FunctionZeros.jl | 18 +++++++++++++----- 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index e1268e1..a59caeb 100644 --- a/README.md +++ b/README.md @@ -21,9 +21,10 @@ of the first values of `1, 2, 3, ...`. See the individual function docstrings f extents of the lookup tables. ### Limitations -The first three zeros (`n` ∈ {1,2,3}) of any of the four functions treated here are always found -correctly for any choice of `nu`. Similarly, when `nu ≤ 93`, any choice for `n` will produce a -correct result. However, for `nu > 93` and `n > 3` some of the zeros may be skipped or found in the wrong order. +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 diff --git a/src/FunctionZeros.jl b/src/FunctionZeros.jl index a0cfe5e..86d822b 100644 --- a/src/FunctionZeros.jl +++ b/src/FunctionZeros.jl @@ -313,13 +313,21 @@ Return the first few negative zeros of the functions `airyai`, `airybi`, `airyai - `(; 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, say 2 or 3. + is a small integer, currently 10. """ @inline function airy_zeros() - ai = (-2.338107410459767, -4.087949444130973, -5.520559828095556) - bi = (-1.173713222709127, -3.2710933028363516, -4.8307378416620095) - aiprime = (-1.0187929716474717, -3.248197582179841, -4.820099211178737) - biprime = (-2.2944396826141227, -4.073155089071816, -5.5123957296635835) + 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 From d585e2da509e4b19eca98a3f2bcf110489aa3d82 Mon Sep 17 00:00:00 2001 From: Peter Simon Date: Mon, 20 Oct 2025 15:12:03 -0700 Subject: [PATCH 8/8] Adjust test for which asymptotic formula to use. --- src/FunctionZeros.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/FunctionZeros.jl b/src/FunctionZeros.jl index 86d822b..0140426 100644 --- a/src/FunctionZeros.jl +++ b/src/FunctionZeros.jl @@ -54,7 +54,7 @@ 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 > 25 && n ≤ length(airy_zeros().ai) + 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 @@ -195,7 +195,7 @@ of order `nu`. `kind == 1 (2)` for Bessel function of the first (second) kind, J """ function bessel_deriv_zero_asymptotic(nu_in::Real, n::Integer, kind=1) nu = abs(nu_in) - if nu > 25 && n ≤ length(airy_zeros().ai) + 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