From dcef20185b6f1db7fe74486724ac6cdbce2d0b4c Mon Sep 17 00:00:00 2001 From: Matthew Wingate Date: Sat, 15 Jul 2017 17:00:22 +0100 Subject: [PATCH 1/3] Correct besseli(nu, 0) & besselix(nu, 0) with nu<0 and nu an Int. An AmosException is thrown for besseli(nu, 0) or besselix(nu, 0) due to the call to _besselk. However, for integer nu, besseli(nu, 0) = besseli(-nu,0) = 0 and similar for besselix. A test is added before calling besselk. --- src/bessel.jl | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/bessel.jl b/src/bessel.jl index 978fd28f..f368237c 100644 --- a/src/bessel.jl +++ b/src/bessel.jl @@ -272,7 +272,11 @@ end function besseli(nu::Float64, z::Complex128) if nu < 0 - return _besseli(-nu,z,Int32(1)) - 2_besselk(-nu,z,Int32(1))*sinpi(nu)/pi + if sinpi(nu) == 0 + return _besseli(-nu,z,Int32(1)) + else + return _besseli(-nu,z,Int32(1)) - 2_besselk(-nu,z,Int32(1))*sinpi(nu)/pi + end else return _besseli(nu,z,Int32(1)) end @@ -280,7 +284,11 @@ end function besselix(nu::Float64, z::Complex128) if nu < 0 - return _besseli(-nu,z,Int32(2)) - 2_besselk(-nu,z,Int32(2))*exp(-abs(real(z))-z)*sinpi(nu)/pi + if sinpi(nu) == 0 + 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 + end else return _besseli(nu,z,Int32(2)) end From 8f7a958eafc09311254f46c018e719f4594706d7 Mon Sep 17 00:00:00 2001 From: Matthew Wingate Date: Sat, 15 Jul 2017 18:33:10 +0100 Subject: [PATCH 2/3] Correct test for AmosException in besselix(nu, z) besselix(nu, z) is divergent for z=0 and negative nu, only if nu is not an integer. --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 747cf724..889a53af 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -282,7 +282,7 @@ 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, 0) + @test_throws SF.AmosException SF.besselix(-1.01, 0) @test_throws SF.AmosException SF.besseljx(-1, 0) @test_throws SF.AmosException SF.besselyx(1, 0) @test_throws DomainError SF.besselix(0.4,-1.0) From 3132a806f8d4b455f2a9da14598f3da9cd033886 Mon Sep 17 00:00:00 2001 From: Matthew Wingate Date: Sun, 16 Jul 2017 11:16:03 +0100 Subject: [PATCH 3/3] Correct besselj and besseljx for nu<0 and nu an Int. * 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)` --- src/bessel.jl | 16 ++++++++++++---- test/runtests.jl | 22 +++++++++++++++++----- 2 files changed, 29 insertions(+), 9 deletions(-) diff --git a/src/bessel.jl b/src/bessel.jl index f368237c..33f7d8ed 100644 --- a/src/bessel.jl +++ b/src/bessel.jl @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/test/runtests.jl b/test/runtests.jl index 889a53af..8b5795ea 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -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) @@ -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.) @@ -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)