diff --git a/src/bessel.jl b/src/bessel.jl index ea65c64a..15c32b3a 100644 --- a/src/bessel.jl +++ b/src/bessel.jl @@ -184,7 +184,7 @@ end # besselj0, besselj1, bessely0, bessely1 for jy in ("j","y"), nu in (0,1) - jynu = Expr(:quote, Symbol(jy,nu)) + jynu = Expr(:quote, Symbol(jy,nu)) jynuf = Expr(:quote, Symbol(jy,nu,"f")) bjynu = Symbol("bessel",jy,nu) if jy == "y" @@ -201,8 +201,10 @@ for jy in ("j","y"), nu in (0,1) end end @eval begin - $bjynu(x::Real) = $bjynu(float(x)) - $bjynu(x::Complex) = $(Symbol("bessel",jy))($nu,x) + $bjynu(x::Real ) = $bjynu(float(x)) + $bjynu(x::Complex{Float64}) = $(Symbol("bessel",jy))($nu,x) + $bjynu(x::Complex{Float32}) = Complex{Float32}($bjynu(Complex{Float64}(x))) + $bjynu(x::Complex{Float16}) = Complex{Float16}($bjynu(Complex{Float64}(x))) end end @@ -389,8 +391,9 @@ function besselj(nu::Float64, z::Complex{Float64}) end end -besselj(nu::Cint, x::Float64) = ccall((:jn, libm), Float64, (Cint, Float64), nu, x) -besselj(nu::Cint, x::Float32) = ccall((:jnf, libm), Float32, (Cint, Float32), nu, x) +besselj(nu::Cint, x::Float64) = ccall((:jn, libm), Float64, (Cint, Float64), nu, x) +besselj(nu::Cint, x::Float32) = ccall((:jnf, libm), Float32, (Cint, Float32), nu, x) +besselj(nu::Cint, x::Float16)::Float16 = besselj(nu, Float32(x)) function besseljx(nu::Float64, z::Complex{Float64}) @@ -421,6 +424,7 @@ function bessely(nu::Cint, x::Float32) end ccall((:ynf, libm), Float32, (Cint, Float32), nu, x) end +bessely(nu::Cint, x::Float16) = Float16(bessely(nu, Float32(x))) function bessely(nu::Float64, z::Complex{Float64}) if nu < 0 @@ -430,14 +434,6 @@ function bessely(nu::Float64, z::Complex{Float64}) end end -function besselyx(nu::Float64, z::Complex{Float64}) - if nu < 0 - return _bessely(-nu,z,Int32(2))*cospi(nu) - _besselj(-nu,z,Int32(2))*sinpi(nu) - else - return _bessely(nu,z,Int32(2)) - end -end - """ besseli(nu, x) @@ -574,13 +570,21 @@ function besselyx(nu::Real, x::AbstractFloat) if x < 0 throw(DomainError(x, "`x` must be nonnegative.")) end - real(besselyx(float(nu), complex(x))) + convert(typeof(x), besselyx(float(nu), complex(x))) +end + +function besselyx(nu::Float64, z::Complex{Float64}) + if nu < 0 + return _bessely(-nu,z,Int32(2))*cospi(nu) - _besselj(-nu,z,Int32(2))*sinpi(nu) + else + return _bessely(nu,z,Int32(2)) + end end for f in ("i", "ix", "j", "jx", "k", "kx", "y", "yx") bfn = Symbol("bessel", f) @eval begin - $bfn(nu::Real, x::Real) = $bfn(nu, float(x)) + $bfn(nu::Real, x::Real) = typeof(float(x))($bfn(float(nu), float(x))) function $bfn(nu::Real, z::Complex) Tf = promote_type(float(typeof(nu)),float(typeof(real(z)))) $bfn(Tf(nu), Complex{Tf}(z)) diff --git a/test/bessel.jl b/test/bessel.jl index 5e853a8c..c3a3af48 100644 --- a/test/bessel.jl +++ b/test/bessel.jl @@ -35,22 +35,58 @@ end @testset "bessel functions" begin - bessel_funcs = [(bessely0, bessely1, bessely), (besselj0, besselj1, besselj)] - @testset "$z, $o" for (z, o, f) in bessel_funcs - @test z(Float32(2.0)) ≈ z(Float64(2.0)) - @test o(Float32(2.0)) ≈ o(Float64(2.0)) - @test z(Float16(2.0)) ≈ z(Float64(2.0)) - @test o(Float16(2.0)) ≈ o(Float64(2.0)) - @test z(2) ≈ z(2.0) - @test o(2) ≈ o(2.0) - @test z(2.0 + im) ≈ f(0, 2.0 + im) - @test o(2.0 + im) ≈ f(1, 2.0 + im) + @testset "bessel{j,y}{0,1}: check return value wrt bessel{j,y}" begin + # fixme can Symbol/eval combinations be simplified??? + for jy in ("j", "y") + bjy = Symbol("bessel", jy) + bjy_ = @eval begin $bjy end + for F in [Float16, Float32] + for nu in (0, 1) + bjynu = Symbol("bessel", jy, nu) + bjynu_ = @eval begin $bjynu end + + @test bjynu_( F(2.0) ) ≈ bjynu_(Float64(2.0) ) + @test bjynu_( 2 ) ≈ bjynu_( 2.0 ) + @test bjynu_( 2.0 ) ≈ bjy_(nu, 2.0 ) + @test bjynu_( 2.0+im ) ≈ bjy_(nu, 2.0+im) + @test bjynu_(Complex{F}(2.0+im)) ≈ bjynu_( 2.0+im) + end + end + end end - @testset "besselj error throwing" begin - @test_throws MethodError besselj(1.2,big(1.0)) - @test_throws MethodError besselj(1,complex(big(1.0))) - @test_throws MethodError besseljx(1,big(1.0)) - @test_throws MethodError besseljx(1,complex(big(1.0))) + + @testset "bessel{j,y}{,0,1}: correct return type" begin + @testset "type stability: $f" for f in [bessely0, bessely1, besselj0, besselj1] + for F in [Float16, Float32, Float64] + @test F == Base.return_types(f, Tuple{ F })[] + @test Complex{F} == Base.return_types(f, Tuple{Complex{F}})[] + end + @test BigFloat == Base.return_types(f, Tuple{BigFloat})[] + end + @testset "type stability: $f" for f in [besselj, bessely] + for F in [Float16, Float32, Float64] + for I in [Int16, Int32, Int64] + @test F == typeof(f(I(2), F( 2))) + @test Complex{promote_type(float(I),F)} == typeof(f(I(2), Complex{F}(2))) + for F2 in [Float16, Float32, Float64] + @test F == typeof(f(F2(2), F( 2))) + @test Complex{promote_type(F,F2)} == typeof(f(F2(2), Complex{F}(2))) + end + end + end + end + end + + @testset "besselj: undefined argument types" begin + @test_throws MethodError besselj( 1.2 , big( 1.0 )) + @test_throws MethodError besselj( 1 , complex(big( 1.0 ))) + @test_throws MethodError besselj(big( 1.0), 3im ) + @test_throws DomainError besselj( 0.1 , -0.4 ) + @test_throws AmosException besselj( 20 , 1000im ) + + @test_throws MethodError besseljx(1, big( 1.0) ) + @test_throws MethodError besseljx(1, complex(big(1.0))) + end @testset "besselh" begin true_h133 = 0.30906272225525164362 - 0.53854161610503161800im @@ -96,49 +132,49 @@ end @test_throws MethodError besselix(1,complex(big(1.0))) end end - @testset "besselj" begin - @test besselj(0,0) == 1 - for i in [-5 -3 -1 1 3 5] - @test besselj(i,0) == 0 - @test besselj(i,Float32(0)) == 0 - @test besselj(i,Complex{Float32}(0)) == 0.0 + @testset "besselj: specific values and (domain) errors" begin + # besselj(nu, 0) == { 1 for nu == 0 + # { 0 else + for nu in [-5 -3 -1 0 1 3 5] + for I in [Int16, Int32, Int64] + @test besselj(nu, zero( I )) == (nu == 0 ? 1 : 0) + end + for F in [Float16, Float32, Float64] + @test besselj(nu, zero( F )) == (nu == 0 ? 1 : 0) + @test besselj(nu, zero(Complex{F})) == (nu == 0 ? 1 : 0) + end end - j33 = besselj(3,3.) - @test besselj(3,3) == j33 - @test besselj(-3,-3) == j33 - @test besselj(-3,3) == -j33 - @test besselj(3,-3) == -j33 - @test besselj(3,3f0) ≈ j33 - @test besselj(3,complex(3.)) ≈ j33 - @test besselj(3,complex(3f0)) ≈ j33 - @test besselj(3,complex(3)) ≈ j33 - - j43 = besselj(4,3.) - @test besselj(4,3) == j43 - @test besselj(-4,-3) == j43 - @test besselj(-4,3) == j43 - @test besselj(4,-3) == j43 - @test besselj(4,3f0) ≈ j43 - @test besselj(4,complex(3.)) ≈ j43 - @test besselj(4,complex(3f0)) ≈ j43 - @test besselj(4,complex(3)) ≈ j43 - - @test j33 ≈ 0.30906272225525164362 - @test j43 ≈ 0.13203418392461221033 - @test besselj(0.1, complex(-0.4)) ≈ 0.820421842809028916 + 0.266571215948350899im - @test besselj(3.2, 1.3+0.6im) ≈ 0.01135309305831220201 + 0.03927719044393515275im - @test besselj(1, 3im) ≈ 3.953370217402609396im - @test besselj(1.0,3im) ≈ besselj(1,3im) - - true_jm3p1_3 = -0.45024252862270713882 - @test besselj(-3.1,3) ≈ true_jm3p1_3 - @test besselj(Float16(-3.1),Float16(3)) ≈ true_jm3p1_3 + j33 = besselj( 3, 3.0) + @test besselj( 3, 3) == j33 + @test besselj(-3, -3) == j33 + @test besselj(-3, 3) == -j33 + @test besselj( 3, -3) == -j33 + for F in [Float16, Float32, Float64] + @test besselj(3, F( 3)) ≈ j33 + @test besselj(3, Complex{F}(3)) ≈ j33 + end - @testset "Error throwing" begin - @test_throws DomainError besselj(0.1, -0.4) - @test_throws AmosException besselj(20,1000im) - @test_throws MethodError besselj(big(1.0),3im) + j43 = besselj( 4, 3.0) + @test besselj( 4, 3) == j43 + @test besselj(-4, -3) == j43 + @test besselj(-4, 3) == j43 + @test besselj( 4, -3) == j43 + for F in [Float16, Float32, Float64] + @test besselj(4, F( 3)) ≈ j43 + @test besselj(4, Complex{F}(3)) ≈ j43 + end + + @test besselj(1.0, 3im) ≈ besselj(1, 3im) + + # specific values + @test j33 ≈ 0.30906272225525164362 + @test j43 ≈ 0.13203418392461221033 + @test besselj(0.1, -0.4+ 0im) ≈ 0.820421842809028916 + 0.266571215948350899im + @test besselj(3.2, 1.3+0.6im) ≈ 0.01135309305831220201 + 0.03927719044393515275im + @test besselj(1, 3im) ≈ 3.953370217402609396im + for F in [Float16, Float32, Float64] + @test besselj(F(-3.1), F(3)) ≈ -0.45024252862270713882 end end @@ -164,29 +200,69 @@ end end @testset "bessely" begin - y33 = bessely(3,3.) - @test bessely(3,3) == y33 - @test bessely(3.,3.) == y33 - @test bessely(3,Float32(3.)) ≈ y33 - @test bessely(-3,3) ≈ -y33 - @test y33 ≈ -0.53854161610503161800 - @test bessely(3,complex(-3)) ≈ 0.53854161610503161800 - 0.61812544451050328724im + y33 = bessely(3, 3.0) + + @testset "same arguments, different data types" begin + @test bessely(3, 3 ) == y33 + @test bessely(3.0, 3.0) == y33 + @test bessely(3.0, 3 ) == y33 + for I in [Int16, Int32, Int64] + for F in [Float16, Float32, Float64] + @test bessely(I(3), F( 3)) ≈ y33 + @test bessely(I(3), Complex{F}(3)) ≈ y33 + end + end + end + + @testset "symmetry" begin + @test bessely(-3, 3) == -y33 + end + + @testset "specific values" begin + @test y33 ≈ -0.53854161610503161800 + @test bessely(3, complex(-3)) ≈ -y33 - 0.61812544451050328724im + end + @testset "Error throwing" begin @test_throws AmosException bessely(200.5,0.1) - @test_throws DomainError bessely(3,-3) - @test_throws DomainError bessely(0.4,-1.0) - @test_throws DomainError bessely(0.4,Float32(-1.0)) - @test_throws DomainError bessely(1,Float32(-1.0)) - @test_throws DomainError bessely(0.4,BigFloat(-1.0)) - @test_throws DomainError bessely(1,BigFloat(-1.0)) - @test_throws DomainError bessely(Cint(3),Float32(-3.)) - @test_throws DomainError bessely(Cint(3),Float64(-3.)) - - @test_throws MethodError bessely(1.2,big(1.0)) - @test_throws MethodError bessely(1,complex(big(1.0))) - @test_throws MethodError besselyx(1,big(1.0)) - @test_throws MethodError besselyx(1,complex(big(1.0))) + + @test_throws DomainError bessely( 3, -3 ) + @test_throws DomainError bessely(Cint(3), -3.0 ) + @test_throws DomainError bessely(Cint(3), Float32(-3.0)) + @test_throws DomainError bessely(0.4, -1.0 ) + @test_throws DomainError bessely(0.4, Float32( -1.0)) + @test_throws DomainError bessely(0.4, BigFloat(-1.0)) + @test_throws DomainError bessely(1, -1.0 ) + @test_throws DomainError bessely(1, Float32( -1.0)) + @test_throws DomainError bessely(1, BigFloat(-1.0)) + + @test_throws MethodError bessely(1.2, big(1.0) ) + @test_throws MethodError bessely(1, complex(big(1.0))) + @test_throws MethodError besselyx(1, big(1.0) ) + @test_throws MethodError besselyx(1.2, big(1.0) ) + @test_throws MethodError besselyx(1, complex(big(1.0))) + end + end + + @testset "besselyx" begin + @testset "return type" begin + for I in [Int16, Int32, Int64], I2 in [Int16, Int32, Int64] + @test Float64 == typeof(besselyx(I(2), I2(2))) + end + + for F in [Float16, Float32, Float64], F2 in [Float16, Float32, Float64] + @test F2 == typeof(besselyx(F(2), F2( 2))) + @test promote_type(F, Complex{F2}) == typeof(besselyx(F(2), Complex{F2}(2))) + end + + for F in [Float16, Float32, Float64], I in [Int16, Int32, Int64] + @test float(I) == typeof(besselyx(F(2), I( 2))) + @test F == typeof(besselyx(I(2), F( 2))) + @test promote_type(float(I), Complex{F}) == typeof(besselyx(I(2), Complex{F}(2))) + end end + + end @testset "sphericalbesselj" begin @@ -199,7 +275,7 @@ end @test sphericalbesselj(0, 0) == 1.0 @test sphericalbesselj(1, 0) == 0.0 @test sphericalbesselj(1, 0.01) ≈ 0.003333300000119047 - + @test_throws DomainError sphericalbesselj(1.25, -5.5) end @@ -211,7 +287,7 @@ end @test sphericalbessely(0, 1e-5) ≈ -99999.9999950000000 @test sphericalbessely(1, 1e-5) ≈ -1e10 - + @test_throws DomainError sphericalbessely(1.25, -5.5) @test_throws AmosException sphericalbessely(1, 0) end