From e9a1fa66b2c8925c62a37b8868957d1e90e5fbdd Mon Sep 17 00:00:00 2001 From: Peter Simon Date: Sun, 10 Aug 2025 11:10:44 -0700 Subject: [PATCH 1/5] precomputation table --- src/FunctionZeros.jl | 55 ++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 50 insertions(+), 5 deletions(-) diff --git a/src/FunctionZeros.jl b/src/FunctionZeros.jl index 2b40a81..cb952d7 100644 --- a/src/FunctionZeros.jl +++ b/src/FunctionZeros.jl @@ -4,6 +4,10 @@ import Roots export besselj_zero, bessely_zero +# Set max order and index of zeros to precompute +const nupre_max = 1 +const npre_max = 500 + besselj_zero_asymptotic(nu, n) = bessel_zero_asymptotic(nu, n, 1) @@ -42,25 +46,66 @@ end # Order 0 is 6 times slower and 50-100 times less accurate # than higher orders, with other parameters constant. """ - besselj_zero(nu, n; order=2) + _besselj_zero(nu, n; order=2) `n`th zero of the Bessel J function of order `nu`, -for `n` = `1,2,...`. +for `n` = `1,2,...`. `order` is passed to the function `Roots.fzero`. """ -besselj_zero(nu, n; order=2) = Roots.find_zero((x) -> SpecialFunctions.besselj(nu, x), +_besselj_zero(nu, n; order=2) = Roots.find_zero((x) -> SpecialFunctions.besselj(nu, x), bessel_zero_asymptotic(nu, n, 1); order=order) + +const jzero_pre = [_besselj_zero(nu, n; order=2) for n in 1:npre_max, nu in 0:nupre_max] + """ - bessely_zero(nu, n; order=2) + besselj_zero(nu, n; order=2) + +`n`th zero of the Bessel J function of order `nu`, +for `n` = `1,2,...`. + +`order` is passed to the function `Roots.fzero` for roots not precomputed in the lookup table. + +For greater speed, table lookup is used when `nu ∈ 0:nupre_max` and `n ∈ 1:npremax`. +""" +function besselj_zero(nu, n; order=2) + if nu in 0:nupre_max && n in 1:npre_max + return jzero_pre[Int(n), Int(nu) + 1] + else + return _besselj_zero(nu, n; order) + end +end + +""" + _bessely_zero(nu, n; order=2) `n`th zero of the Bessel Y function of order `nu`, for `n` = `1,2,...`. `order` is passed to the function `Roots.fzero`. """ -bessely_zero(nu, n; order=2) = Roots.find_zero((x) -> SpecialFunctions.bessely(nu, x), +_bessely_zero(nu, n; order=2) = Roots.find_zero((x) -> SpecialFunctions.bessely(nu, x), bessel_zero_asymptotic(nu, n, 2); order=order) +const yzero_pre = [_bessely_zero(nu, n; order=2) for n in 1:npre_max, nu in 0:nupre_max] + +""" + bessely_zero(nu, n; order=2) + +`n`th zero of the Bessel Y function of order `nu`, +for `n` = `1,2,...`. + +`order` is passed to the function `Roots.fzero` for roots not precomputed in the lookup table. + +For greater speed, table lookup is used when `nu ∈ 0:nupre_max` and `n ∈ 1:npremax`. +""" +function bessely_zero(nu, n; order=2) + if nu in 0:nupre_max && n in 1:npre_max + return yzero_pre[Int(n), Int(nu) + 1] + else + return _bessely_zero(nu, n; order) + end +end + end # module FunctionZeros From d46c6216d7a28a686151c8f33878e321be13ad58 Mon Sep 17 00:00:00 2001 From: Peter Simon Date: Sun, 10 Aug 2025 13:46:17 -0700 Subject: [PATCH 2/5] Add precomputed tables (Issue 7) and fix Issue #21 --- src/FunctionZeros.jl | 24 ++++++++++++++++-------- test/runtests.jl | 15 +++++++++++++++ 2 files changed, 31 insertions(+), 8 deletions(-) diff --git a/src/FunctionZeros.jl b/src/FunctionZeros.jl index cb952d7..210bff1 100644 --- a/src/FunctionZeros.jl +++ b/src/FunctionZeros.jl @@ -4,7 +4,7 @@ import Roots export besselj_zero, bessely_zero -# Set max order and index of zeros to precompute +# Set max order and index of zeros to precompute and tabulate const nupre_max = 1 const npre_max = 500 @@ -56,8 +56,8 @@ for `n` = `1,2,...`. _besselj_zero(nu, n; order=2) = Roots.find_zero((x) -> SpecialFunctions.besselj(nu, x), bessel_zero_asymptotic(nu, n, 1); order=order) - -const jzero_pre = [_besselj_zero(nu, n; order=2) for n in 1:npre_max, nu in 0:nupre_max] +# tabulation of selected besselj_zero values +const jzero_pre = [_besselj_zero(nu, n; order=2) for nu in 0:nupre_max, n in 1:npre_max] """ besselj_zero(nu, n; order=2) @@ -71,7 +71,7 @@ For greater speed, table lookup is used when `nu ∈ 0:nupre_max` and `n ∈ 1:n """ function besselj_zero(nu, n; order=2) if nu in 0:nupre_max && n in 1:npre_max - return jzero_pre[Int(n), Int(nu) + 1] + return jzero_pre[Int(nu) + 1, Int(n)] else return _besselj_zero(nu, n; order) end @@ -85,10 +85,18 @@ for `n` = `1,2,...`. `order` is passed to the function `Roots.fzero`. """ -_bessely_zero(nu, n; order=2) = Roots.find_zero((x) -> SpecialFunctions.bessely(nu, x), - bessel_zero_asymptotic(nu, n, 2); order=order) +function _bessely_zero(nu, n; order=2) + if isone(n) && abs(nu) < 0.1587 # See Issue 21 + return Roots.find_zero((x) -> SpecialFunctions.bessely(nu, x), + 0.5 * (nu + besselj_zero(nu, n)); order) + else + return Roots.find_zero((x) -> SpecialFunctions.bessely(nu, x), + bessel_zero_asymptotic(nu, n, 2); order) + end +end -const yzero_pre = [_bessely_zero(nu, n; order=2) for n in 1:npre_max, nu in 0:nupre_max] +# tabulation of selected bessely_zero values +const yzero_pre = [_bessely_zero(nu, n; order=2) for nu in 0:nupre_max, n in 1:npre_max] """ bessely_zero(nu, n; order=2) @@ -102,7 +110,7 @@ For greater speed, table lookup is used when `nu ∈ 0:nupre_max` and `n ∈ 1:n """ function bessely_zero(nu, n; order=2) if nu in 0:nupre_max && n in 1:npre_max - return yzero_pre[Int(n), Int(nu) + 1] + return yzero_pre[Int(nu) + 1, Int(n)] else return _bessely_zero(nu, n; order) end diff --git a/test/runtests.jl b/test/runtests.jl index e2aa123..980df60 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -138,3 +138,18 @@ end @testset "asymptotic" begin FunctionZeros.bessel_zero_asymptotic(4, 2, 1) == FunctionZeros.besselj_zero_asymptotic(4, 2) end + +@testset "Issue 21" begin + @test bessely_zero(-0.1586, 1) ≈ 0.6559631635143002 + @test bessely_zero(0, 1) ≈ 0.8935769662791675 + @test bessely_zero(0.1586, 1) ≈ 1.117167411163268 +end + +@testset "precomputed" begin + for nu in 0:FunctionZeros.nupre_max + for n in union(1:27:FunctionZeros.npre_max, FunctionZeros.npre_max) + @test isapprox(besselj_zero(nu, n), FunctionZeros._besselj_zero(nu, n)) + @test isapprox(bessely_zero(nu, n), FunctionZeros._bessely_zero(nu, n)) + end + end +end From 35090224502d06ebf0604d52ef7bcc81559812b2 Mon Sep 17 00:00:00 2001 From: Peter Simon Date: Sun, 10 Aug 2025 13:56:04 -0700 Subject: [PATCH 3/5] Use syntax for keyword arguments consistent with older Julia versions --- src/FunctionZeros.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/FunctionZeros.jl b/src/FunctionZeros.jl index 210bff1..b152e23 100644 --- a/src/FunctionZeros.jl +++ b/src/FunctionZeros.jl @@ -73,7 +73,7 @@ function besselj_zero(nu, n; order=2) if nu in 0:nupre_max && n in 1:npre_max return jzero_pre[Int(nu) + 1, Int(n)] else - return _besselj_zero(nu, n; order) + return _besselj_zero(nu, n; order=order) end end @@ -88,10 +88,10 @@ for `n` = `1,2,...`. function _bessely_zero(nu, n; order=2) if isone(n) && abs(nu) < 0.1587 # See Issue 21 return Roots.find_zero((x) -> SpecialFunctions.bessely(nu, x), - 0.5 * (nu + besselj_zero(nu, n)); order) + 0.5 * (nu + besselj_zero(nu, n)); order=order) else return Roots.find_zero((x) -> SpecialFunctions.bessely(nu, x), - bessel_zero_asymptotic(nu, n, 2); order) + bessel_zero_asymptotic(nu, n, 2); order=order) end end From f6fcd070639bf0ec770fde92a44059a736e20126 Mon Sep 17 00:00:00 2001 From: Peter Simon Date: Sun, 10 Aug 2025 13:56:04 -0700 Subject: [PATCH 4/5] Use syntax for keyword arguments consistent with older Julia versions --- src/FunctionZeros.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/FunctionZeros.jl b/src/FunctionZeros.jl index 210bff1..a6c8d0a 100644 --- a/src/FunctionZeros.jl +++ b/src/FunctionZeros.jl @@ -73,7 +73,7 @@ function besselj_zero(nu, n; order=2) if nu in 0:nupre_max && n in 1:npre_max return jzero_pre[Int(nu) + 1, Int(n)] else - return _besselj_zero(nu, n; order) + return _besselj_zero(nu, n; order=order) end end @@ -88,10 +88,10 @@ for `n` = `1,2,...`. function _bessely_zero(nu, n; order=2) if isone(n) && abs(nu) < 0.1587 # See Issue 21 return Roots.find_zero((x) -> SpecialFunctions.bessely(nu, x), - 0.5 * (nu + besselj_zero(nu, n)); order) + 0.5 * (nu + besselj_zero(nu, n)); order=order) else return Roots.find_zero((x) -> SpecialFunctions.bessely(nu, x), - bessel_zero_asymptotic(nu, n, 2); order) + bessel_zero_asymptotic(nu, n, 2); order=order) end end @@ -112,7 +112,7 @@ function bessely_zero(nu, n; order=2) if nu in 0:nupre_max && n in 1:npre_max return yzero_pre[Int(nu) + 1, Int(n)] else - return _bessely_zero(nu, n; order) + return _bessely_zero(nu, n; order=order) end end From a599865ad37dbb27eed292dc8901b984481fc81f Mon Sep 17 00:00:00 2001 From: Peter Simon Date: Mon, 11 Aug 2025 19:56:10 -0700 Subject: [PATCH 5/5] Use precomputed tables only for Float64 outputs --- src/FunctionZeros.jl | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/FunctionZeros.jl b/src/FunctionZeros.jl index a6c8d0a..5f8f86d 100644 --- a/src/FunctionZeros.jl +++ b/src/FunctionZeros.jl @@ -56,7 +56,7 @@ for `n` = `1,2,...`. _besselj_zero(nu, n; order=2) = Roots.find_zero((x) -> SpecialFunctions.besselj(nu, x), bessel_zero_asymptotic(nu, n, 1); order=order) -# tabulation of selected besselj_zero values +# Float64 tabulation of selected besselj_zero values const jzero_pre = [_besselj_zero(nu, n; order=2) for nu in 0:nupre_max, n in 1:npre_max] """ @@ -67,9 +67,9 @@ for `n` = `1,2,...`. `order` is passed to the function `Roots.fzero` for roots not precomputed in the lookup table. -For greater speed, table lookup is used when `nu ∈ 0:nupre_max` and `n ∈ 1:npremax`. +For greater speed in Float64, table lookup is used when `nu ∈ 0:nupre_max` and `n ∈ 1:npremax`. """ -function besselj_zero(nu, n; order=2) +function besselj_zero(nu::Union{Int,Float64}, n::Integer; order=2) if nu in 0:nupre_max && n in 1:npre_max return jzero_pre[Int(nu) + 1, Int(n)] else @@ -77,6 +77,8 @@ function besselj_zero(nu, n; order=2) end end +besselj_zero(nu::Real, n::Integer; order=2) = _besselj_zero(nu::Real, n::Integer; order=order) + """ _bessely_zero(nu, n; order=2) @@ -88,14 +90,14 @@ for `n` = `1,2,...`. function _bessely_zero(nu, n; order=2) if isone(n) && abs(nu) < 0.1587 # See Issue 21 return Roots.find_zero((x) -> SpecialFunctions.bessely(nu, x), - 0.5 * (nu + besselj_zero(nu, n)); order=order) + (nu + besselj_zero(nu, n)) / 2; order=order) else return Roots.find_zero((x) -> SpecialFunctions.bessely(nu, x), bessel_zero_asymptotic(nu, n, 2); order=order) end end -# tabulation of selected bessely_zero values +# tabulation of selected bessely_zero values in Float64 const yzero_pre = [_bessely_zero(nu, n; order=2) for nu in 0:nupre_max, n in 1:npre_max] """ @@ -108,7 +110,7 @@ for `n` = `1,2,...`. For greater speed, table lookup is used when `nu ∈ 0:nupre_max` and `n ∈ 1:npremax`. """ -function bessely_zero(nu, n; order=2) +function bessely_zero(nu::Union{Int,Float64}, n::Integer; order=2) if nu in 0:nupre_max && n in 1:npre_max return yzero_pre[Int(nu) + 1, Int(n)] else @@ -116,4 +118,6 @@ function bessely_zero(nu, n; order=2) end end +bessely_zero(nu::Real, n::Integer; order=2) = _bessely_zero(nu, n; order=order) + end # module FunctionZeros