Skip to content

Commit

Permalink
Correct besselj and besseljx for nu<0 and nu an Int.
Browse files Browse the repository at this point in the history
* besselj(nu, 0) and besseljx(nu, 0) should equal 0 for integer nu!=0.
  AmosException was thrown for z=complex(0).
* Tests corrected and more added to runtest.jl
* Changed `if sinpi(nu) == 0` to `if isinteger(nu)`
  • Loading branch information
Matthew Wingate authored and Matthew Wingate committed Jul 16, 2017
1 parent 8f7a958 commit 3132a80
Show file tree
Hide file tree
Showing 2 changed files with 29 additions and 9 deletions.
16 changes: 12 additions & 4 deletions src/bessel.jl
Original file line number Diff line number Diff line change
Expand Up @@ -272,7 +272,7 @@ end

function besseli(nu::Float64, z::Complex128)
if nu < 0
if sinpi(nu) == 0
if isinteger(nu)
return _besseli(-nu,z,Int32(1))
else
return _besseli(-nu,z,Int32(1)) - 2_besselk(-nu,z,Int32(1))*sinpi(nu)/pi
Expand All @@ -284,7 +284,7 @@ end

function besselix(nu::Float64, z::Complex128)
if nu < 0
if sinpi(nu) == 0
if isinteger(nu)
return _besseli(-nu,z,Int32(2))
else
return _besseli(-nu,z,Int32(2)) - 2_besselk(-nu,z,Int32(2))*exp(-abs(real(z))-z)*sinpi(nu)/pi
Expand All @@ -296,7 +296,11 @@ end

function besselj(nu::Float64, z::Complex128)
if nu < 0
return _besselj(-nu,z,Int32(1))*cospi(nu) + _bessely(-nu,z,Int32(1))*sinpi(nu)
if isinteger(nu)
return _besselj(-nu,z,Int32(1))*cospi(nu)
else
return _besselj(-nu,z,Int32(1))*cospi(nu) + _bessely(-nu,z,Int32(1))*sinpi(nu)
end
else
return _besselj(nu,z,Int32(1))
end
Expand All @@ -308,7 +312,11 @@ besselj(nu::Cint, x::Float32) = ccall((:jnf, libm), Float32, (Cint, Float32), nu

function besseljx(nu::Float64, z::Complex128)
if nu < 0
return _besselj(-nu,z,Int32(2))*cospi(nu) + _bessely(-nu,z,Int32(2))*sinpi(nu)
if isinteger(nu)
return _besselj(-nu,z,Int32(2))*cospi(nu)
else
return _besselj(-nu,z,Int32(2))*cospi(nu) + _bessely(-nu,z,Int32(2))*sinpi(nu)
end
else
return _besselj(nu,z,Int32(2))
end
Expand Down
22 changes: 17 additions & 5 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -161,6 +161,11 @@ end
@test SF.besseli(3,-3) -true_i33
@test SF.besseli(-3,-3) -true_i33
@test SF.besseli(Float32(-3),Complex64(-3,0)) -true_i33
for i in [-5 -3 -1 1 3 5]
@test SF.besseli(i,0) == 0.0
@test SF.besseli(i,Float32(0)) == 0
@test SF.besseli(i,Complex64(0)) == 0
end
@testset "Error throwing" begin
@test_throws SF.AmosException SF.besseli(1,1000)
@test_throws DomainError SF.besseli(0.4,-1.0)
Expand All @@ -172,11 +177,10 @@ end
end
@testset "besselj" begin
@test SF.besselj(0,0) == 1
for i = 1:5
for i in [-5 -3 -1 1 3 5]
@test SF.besselj(i,0) == 0
@test SF.besselj(-i,0) == 0
@test SF.besselj(-i,Float32(0)) == 0
@test SF.besselj(-i,Float32(0)) == 0
@test SF.besselj(i,Float32(0)) == 0
@test SF.besselj(i,Complex64(0)) == 0
end

j33 = SF.besselj(3,3.)
Expand Down Expand Up @@ -279,11 +283,19 @@ end
z == zero(z) || @test SF.besselyx(nu, z) SF.bessely(nu, z) * exp(-abs(imag(z)))
end
@test SF.besselkx(1, 0) == Inf
for i = [-5 -3 -1 1 3 5]
@test SF.besseljx(i,0) == 0
@test SF.besselix(i,0) == 0
@test SF.besseljx(i,Float32(0)) == 0
@test SF.besselix(i,Float32(0)) == 0
@test SF.besseljx(i,Complex64(0)) == 0
@test SF.besselix(i,Complex64(0)) == 0
end
@testset "Error throwing" begin
@test_throws SF.AmosException SF.hankelh1x(1, 0)
@test_throws SF.AmosException SF.hankelh2x(1, 0)
@test_throws SF.AmosException SF.besselix(-1.01, 0)
@test_throws SF.AmosException SF.besseljx(-1, 0)
@test_throws SF.AmosException SF.besseljx(-1.01, 0)
@test_throws SF.AmosException SF.besselyx(1, 0)
@test_throws DomainError SF.besselix(0.4,-1.0)
@test_throws DomainError SF.besseljx(0.4, -1.0)
Expand Down

0 comments on commit 3132a80

Please sign in to comment.