Merge pull request #71 from JuliaMath/fixes

Fix error in recurrence when length is 2

Author: heltonmc

NEWS.md

@@ -10,8 +10,14 @@ For bug fixes, performance enhancements, or fixes to unexported functions we wil
 
 # Unreleased
 
+# Version 0.2.8
+
 ### Added
 - Add support for complex numbers in `besseli0`, `besseli1`, `besselj0`, `besselj1` ([PR #68](https://github.com/JuliaMath/Bessels.jl/pull/68))
+- Add separate documentation page with API ([PR #69](https://github.com/JuliaMath/Bessels.jl/pull/69)). This currently fails to build see ([Issue #70](https://github.com/JuliaMath/Bessels.jl/issues/70))
+
+### Fixed
+- Fixed wrong return when iterating over a range of `nu` values when `length(nu) == 2` ([PR #71](https://github.com/JuliaMath/Bessels.jl/pull/71))
 
 # Version 0.2.7
 
@@ -34,7 +40,7 @@ For bug fixes, performance enhancements, or fixes to unexported functions we wil
 # Version 0.2.5
 
 ### Fixed
-- Fix bug for very large inputs (x>1e16) in besselj0 and friends routines. A particularly thank you to @jwscook for reporting the bug in ([Issue #53](https://github.com/JuliaMath/Bessels.jl/pull/56)) and providing a very detailed analysis.
+- Fix bug for very large inputs (x>1e16) in besselj0 and friends routines. A particularly thank you to @jwscook for reporting the bug in ([Issue #56](https://github.com/JuliaMath/Bessels.jl/issues/56)) and providing a very detailed analysis.
 
 
 # Version 0.2.4

Project.toml

@@ -1,7 +1,7 @@
 name = "Bessels"
 uuid = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
 authors = ["Michael Helton <[email protected]> and contributors"]
-version = "0.2.7"
+version = "0.2.8"
 
 [compat]
 julia = "1.6"

README.md

@@ -11,7 +11,7 @@ Numerical routines for computing Bessel, Airy, and Hankel functions for real arg
 
 The goal of the library is to provide high quality numerical implementations of Bessel functions with high accuracy without comprimising on computational time. In general, we try to match (and often exceed) the accuracy of other open source routines such as those provided by [SpecialFunctions.jl](https://github.com/JuliaMath/SpecialFunctions.jl). There are instances where we don't quite match that desired accuracy (within a digit or two) but in general will provide implementations that are 5-10x faster (see [benchmarks](https://github.com/JuliaMath/Bessels.jl#benchmarks)).
 
-The library currently supports Bessel functions, modified Bessel functions, Hankel functions, spherical Bessel functions, and Airy functions of the first and second kind for positive real arguments and integer and noninteger orders. Negative arguments are also supported only if the return value is real. We plan to support complex arguments in the future. An unexported gamma function is also provided.
+The library currently supports Bessel functions, modified Bessel functions, Hankel functions, spherical Bessel functions, and Airy functions of the first and second kind for positive real arguments and integer and noninteger orders. Negative arguments are also supported only if the return value is real. [Limited support](https://github.com/JuliaMath/Bessels.jl#complex-numbers) is provided for complex arguments. An unexported gamma function is also provided.
 
 # Quick start
 
@@ -152,6 +152,10 @@ a = zeros(10)
 out = Bessels.besselj!(a, 1:10, 1.0)
 ```
 
+### Complex numbers
+
+Support for complex numbers is only provided for the Airy functions (`airyai`, `airyaiprime`, `airybi`, `airybiprime`) and the Bessel functions of the first kind with orders 0 and 1 (`besselj0`, `besselj1`, `besseli0`, `besseli1`).
+
 ### Support for negative arguments
 
 Support is provided for negative arguments and orders only if the return value is real. A domain error will be thrown if the return value is complex. See https://github.com/heltonmc/Bessels.jl/issues/30 for more details.

src/besseli.jl

@@ -90,7 +90,7 @@ function besseli0x(x::T) where T <: Union{Float32, Float64}
 end
 
 """
-    besseli1(x::T) where T <: Union{Float32, Float64}
+    besseli1(x::T) where T <: Union{Float32, Float64, ComplexF32, ComplexF64}
 
 Modified Bessel function of the first kind of order one, ``I_1(x)``.
 
@@ -489,7 +489,7 @@ function _besseli!(out::DenseVector{T}, nu::AbstractRange, x::T) where T
         k -= 1
         k < 1 && break
     end
-    if k > 1
+    if k >= 1
         out[k] = _besseli(nu[k], x)
         tmp = @view out[begin:k+1]
         besselk_down_recurrence!(tmp, x, nu[begin:k+1])

src/besselj.jl

@@ -102,7 +102,7 @@ function _besselj0(x::Float32)
 end
 
 """
-    besselj1(x::T) where T <: Union{Float32, Float64}
+    besselj1(x::T) where T <: Union{Float32, Float64, ComplexF32, ComplexF64}
 
 Bessel function of the first kind of order one, ``J_1(x)``.
 
@@ -358,7 +358,7 @@ function _besselj!(out::DenseVector{T}, nu::AbstractVector, x::T) where T <: Uni
             k -= 1
             k < 1 && break
         end
-        if k > 1
+        if k >= 1
             out[k] = _besselj(nu[k], x)
             tmp = @view out[begin:k+1]
             besselj_down_recurrence!(tmp, x, nu[begin:k+1])

src/besselk.jl

@@ -297,7 +297,7 @@ function _besselk!(out::DenseVector{T}, nu::AbstractRange, x::T) where T
         k += 1
         k == len && break
     end
-    if k < len
+    if k <= len
         out[k] = _besselk(nu[k], x)
         tmp = @view out[k-1:end]
         besselk_up_recurrence!(tmp, x, nu[k-1:end])

test/besseli_test.jl

@@ -118,6 +118,7 @@ end
 # test nu_range
 @test besseli(0:250, 2.0) ≈ SpecialFunctions.besseli.(0:250, 2.0) rtol=1e-13
 @test besseli(0.5:1:10.5, 2.0) ≈ SpecialFunctions.besseli.(0.5:1:10.5, 2.0) rtol=1e-13
+@test besseli(3:4, 1.2) ≈ SpecialFunctions.besseli.(3:4, 1.2)
 @test Bessels.besseli!(zeros(Float64, 10), 1:10, 1.0) ≈ besseli(1:10, 1.0)
 
 ### need to fix method ambiguities for other functions ###### 

test/besselj_test.jl

@@ -134,6 +134,7 @@ end
 @test besselj(0:95, 100.0) ≈ SpecialFunctions.besselj.(0:95, 100.0) rtol=1e-11
 @test besselj(0.5:1:150.5, 2.0) ≈ SpecialFunctions.besselj.(0.5:1:150.5, 2.0) rtol=1e-11
 @test besselj(0.5:1:10.5, 40.0) ≈ SpecialFunctions.besselj.(0.5:1:10.5, 40.0) rtol=1e-11
+@test besselj(3:4, 1.2) ≈ SpecialFunctions.besselj.(3:4, 1.2)
 @test Bessels.besselj!(zeros(Float64, 10), 1:10, 1.0) ≈ besselj(1:10, 1.0)
 
 # test Float16 and Float32

test/besselk_test.jl

@@ -130,6 +130,7 @@ end
 @test besselk(0:50, 2.0) ≈ SpecialFunctions.besselk.(0:50, 2.0) rtol=1e-13
 @test besselk(0.5:1:10.5, 12.0) ≈ SpecialFunctions.besselk.(0.5:1:10.5, 12.0) rtol=1e-13
 @test besselk(1:700, 800.0) ≈ SpecialFunctions.besselk.(1:700, 800.0)
+@test besselk(3:4, 1.2) ≈ SpecialFunctions.besselk.(3:4, 1.2)
 @test Bessels.besselk!(zeros(Float64, 10), 1:10, 1.0) ≈ besselk(1:10, 1.0)
 
 # test Float16