Skip to content

Type stability of Bessel functions #234

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 12 commits into
base: master
Choose a base branch
from
Draft

Type stability of Bessel functions #234

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
554f4c6
bessel j0,j1,y0,y1: test for type stability
Jun 19, 2020
d25b7b4
bessel j0,j1,y0,y1: type stable for Complex{FloatXY}
Jun 19, 2020
d107fa0
Merge remote-tracking branch 'upstream/master' into pc/bessel
Jun 19, 2020
74204c7
added: besselj(Int32, Float16)
Jun 21, 2020
077e5e2
type stability tests for besselj(Int32, FloatXY)
Jun 21, 2020
4b7a513
bessel j0,j1,y0,y1: testing Complex{Float16/32}
Jun 21, 2020
2f5eef4
besselj,y: tests
Jun 26, 2020
d105427
bessel{j,y} tests
Sep 22, 2020
e5d4d32
bessel{j,y} tests
Sep 22, 2020
df24280
code alignment
Sep 22, 2020
b131e8c
bessely tests
Sep 24, 2020
5e23d1c
besselyx: first tests for type stability
Oct 7, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 19 additions & 15 deletions src/bessel.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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

Expand Down Expand Up @@ -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})
Expand Down Expand Up @@ -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
Expand All @@ -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)

Expand Down Expand Up @@ -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))
Expand Down
230 changes: 153 additions & 77 deletions test/bessel.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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
Expand All @@ -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

Expand All @@ -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
Expand Down